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 Helper routines for the Floquet-Bloch band-structure calculation (floquet_main).
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_blacs_env, ONLY: cp_blacs_env_create,&
17 : cp_blacs_env_type
18 : USE cp_cfm_types, ONLY: cp_cfm_get_submatrix,&
19 : cp_cfm_set_submatrix,&
20 : cp_cfm_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
22 : USE cp_files, ONLY: close_file,&
23 : open_file
24 : USE floquet_types, ONLY: floquet_env_type
25 : USE input_constants, ONLY: small_cell_full_kp
26 : USE kinds, ONLY: default_string_length,&
27 : dp,&
28 : int_8
29 : USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices,&
30 : rs_to_kp
31 : USE kpoint_methods, ONLY: kpoint_init_cell_index
32 : USE kpoint_types, ONLY: get_kpoint_info,&
33 : kpoint_create,&
34 : kpoint_release,&
35 : kpoint_type
36 : USE machine, ONLY: m_memory_details
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_comm_split_type_shared,&
44 : mp_comm_type,&
45 : mp_para_env_type
46 : USE physcon, ONLY: a_bohr,&
47 : evolt,&
48 : kelvin
49 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
50 : USE qs_environment_types, ONLY: get_qs_env,&
51 : qs_environment_type
52 : USE qs_mo_types, ONLY: get_mo_set,&
53 : mo_set_type
54 : USE qs_moments, ONLY: qs_moment_kpoints_deep
55 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
56 : USE util, ONLY: sort
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'floquet_utils'
64 :
65 : ! Number of full Floquet-Hamiltonian copies (each n_f_size^2 complex(dp), 16 bytes per
66 : ! element) that the subgroup must be able to hold, distributed across its ranks.
67 : ! We size for 8 to leave headroom.
68 : INTEGER, PARAMETER, PRIVATE :: n_fm_work_copies = 8
69 :
70 : ! Public subroutines
71 : PUBLIC :: build_floquet_matrix, &
72 : compute_e_k_de_dk_dipole, &
73 : make_floquet_subgroups, &
74 : distribute_floquet_kp_data, &
75 : floquet_sector_weights, &
76 : check_floquet_convergence, &
77 : calculate_floquet_observables, &
78 : write_floquet_header, &
79 : write_floquet_results
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief ...
85 : !> \param qs_env ...
86 : !> \param xkp ...
87 : !> \param e_k ...
88 : !> \param de_dk ...
89 : !> \param do_parallel the option to distribute the results (e_k/de_dk) across MPI ranks.
90 : !> When .TRUE. k-point ikp is computed and stored only on rank MOD(ikp-1,num_pe)
91 : !> Default .FALSE. -> every rank computes/stores all k-points (replicated).
92 : ! **************************************************************************************************
93 2 : SUBROUTINE calculate_epsilon_derivative(qs_env, xkp, e_k, de_dk, do_parallel)
94 : TYPE(qs_environment_type), POINTER :: qs_env
95 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
96 : INTENT(IN) :: xkp
97 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
98 : INTENT(OUT), OPTIONAL :: e_k
99 : REAL(KIND=dp), ALLOCATABLE, &
100 : DIMENSION(:, :, :, :), INTENT(OUT), OPTIONAL :: de_dk
101 : LOGICAL, INTENT(IN), OPTIONAL :: do_parallel
102 :
103 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_epsilon_derivative'
104 :
105 2 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: C_dH_C, C_dS_C, C_k, dH_dk_i, dS_dk_i, &
106 2 : H_k, S_k
107 : INTEGER :: handle, i_dir, ikp, ispin, mepos, n, &
108 : n_img_all, n_spin, nao, nkp, num_copy, &
109 : num_pe
110 2 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all
111 2 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
112 : LOGICAL :: my_do_parallel, present_dedk, present_ek
113 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvals
114 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: H_rs, S_rs
115 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
116 : TYPE(cell_type), POINTER :: cell
117 2 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
118 : TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
119 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
120 : TYPE(mp_para_env_type), POINTER :: para_env
121 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
122 2 : POINTER :: sab_all
123 :
124 2 : CALL timeset(routineN, handle)
125 :
126 2 : present_ek = PRESENT(e_k) ! calculate band energies, ε_k for all kpoints xkp
127 2 : present_dedk = PRESENT(de_dk) ! calculate derivative, ∇_k ε_k of band energies for all kpoints xkp
128 :
129 2 : IF (.NOT. (present_ek .OR. present_dedk)) CPABORT("Subroutine needs either e_k or de_dk")
130 :
131 2 : my_do_parallel = .FALSE.
132 2 : IF (PRESENT(do_parallel)) my_do_parallel = do_parallel
133 :
134 : CALL get_qs_env(qs_env, &
135 : matrix_ks_kp=matrix_ks_kp, &
136 : matrix_s_kp=matrix_s_kp, &
137 : sab_all=sab_all, &
138 : cell=cell, &
139 : kpoints=kpoints_scf, &
140 : para_env=para_env, &
141 2 : mos=mos)
142 :
143 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
144 2 : CALL get_cell(cell=cell, h=hmat)
145 :
146 2 : n_spin = SIZE(matrix_ks_kp, 1)
147 2 : nkp = SIZE(xkp, 2)
148 :
149 : ! Distribution of k-points across ranks (mirrors qs_moment_kpoints_deep)
150 : ! do_parallel = .TRUE. => ikp stored in mepos==MOD(ikp-1,num_pe)
151 : ! do_parallel = .FALSE. => every rank computes all k-points (replicated)
152 2 : mepos = 0
153 2 : num_pe = 1
154 2 : num_copy = nkp
155 2 : IF (my_do_parallel) THEN
156 2 : mepos = para_env%mepos
157 2 : num_pe = para_env%num_pe
158 2 : num_copy = CEILING(REAL(nkp)/num_pe)
159 : END IF
160 :
161 : ! create kpoint environment kpoints_all which contains all neighbor cells R
162 : ! without considering any lattice symmetry
163 2 : NULLIFY (kpoints_all)
164 2 : CALL kpoint_create(kpoints_all)
165 2 : CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_all)
166 2 : CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index_all, index_to_cell=index_to_cell_all)
167 :
168 20 : ALLOCATE (S_rs(1, nao, nao, n_img_all), H_rs(n_spin, nao, nao, n_img_all), source=0.0_dp)
169 :
170 : ! Convert real-space dbcsr matrices into arrays
171 2 : CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, S_rs, cell_to_index_all)
172 2 : CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, H_rs, cell_to_index_all)
173 :
174 12 : IF (present_dedk) ALLOCATE (de_dk(n_spin, num_copy, 3, nao), source=0.0_dp)
175 10 : IF (present_ek) ALLOCATE (e_k(n_spin, num_copy, nao), source=0.0_dp)
176 :
177 : !$OMP PARALLEL DEFAULT(NONE) &
178 : !$OMP PRIVATE(ikp, ispin, S_k, H_k, eigenvals, C_k, dS_dk_i, dH_dk_i, C_dS_C, C_dH_C) &
179 : !$OMP SHARED(nao, n_spin, de_dk, e_k, present_ek, present_dedk, mepos, num_pe, &
180 2 : !$OMP nkp, xkp, S_rs, H_rs, index_to_cell_all, hmat)
181 : IF (present_dedk) ALLOCATE (dS_dk_i(nao, nao), C_dS_C(nao, nao), &
182 : dH_dk_i(nao, nao), C_dH_C(nao, nao), source=z_zero)
183 : ALLOCATE (C_k(nao, nao), S_k(nao, nao), H_k(nao, nao), source=z_zero)
184 : ALLOCATE (eigenvals(nao), source=0.0_dp)
185 : !$OMP DO COLLAPSE(2)
186 : DO ispin = 1, n_spin
187 : DO ikp = 1, nkp
188 : IF (MOD(ikp - 1, num_pe) /= mepos) CYCLE
189 :
190 : ! S^R -> S(k), H^R -> H(k)
191 : S_k = 0
192 : H_k = 0
193 : CALL rs_to_kp(S_rs(1, :, :, :), S_k, index_to_cell_all, xkp(:, ikp))
194 : CALL rs_to_kp(H_rs(ispin, :, :, :), H_k, index_to_cell_all, xkp(:, ikp))
195 :
196 : ! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
197 : CALL geeig_right(H_k, S_k, eigenvals, C_k)
198 : IF (present_ek) e_k(ispin, CEILING(REAL(ikp)/num_pe), :) = eigenvals(:)
199 :
200 : IF (present_dedk) THEN
201 : ! Evaluate the derivatives using
202 : ! ∇ ε_k = C^H(k) ∇ H_k C(k) - ε_k C^H(k) ∇ S_k C(k)
203 : DO i_dir = 1, 3
204 : CALL rs_to_kp(S_rs(1, :, :, :), dS_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
205 : CALL rs_to_kp(H_rs(ispin, :, :, :), dH_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
206 :
207 : CALL gemm_square(C_k, 'C', dS_dk_i, 'N', C_k, 'N', C_dS_C)
208 : CALL gemm_square(C_k, 'C', dH_dk_i, 'N', C_k, 'N', C_dH_C)
209 :
210 : DO n = 1, nao
211 : de_dk(ispin, CEILING(REAL(ikp)/num_pe), i_dir, n) = &
212 : DBLE(C_dH_C(n, n)) - DBLE(eigenvals(n)*C_dS_C(n, n))
213 : END DO
214 : END DO
215 : END IF
216 : END DO
217 : END DO
218 : !$OMP END DO
219 : IF (present_dedk) DEALLOCATE (dS_dk_i, C_dS_C, dH_dk_i, C_dH_C)
220 : DEALLOCATE (S_k, H_k, C_k, eigenvals)
221 : !$OMP END PARALLEL
222 2 : DEALLOCATE (S_rs, H_rs)
223 2 : CALL kpoint_release(kpoints_all)
224 :
225 2 : CALL timestop(handle)
226 :
227 8 : END SUBROUTINE calculate_epsilon_derivative
228 :
229 : ! **************************************************************************************************
230 : !> \brief Momentum matrix elements p_nm(k) for one k-point and one spin.
231 : !> \param e_k_kp_spin ...
232 : !> \param de_dk_kp_spin ...
233 : !> \param dipole_kp_spin ...
234 : !> \param momentum ...
235 : ! **************************************************************************************************
236 2 : SUBROUTINE build_momentum_matrix(e_k_kp_spin, de_dk_kp_spin, dipole_kp_spin, momentum)
237 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e_k_kp_spin
238 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: de_dk_kp_spin
239 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dipole_kp_spin
240 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: momentum
241 :
242 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_momentum_matrix'
243 :
244 : INTEGER :: handle, i, i_dir, j, nao
245 :
246 2 : CALL timeset(routineN, handle)
247 :
248 2 : nao = SIZE(e_k_kp_spin)
249 :
250 : ! We calculate momentum matrix elements p_nm = <ψ_n|-iħ∇_r|ψ_m>
251 : ! p_nm = <u_n|(ħk - iħ∇_r)|u_m>
252 : ! p_nm = i d_nm (ε_n - ε_m) + ∇_k ε_n δ_nm
253 :
254 : !$OMP PARALLEL DEFAULT(NONE) &
255 : !$OMP PRIVATE(i_dir, i, j) &
256 2 : !$OMP SHARED(nao, momentum, dipole_kp_spin, e_k_kp_spin, de_dk_kp_spin)
257 : !$OMP DO COLLAPSE(3)
258 : DO i_dir = 1, 3
259 : DO i = 1, nao
260 : DO j = 1, nao
261 : IF (j == i) THEN
262 : momentum(i_dir, i, j) = de_dk_kp_spin(i_dir, i)
263 : ELSE
264 : momentum(i_dir, i, j) = gaussi*(e_k_kp_spin(i) - e_k_kp_spin(j))* &
265 : dipole_kp_spin(i_dir, i, j)
266 : END IF
267 : END DO
268 : END DO
269 : END DO
270 : !$OMP END DO
271 : !$OMP END PARALLEL
272 :
273 2 : CALL timestop(handle)
274 2 : END SUBROUTINE build_momentum_matrix
275 :
276 : ! **************************************************************************************************
277 : !> \brief Off-diagonal coupling block of H_F(k) for one k-point and spin, from precomputed data.
278 : !> \param bs_env ...
279 : !> \param e_k_kp_spin ...
280 : !> \param de_dk_kp_spin ...
281 : !> \param dipole_kp_spin ...
282 : !> \param off_diag_m ...
283 : ! **************************************************************************************************
284 2 : SUBROUTINE build_off_diagonal_matrix(bs_env, e_k_kp_spin, de_dk_kp_spin, dipole_kp_spin, off_diag_m)
285 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
286 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e_k_kp_spin
287 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: de_dk_kp_spin
288 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dipole_kp_spin
289 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: off_diag_m
290 :
291 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_off_diagonal_matrix'
292 :
293 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: momentum
294 : COMPLEX(KIND=dp), DIMENSION(3) :: efactor
295 : INTEGER :: handle, i, i_dir, j, nao
296 : REAL(KIND=dp) :: amplitude, omega
297 : REAL(KIND=dp), DIMENSION(3) :: e_vec, phi, polarisation
298 :
299 2 : CALL timeset(routineN, handle)
300 :
301 : ! Builds the matrix that occupies the off diagonal blocks in the Floquet Matrix H_F
302 :
303 2 : nao = SIZE(e_k_kp_spin)
304 :
305 8 : polarisation(:) = bs_env%floquet_polarisation(:)
306 8 : IF (SQRT(SUM(polarisation**2)) < EPSILON(0.0_dp)) THEN
307 0 : CPABORT("Invalid (too small) polarisation vector specified for POLARISATION")
308 : END IF
309 :
310 2 : amplitude = bs_env%floquet_amplitude
311 8 : e_vec(:) = amplitude*polarisation
312 8 : phi(:) = pi*bs_env%floquet_phi(:)
313 2 : omega = bs_env%floquet_omega
314 :
315 8 : ALLOCATE (momentum(3, nao, nao), source=z_zero)
316 2 : CALL build_momentum_matrix(e_k_kp_spin, de_dk_kp_spin, dipole_kp_spin, momentum)
317 :
318 : ! E_factor(α) = (i E(α) exp(i·φ(α)))/(2ω), where α = x,y,z
319 8 : DO i_dir = 1, 3
320 8 : efactor(i_dir) = gaussi*e_vec(i_dir)*CMPLX(COS(phi(i_dir)), SIN(phi(i_dir)), KIND=dp)/(2*omega)
321 : END DO
322 :
323 : ! off_diag_uv = Σ_α [p_uv^α(k) · E_factor(α)], where α = x,y,z
324 146 : off_diag_m(:, :) = z_zero
325 8 : DO i_dir = 1, 3
326 56 : DO i = 1, nao
327 438 : DO j = 1, nao
328 : off_diag_m(i, j) = off_diag_m(i, j) + &
329 432 : momentum(i_dir, i, j)*efactor(i_dir)
330 : END DO
331 : END DO
332 : END DO
333 2 : DEALLOCATE (momentum)
334 :
335 2 : CALL timestop(handle)
336 :
337 2 : END SUBROUTINE build_off_diagonal_matrix
338 :
339 : ! **************************************************************************************************
340 : !> \brief Diagonal block of H_F(k) for Floquet sector f_index.
341 : !> \param bs_env ...
342 : !> \param e_k_kp_spin ...
343 : !> \param f_index Floquet sector index
344 : !> \param diag_e ...
345 : ! **************************************************************************************************
346 202 : SUBROUTINE build_diagonal_matrix(bs_env, e_k_kp_spin, f_index, diag_e)
347 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
348 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e_k_kp_spin
349 : INTEGER, INTENT(IN) :: f_index
350 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: diag_e
351 :
352 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_diagonal_matrix'
353 :
354 : INTEGER :: handle, i, nao
355 : REAL(KIND=dp) :: omega
356 :
357 202 : CALL timeset(routineN, handle)
358 : ! Builds the diagonal block of the H_F for Floquet sector index f_index
359 : ! Equilibrium band energies shifted by f_index·ħΩ, i.e., diag_e(n,n) = ε_{nk} + f_index·ħΩ.
360 202 : nao = SIZE(e_k_kp_spin)
361 202 : omega = bs_env%floquet_omega
362 14746 : diag_e(:, :) = z_zero
363 1818 : DO i = 1, nao
364 1818 : diag_e(i, i) = e_k_kp_spin(i) + f_index*omega
365 : END DO
366 :
367 202 : CALL timestop(handle)
368 :
369 202 : END SUBROUTINE build_diagonal_matrix
370 :
371 : ! **************************************************************************************************
372 : !> \brief Builds the Floquet-Bloch Hamiltonian H_F(k) for a single k-point and spin channel from
373 : !> the band data currently held in floquet_env.
374 : !> \param bs_env ...
375 : !> \param floquet_env ...
376 : !> \param floquet_matrix ...
377 : ! **************************************************************************************************
378 2 : SUBROUTINE build_floquet_matrix(bs_env, floquet_env, floquet_matrix)
379 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
380 : TYPE(floquet_env_type), INTENT(IN) :: floquet_env
381 : TYPE(cp_cfm_type), INTENT(IN) :: 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, nao
388 :
389 2 : CALL timeset(routineN, handle)
390 :
391 2 : nao = floquet_env%nao
392 2 : max_f_index = floquet_env%max_f_index
393 2 : n_fbands = 1 + 2*max_f_index
394 :
395 : ! Creates the floquet matrix H_F by placing the diagonal and off diagonal blocks
396 8 : ALLOCATE (diag_e(nao, nao), source=z_zero)
397 6 : ALLOCATE (off_diag_m(nao, nao), source=z_zero)
398 6 : ALLOCATE (conj_off_diag_m(nao, nao), source=z_zero)
399 :
400 : CALL build_off_diagonal_matrix(bs_env, floquet_env%e_k_kp_spin, floquet_env%de_dk_kp_spin, &
401 2 : floquet_env%dipole_kp_spin, off_diag_m)
402 146 : conj_off_diag_m(:, :) = CONJG(TRANSPOSE(off_diag_m(:, :)))
403 :
404 654482 : floquet_matrix%local_data(:, :) = z_zero
405 204 : DO i = 1, n_fbands
406 202 : i_f = 1 + (i - 1)*nao
407 202 : f_index = i - max_f_index - 1
408 202 : CALL build_diagonal_matrix(bs_env, floquet_env%e_k_kp_spin, f_index, diag_e)
409 202 : CALL cp_cfm_set_submatrix(floquet_matrix, diag_e, i_f, i_f)
410 204 : IF (i > 1) THEN
411 200 : CALL cp_cfm_set_submatrix(floquet_matrix, off_diag_m, i_f, i_f - nao)
412 200 : CALL cp_cfm_set_submatrix(floquet_matrix, conj_off_diag_m, i_f - nao, i_f)
413 : END IF
414 : END DO
415 2 : DEALLOCATE (diag_e, off_diag_m, conj_off_diag_m)
416 :
417 2 : CALL timestop(handle)
418 2 : END SUBROUTINE build_floquet_matrix
419 :
420 : ! **************************************************************************************************
421 : !> \brief Make this k-point/spin's band data available to every rank of the owning subgroup.
422 : !> \param para_env the global parallel environment
423 : !> \param para_env_sub the subgroup parallel environment
424 : !> \param ispin the spin channel
425 : !> \param ikp the DOS k-point index
426 : !> \param e_k band energies for all DOS k-points, distributed
427 : !> \param de_dk band-energy k-derivatives, distributed
428 : !> \param dipole dipole matrix elements, distributed
429 : !> \param floquet_env ...
430 : ! **************************************************************************************************
431 2 : SUBROUTINE distribute_floquet_kp_data(para_env, para_env_sub, ispin, ikp, e_k, de_dk, dipole, &
432 : floquet_env)
433 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
434 : INTEGER, INTENT(IN) :: ispin, ikp
435 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: e_k
436 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: de_dk
437 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
438 : INTENT(IN) :: dipole
439 : TYPE(floquet_env_type), INTENT(INOUT) :: floquet_env
440 :
441 : CHARACTER(LEN=*), PARAMETER :: routineN = 'distribute_floquet_kp_data'
442 :
443 : INTEGER :: handle, loc_idx, local_src, owner
444 :
445 2 : CALL timeset(routineN, handle)
446 :
447 : ! Identify the owner's rank within the subgroup (all subgroup ranks agree on local_src).
448 2 : owner = MOD(ikp - 1, para_env%num_pe)
449 :
450 2 : local_src = -1
451 2 : IF (para_env%mepos == owner) local_src = para_env_sub%mepos
452 2 : CALL para_env_sub%max(local_src)
453 :
454 : ! The owner copies its single-spin slice of this k-point's distributed band data into the
455 : ! broadcast buffers, then broadcasts them to the whole subgroup.
456 2 : IF (para_env%mepos == owner) THEN
457 1 : loc_idx = CEILING(REAL(ikp)/para_env%num_pe)
458 9 : floquet_env%e_k_kp_spin(:) = e_k(ispin, loc_idx, :)
459 33 : floquet_env%de_dk_kp_spin(:, :) = de_dk(ispin, loc_idx, :, :)
460 265 : floquet_env%dipole_kp_spin(:, :, :) = dipole(ispin, loc_idx, :, :, :)
461 : END IF
462 2 : CALL para_env_sub%bcast(floquet_env%e_k_kp_spin, local_src)
463 2 : CALL para_env_sub%bcast(floquet_env%de_dk_kp_spin, local_src)
464 2 : CALL para_env_sub%bcast(floquet_env%dipole_kp_spin, local_src)
465 :
466 2 : CALL timestop(handle)
467 :
468 2 : END SUBROUTINE distribute_floquet_kp_data
469 :
470 : ! **************************************************************************************************
471 : !> \brief Precompute, on the global communicator, the band quantities needed to assemble the
472 : !> Floquet-Bloch Hamiltonian for all DOS k-points (all spins), distributed across ranks.
473 : !> \param qs_env ...
474 : !> \param bs_env ...
475 : !> \param e_k ...
476 : !> \param de_dk ...
477 : !> \param dipole ...
478 : ! **************************************************************************************************
479 2 : SUBROUTINE compute_e_k_de_dk_dipole(qs_env, bs_env, e_k, de_dk, dipole)
480 : TYPE(qs_environment_type), POINTER :: qs_env
481 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
482 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
483 : INTENT(OUT) :: e_k
484 : REAL(KIND=dp), ALLOCATABLE, &
485 : DIMENSION(:, :, :, :), INTENT(OUT) :: de_dk
486 : COMPLEX(KIND=dp), ALLOCATABLE, &
487 : DIMENSION(:, :, :, :, :), INTENT(OUT) :: dipole
488 :
489 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_e_k_de_dk_dipole'
490 :
491 : INTEGER :: handle, ikp, nkp_only_bs, nkp_start
492 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp_all
493 :
494 2 : CALL timeset(routineN, handle)
495 :
496 : ! The band-structure k-points are the last nkp_only_bs entries of kpoints_DOS%xkp, stored
497 : ! after the nkp_only_DOS DOS-only k-points.
498 2 : nkp_start = bs_env%nkp_only_DOS
499 2 : nkp_only_bs = bs_env%nkp_only_bs
500 :
501 : ! Collect the coordinates of all band-structure k-points (global index nkp_start + local).
502 6 : ALLOCATE (xkp_all(3, nkp_only_bs))
503 4 : DO ikp = 1, nkp_only_bs
504 10 : xkp_all(:, ikp) = bs_env%kpoints_DOS%xkp(1:3, nkp_start + ikp)
505 : END DO
506 :
507 : ! Calculate and distribute the results across ranks
508 : ! One k-point per rank, round-robin by global rank
509 : ! So the results for ikp are stored in mepos==MOD(ikp-1,num_pe)
510 2 : CALL calculate_epsilon_derivative(qs_env, xkp_all, e_k=e_k, de_dk=de_dk, do_parallel=.TRUE.)
511 2 : CALL qs_moment_kpoints_deep(qs_env, xkp_all, dipole, do_parallel=.TRUE.)
512 :
513 2 : DEALLOCATE (xkp_all)
514 :
515 2 : CALL timestop(handle)
516 2 : END SUBROUTINE compute_e_k_de_dk_dipole
517 :
518 : ! **************************************************************************************************
519 : !> \brief Finds the number of MPI ranks that share a physical node, determined by splitting the
520 : !> given communicator into node-local communicators via MPI's shared-memory split type.
521 : !> \param para_env ...
522 : !> \return the number of ranks on the calling rank's node
523 : ! **************************************************************************************************
524 2 : FUNCTION find_ranks_per_node(para_env) RESULT(ranks_per_node)
525 : TYPE(mp_para_env_type), POINTER :: para_env
526 : INTEGER :: ranks_per_node
527 :
528 : TYPE(mp_comm_type) :: node_comm
529 :
530 2 : CALL node_comm%from_split_type(para_env, mp_comm_split_type_shared, key=para_env%mepos)
531 2 : ranks_per_node = MAX(1, node_comm%num_pe)
532 2 : CALL node_comm%free()
533 :
534 2 : END FUNCTION find_ranks_per_node
535 :
536 : ! **************************************************************************************************
537 : !> \brief Choose the optimal number of MPI ranks per subgroup for Floquet calculations.
538 : !> \param qs_env ...
539 : !> \param bs_env ...
540 : !> \return the chosen number of ranks per subgroup (1 .. num_pe)
541 : ! **************************************************************************************************
542 4 : FUNCTION floquet_determine_subgroup_size(qs_env, bs_env) RESULT(group_size)
543 : TYPE(qs_environment_type), POINTER :: qs_env
544 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
545 : INTEGER :: group_size
546 :
547 : CHARACTER(LEN=*), PARAMETER :: routineN = 'floquet_determine_subgroup_size'
548 :
549 : INTEGER :: g_load, g_mem, handle, n_f_size, nao, &
550 : nkp_only_bs, ranks_per_node, unit_nr
551 : INTEGER(KIND=int_8) :: Buffers, Cached, MemFree, MemLikelyFree, &
552 : MemTotal, needed_bytes, Slab, &
553 : SReclaimable, usable_per_rank
554 : REAL(KIND=dp) :: mem_fill_fraction
555 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
556 : TYPE(mp_para_env_type), POINTER :: para_env
557 :
558 2 : CALL timeset(routineN, handle)
559 :
560 2 : CALL get_qs_env(qs_env, para_env=para_env, mos=mos)
561 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
562 :
563 2 : unit_nr = bs_env%unit_nr
564 2 : n_f_size = nao*(1 + 2*bs_env%max_floquet_index)
565 2 : nkp_only_bs = bs_env%nkp_only_bs
566 2 : mem_fill_fraction = bs_env%floquet_mem_fill_fraction
567 :
568 : ! We choose the subgroup size to be the larger of two floors:
569 : ! 1. Memory floor G_mem: the subgroup must jointly hold the Floquet working set,
570 : ! G_mem = ceil(needed_bytes/usable per-rank memory), where the usable memory is
571 : ! the mem_fill_fraction of the per-rank FREE memory.
572 : ! 2. K-point floor G_load = floor(num_pe / nkp_only_bs): the LARGEST G that still
573 : ! yields at least nkp_only_bs subgroups (floor(num_pe/G) >= nkp_only_bs). When ranks
574 : ! outnumber k-points this hands every k-point its own (larger, hence faster)
575 : ! subgroup instead of leaving spare ranks idle.
576 :
577 : ! usable_per_rank = a fraction of the per-rank FREE memory, not including SCF and other data
578 2 : CALL m_memory_details(MemTotal, MemFree, Buffers, Cached, Slab, SReclaimable, MemLikelyFree)
579 2 : ranks_per_node = find_ranks_per_node(para_env)
580 2 : usable_per_rank = INT(mem_fill_fraction*REAL(MemFree, dp), int_8)/INT(ranks_per_node, int_8)
581 :
582 : ! Total memory the subgroup must hold: n_fm_work_copies copies of the Floquet
583 : ! Hamiltonian, each n_f_size^2 complex(dp) numbers at 16 bytes each.
584 2 : needed_bytes = 16_int_8*INT(n_fm_work_copies, int_8)*INT(n_f_size, int_8)**2
585 :
586 2 : IF (usable_per_rank <= 0_int_8) THEN
587 : ! Memory info unavailable (e.g. no /proc found) use one subgroup spanning all ranks
588 0 : g_mem = para_env%num_pe
589 : ELSE
590 : ! Memory floor: smallest G with needed_bytes/G <= usable_per_rank
591 : g_mem = INT(MIN((needed_bytes + usable_per_rank - 1_int_8)/usable_per_rank, &
592 2 : INT(para_env%num_pe, int_8)))
593 : END IF
594 :
595 2 : IF (g_mem > para_env%num_pe) CPWARN("Total Memory Likely Insufficient, process may be killed")
596 :
597 : ! K-point floor: floor(num_pe / nkp_only_bs) = largest G that gives at least nkp_only_bs subgroups
598 2 : g_load = para_env%num_pe/MAX(nkp_only_bs, 1)
599 :
600 2 : group_size = MAX(g_mem, g_load)
601 2 : group_size = MIN(group_size, para_env%num_pe)
602 2 : group_size = MAX(group_size, 1)
603 2 : CALL para_env%max(group_size)
604 :
605 2 : IF (unit_nr > 0) THEN
606 1 : WRITE (unit_nr, '(T2,A)') ""
607 1 : WRITE (unit_nr, '(T2,A,T66,I15)') "FLOQUET MEMORY | Detected MPI ranks per node:", ranks_per_node
608 1 : WRITE (unit_nr, '(T2,A,T66,I12,A)') "FLOQUET MEMORY | Free memory per rank:", &
609 2 : (MemLikelyFree/INT(ranks_per_node, int_8))/1048576_int_8, " MB"
610 1 : WRITE (unit_nr, '(T2,A,T66,I12,A)') "FLOQUET MEMORY | Usable per rank (reserve applied):", &
611 2 : usable_per_rank/1048576_int_8, " MB"
612 1 : WRITE (unit_nr, '(T2,A,T66,I15)') "FLOQUET MEMORY | Floquet copies needed:", n_fm_work_copies
613 1 : WRITE (unit_nr, '(T2,A,T66,I12,A)') "FLOQUET MEMORY | Total Floquet working set:", &
614 2 : needed_bytes/1048576_int_8, " MB"
615 1 : WRITE (unit_nr, '(T2,A,T66,I12,A)') "FLOQUET MEMORY | Floquet working set per rank:", &
616 2 : (needed_bytes/INT(MAX(group_size, 1), int_8))/1048576_int_8, " MB"
617 1 : WRITE (unit_nr, '(T2,A,T66,I15)') "FLOQUET MEMORY | Memory floor (ranks/subgroup):", g_mem
618 1 : WRITE (unit_nr, '(T2,A,T66,I15)') "FLOQUET MEMORY | K-point floor (ranks/subgroup):", g_load
619 1 : WRITE (unit_nr, '(T2,A,T66,I15)') "FLOQUET MEMORY | Chosen MPI ranks per subgroup:", group_size
620 : END IF
621 :
622 2 : CALL timestop(handle)
623 2 : END FUNCTION floquet_determine_subgroup_size
624 :
625 : ! **************************************************************************************************
626 : !> \brief Split the para_env into subgroups so that the Floquet Hamiltonian is distributed across
627 : !> the ranks of each subgroup. Also creates the subgroup BLACS context.
628 : !> \param qs_env ...
629 : !> \param bs_env ...
630 : !> \param para_env_sub the created subgroup parallel environment
631 : !> \param blacs_env_sub the created subgroup BLACS context
632 : !> \param group_distribution subgroup index of every global rank, 0:num_pe-1
633 : !> \param ngroups the number of subgroups created
634 : ! **************************************************************************************************
635 2 : SUBROUTINE make_floquet_subgroups(qs_env, bs_env, para_env_sub, blacs_env_sub, &
636 : group_distribution, ngroups)
637 : TYPE(qs_environment_type), POINTER :: qs_env
638 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
639 : TYPE(mp_para_env_type), POINTER :: para_env_sub
640 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
641 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: group_distribution
642 : INTEGER, INTENT(OUT) :: ngroups
643 :
644 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_floquet_subgroups'
645 :
646 : INTEGER :: group_size, handle, n_spin, nkp_only_bs, &
647 : stride_kp, unit_nr
648 : TYPE(mp_para_env_type), POINTER :: para_env
649 :
650 2 : CALL timeset(routineN, handle)
651 :
652 2 : CALL get_qs_env(qs_env, para_env=para_env)
653 2 : unit_nr = bs_env%unit_nr
654 2 : nkp_only_bs = bs_env%nkp_only_bs
655 2 : n_spin = bs_env%n_spin
656 :
657 : ! We split the global communicator into subgroups with the Floquet Hamiltonian distributed
658 : ! across the processes in each subgroup and the k-points distributed across subgroups.
659 2 : group_size = floquet_determine_subgroup_size(qs_env, bs_env)
660 :
661 2 : IF (nkp_only_bs < para_env%num_pe) THEN
662 2 : stride_kp = para_env%num_pe/group_size
663 : ELSE
664 0 : stride_kp = 1
665 : END IF
666 6 : ALLOCATE (group_distribution(0:para_env%num_pe - 1))
667 2 : ALLOCATE (para_env_sub)
668 : CALL para_env_sub%from_split(comm=para_env, ngroups=ngroups, &
669 : group_distribution=group_distribution, &
670 2 : subgroup_min_size=group_size, stride=stride_kp)
671 2 : IF (unit_nr > 0) THEN
672 1 : WRITE (unit_nr, '(T2,A,T66,I15)') "FLOQUET MEMORY | Subgroup rank stride:", stride_kp
673 1 : WRITE (unit_nr, '(T2,A,T66,I15)') "FLOQUET MEMORY | Number of subgroups:", ngroups
674 1 : WRITE (unit_nr, '(/,T2,A,I5,A,I1,A)') "FLOQUET CALCULATIONS PROGRESS OUT OF", nkp_only_bs, &
675 2 : " K-POINTS AND ", n_spin, " SPINS"
676 : END IF
677 :
678 2 : NULLIFY (blacs_env_sub)
679 2 : CALL cp_blacs_env_create(blacs_env_sub, para_env_sub)
680 :
681 2 : CALL timestop(handle)
682 :
683 2 : END SUBROUTINE make_floquet_subgroups
684 :
685 : ! **************************************************************************************************
686 : !> \brief Central and boundary-sector weights of every Floquet eigenvector, stored into
687 : !> floquet_env%w0 and floquet_env%wE.
688 : !> \param floquet_env ...
689 : !> \param cfm_eigenvectors Floquet eigenvectors
690 : ! **************************************************************************************************
691 2 : SUBROUTINE floquet_sector_weights(floquet_env, cfm_eigenvectors)
692 : TYPE(floquet_env_type), INTENT(INOUT) :: floquet_env
693 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_eigenvectors
694 :
695 : CHARACTER(LEN=*), PARAMETER :: routineN = 'floquet_sector_weights'
696 :
697 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: v_block
698 : INTEGER :: handle, max_f_index, n_f_size, nao
699 :
700 2 : CALL timeset(routineN, handle)
701 :
702 2 : nao = floquet_env%nao
703 2 : max_f_index = floquet_env%max_f_index
704 2 : n_f_size = floquet_env%n_f_size
705 :
706 8 : ALLOCATE (v_block(nao, n_f_size))
707 :
708 : ! Evaluate w0_α = Σ_{n=1..nao} |<n,m=0|α>|^2 (central sector)
709 2 : CALL cp_cfm_get_submatrix(cfm_eigenvectors, v_block, 1 + nao*max_f_index, 1)
710 14546 : floquet_env%w0(:) = SUM(ABS(v_block)**2, DIM=1)
711 :
712 : ! Evaluate wE_α = Σ_{n=1..nao} (|<n,m=-M|α>|^2 + |<n,m=+M|α>|^2) (outermost rungs)
713 1618 : floquet_env%wE(:) = 0.0_dp
714 2 : IF (max_f_index >= 1) THEN
715 2 : CALL cp_cfm_get_submatrix(cfm_eigenvectors, v_block, 1, 1)
716 14546 : floquet_env%wE(:) = SUM(ABS(v_block)**2, DIM=1)
717 :
718 2 : CALL cp_cfm_get_submatrix(cfm_eigenvectors, v_block, 1 + nao*2*max_f_index, 1)
719 14546 : floquet_env%wE(:) = floquet_env%wE(:) + SUM(ABS(v_block)**2, DIM=1)
720 : END IF
721 :
722 2 : DEALLOCATE (v_block)
723 :
724 : ! The orthonormal rows of a unitary matrix are such that Σ_α w0_α = nao.
725 1618 : CPASSERT(ABS(SUM(floquet_env%w0) - REAL(nao, dp)) < 1.0E-6_dp*REAL(nao, dp))
726 :
727 2 : CALL timestop(handle)
728 :
729 2 : END SUBROUTINE floquet_sector_weights
730 :
731 : ! **************************************************************************************************
732 : !> \brief Checks that MAX_FLOQUET_INDEX is large enough and the Floquet Hamiltonian was
733 : !> truncated far enough away from the central sector to leave it unperturbed.
734 : !> \param bs_env ...
735 : !> \param floquet_env holds the central- (w0) and outermost-sector (wE) weights
736 : ! **************************************************************************************************
737 2 : SUBROUTINE check_floquet_convergence(bs_env, floquet_env)
738 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
739 : TYPE(floquet_env_type), INTENT(IN) :: floquet_env
740 :
741 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_floquet_convergence'
742 :
743 : CHARACTER(LEN=default_string_length) :: msg
744 : INTEGER :: handle
745 : REAL(KIND=dp) :: leak
746 :
747 : ! The truncation is converged if m=0 states have small weight on outermost rungs m = ±M
748 2 : CALL timeset(routineN, handle)
749 :
750 : ! MAX_FLOQUET_INDEX = 0 means no sidebands at all
751 2 : IF (bs_env%max_floquet_index < 1) THEN
752 0 : CALL timestop(handle)
753 0 : RETURN
754 : END IF
755 :
756 : ! No cleanly m=0-dominated state exists: the drive has hybridised every band with
757 : ! its sidebands, strong field but not a truncation failure.
758 630 : IF (.NOT. ANY(floquet_env%w0 > 0.5_dp)) THEN
759 0 : CPWARN("Floquet: no m=0-dominated state; cannot assess MAX_FLOQUET_INDEX convergence.")
760 0 : CALL timestop(handle)
761 0 : RETURN
762 : END IF
763 :
764 1618 : leak = MAXVAL(floquet_env%wE, MASK=(floquet_env%w0 > 0.5_dp))
765 :
766 2 : IF (bs_env%eps_floquet > 0.0_dp .AND. leak > bs_env%eps_floquet) THEN
767 : WRITE (msg, '(A,ES10.2E2,A,ES10.2E2)') &
768 0 : "MAX_FLOQUET_INDEX is too small. Leak: ", leak, "exceeds EPS_FLOQUET: ", bs_env%eps_floquet
769 0 : CPABORT(TRIM(msg))
770 : END IF
771 :
772 2 : CALL timestop(handle)
773 :
774 : END SUBROUTINE check_floquet_convergence
775 :
776 : ! **************************************************************************************************
777 : !> \brief Energy origin for the Floquet calculations, chosen to be the valence band maximum.
778 : !> \param bs_env ...
779 : !> \return the valence band maximum
780 : ! **************************************************************************************************
781 3 : FUNCTION floquet_reference_energy(bs_env) RESULT(mu)
782 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
783 : REAL(KIND=dp) :: mu
784 :
785 3 : IF (bs_env%do_gw .OR. &
786 : bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
787 3 : mu = bs_env%band_edges_scf%VBM
788 0 : ELSE IF (bs_env%n_spin == 1) THEN
789 0 : mu = bs_env%eigenval_scf_Gamma(bs_env%n_occ(1), 1)
790 : ELSE
791 : mu = MAX(bs_env%eigenval_scf_Gamma(bs_env%n_occ(1), 1), &
792 0 : bs_env%eigenval_scf_Gamma(bs_env%n_occ(2), 2))
793 : END IF
794 :
795 3 : END FUNCTION floquet_reference_energy
796 :
797 : ! **************************************************************************************************
798 : !> \brief Calculate all Floquet observables for one k-point and spin: the DOS (a_k), the m=0 bands,
799 : !> the quasi-energies, and store in floquet_env.
800 : !> \param bs_env ...
801 : !> \param para_env_sub ...
802 : !> \param ispin ...
803 : !> \param ikp ...
804 : !> \param floquet_env holds eigenvalues, w0 (input) and a_k, m0_energies, m0_weights,
805 : !> quasi_energies (output)
806 : ! **************************************************************************************************
807 2 : SUBROUTINE calculate_floquet_observables(bs_env, para_env_sub, ispin, ikp, floquet_env)
808 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
809 : TYPE(mp_para_env_type), POINTER :: para_env_sub
810 : INTEGER, INTENT(IN) :: ispin, ikp
811 : TYPE(floquet_env_type), INTENT(INOUT) :: floquet_env
812 :
813 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_floquet_observables'
814 :
815 : INTEGER :: handle, i, i_E, j, n_E, n_f_size, nao
816 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: work
817 : REAL(KIND=dp) :: broad, cum, E_min, energy, energy_step, &
818 : mu, omega, target_level
819 :
820 2 : CALL timeset(routineN, handle)
821 :
822 2 : n_f_size = floquet_env%n_f_size
823 2 : nao = floquet_env%nao
824 2 : n_E = floquet_env%n_E
825 :
826 2 : mu = floquet_reference_energy(bs_env)
827 2 : omega = bs_env%floquet_omega
828 2 : broad = bs_env%broadening_floquet
829 2 : energy_step = bs_env%energy_step_floquet
830 2 : E_min = mu - bs_env%energy_window_floquet
831 :
832 : ! (1) Spectral function / DOS, A(k,E) = -(1/π) Im Tr_KS G^R_00(k,E),
833 : ! where G^R_00(k,E) = (E + iη - H_F)^(-1) is the retarded Green's function
834 : ! With H_F = Σ_α λ_α |α><α|, we simplify
835 : ! A(k,E) = -(1/π) Σ_α w0_α Im[1/(E + iη - λ_α)], η = broadening/2.
836 : ! This avoids the per-energy dense inversion of the full H_F.
837 :
838 6 : DO i_E = 1, n_E
839 4 : energy = E_min + i_E*energy_step
840 : floquet_env%a_k(i_E) = -SUM(floquet_env%w0(:)* &
841 : AIMAG(z_one/(energy + gaussi*broad/2.0_dp - &
842 3238 : floquet_env%eigenvalues(:))))/pi
843 : END DO
844 :
845 : ! (2) m=0 band structure by weighted count
846 : cum = 0.0_dp
847 : j = 0
848 18 : DO i = 1, nao
849 16 : target_level = REAL(i, dp) - 0.5_dp
850 1052 : DO WHILE (cum < target_level .AND. j < n_f_size)
851 1036 : j = j + 1
852 1036 : cum = cum + floquet_env%w0(j)
853 : END DO
854 16 : floquet_env%m0_energies(i) = floquet_env%eigenvalues(j) - mu
855 18 : floquet_env%m0_weights(i) = floquet_env%w0(j)
856 : END DO
857 :
858 : ! (3) Fold the m=0 bands (already relative to the VBM) into the first Floquet Brillouin zone.
859 : floquet_env%quasi_energies(:) = floquet_env%m0_energies &
860 18 : - omega*REAL(CEILING(floquet_env%m0_energies/omega - 0.5_dp), dp)
861 :
862 : ! Sort the quasi-energies
863 6 : ALLOCATE (work(nao))
864 2 : CALL sort(floquet_env%quasi_energies, nao, work)
865 2 : DEALLOCATE (work)
866 :
867 : ! Store this result on the subgroup source only, so the global sum below picks up
868 : ! each work item only once
869 2 : IF (para_env_sub%is_source()) THEN
870 9 : floquet_env%all_quasi_energies(:, ispin, ikp) = floquet_env%quasi_energies(:)
871 3 : floquet_env%all_a_k(:, ispin, ikp) = floquet_env%a_k(:)
872 9 : floquet_env%all_m0_energies(:, ispin, ikp) = floquet_env%m0_energies(:)
873 9 : floquet_env%all_m0_weights(:, ispin, ikp) = floquet_env%m0_weights(:)
874 : END IF
875 :
876 2 : CALL timestop(handle)
877 :
878 2 : END SUBROUTINE calculate_floquet_observables
879 :
880 : ! **************************************************************************************************
881 : !> \brief Print the Floquet header with the input parameters and a short description of the output.
882 : !> \param bs_env ...
883 : ! **************************************************************************************************
884 2 : SUBROUTINE write_floquet_header(bs_env)
885 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
886 :
887 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_floquet_header'
888 :
889 : INTEGER :: handle, n_f_size, unit_nr
890 :
891 2 : CALL timeset(routineN, handle)
892 :
893 2 : unit_nr = bs_env%unit_nr
894 2 : n_f_size = bs_env%n_ao*(1 + 2*bs_env%max_floquet_index)
895 :
896 2 : IF (unit_nr > 0) THEN
897 :
898 1 : WRITE (unit_nr, '(T2,A)') ' '
899 1 : WRITE (unit_nr, '(T2,A)') REPEAT('-', 79)
900 1 : WRITE (unit_nr, '(T2,A,A78)') '-', '-'
901 1 : WRITE (unit_nr, '(T2,A,A51,A27)') '-', 'FLOQUET BANDSTRUCTURE CALCULATION', '-'
902 1 : WRITE (unit_nr, '(T2,A,A78)') '-', '-'
903 1 : WRITE (unit_nr, '(T2,A)') REPEAT('-', 79)
904 1 : WRITE (unit_nr, '(T2,A)') ' '
905 :
906 1 : WRITE (unit_nr, '(T2,A,T37,A,T67,ES12.2E2)') "FLOQUET PARAMETERS", "Amplitude [V/m]:", &
907 2 : bs_env%floquet_amplitude*evolt/a_bohr
908 1 : WRITE (unit_nr, '(T37,A,T67,F12.4)') "Frequency [eV]:", bs_env%floquet_omega*evolt
909 1 : WRITE (unit_nr, '(T37,A)') REPEAT("-", 44)
910 4 : WRITE (unit_nr, '(T37,A,T55,3F8.4)') "Polarisation:", bs_env%floquet_polarisation(1:3)
911 4 : WRITE (unit_nr, '(T37,A,T55,3F8.4)') "Phase offsets:", pi*bs_env%floquet_phi(1:3)
912 1 : WRITE (unit_nr, '(T37,A)') REPEAT("-", 44)
913 1 : WRITE (unit_nr, '(T37,A,T67,I12)') "Max Floquet index:", bs_env%max_floquet_index
914 1 : WRITE (unit_nr, '(T37,A,T67,I12)') "Floquet Hamiltonian Size:", n_f_size
915 1 : WRITE (unit_nr, '(T37,A)') REPEAT("-", 44)
916 : WRITE (unit_nr, '(T37,A,T67,F12.4)') &
917 1 : "Energy window [eV]:", bs_env%energy_window_floquet*evolt
918 1 : WRITE (unit_nr, '(T37,A,T67,F12.4)') "Energy step [eV]:", bs_env%energy_step_floquet*evolt
919 1 : WRITE (unit_nr, '(T37,A,T67,F12.4)') "Broadening [eV]:", bs_env%broadening_floquet*evolt
920 1 : WRITE (unit_nr, '(T37,A)') REPEAT("-", 44)
921 1 : WRITE (unit_nr, '(A)') ""
922 :
923 : WRITE (unit_nr, '(T2,A)') &
924 1 : "We construct the Floquet-Bloch Hamiltonian and diagonalise it. Projecting the"
925 : WRITE (unit_nr, '(T2,A)') &
926 1 : "eigenvectors onto the Floquet sectors gives the weights w = Σ_n |<n,m|α>|²,"
927 : WRITE (unit_nr, '(T2,A)') &
928 1 : "from which all of the following are obtained."
929 1 : WRITE (unit_nr, '(A)') ""
930 : WRITE (unit_nr, '(T2,A)') &
931 1 : "The m=0 bands are the eigenvectors with the largest central-sector weight; they"
932 : WRITE (unit_nr, '(T2,A)') &
933 1 : "reduce to the equilibrium bands at zero field and are stored, with their"
934 : WRITE (unit_nr, '(T2,A)') &
935 1 : " weights, in FLOQUET_BANDSTRUCTURE.bs"
936 1 : WRITE (unit_nr, '(A)') ""
937 : WRITE (unit_nr, '(T2,A)') &
938 1 : "Folding those bands into the first Floquet Brillouin zone, relative to the VBM,"
939 : WRITE (unit_nr, '(T2,A)') &
940 1 : "gives the quasi-energies stored in QUASI_ENERGIES.bs"
941 1 : WRITE (unit_nr, '(A)') ""
942 : WRITE (unit_nr, '(T2,A)') &
943 1 : "The k-resolved density of states is obtained by computing the trace of"
944 : WRITE (unit_nr, '(T2,A)') &
945 1 : "the retarded Green's function and stored in FLOQUET_DOS.out"
946 : WRITE (unit_nr, '(T2,A)') &
947 1 : "DOS(ω,k) = -1/π*Im[Tr_KS(G^R(ω,k))]"
948 1 : WRITE (unit_nr, '(A)') ""
949 : END IF
950 :
951 2 : CALL timestop(handle)
952 :
953 2 : END SUBROUTINE write_floquet_header
954 :
955 : ! **************************************************************************************************
956 : !> \brief Sum the accumulated results across MPI ranks and write the m=0 band structure,
957 : !> the quasi-energies and the Floquet DOS to their files.
958 : !> \param bs_env ...
959 : !> \param floquet_env ...
960 : ! **************************************************************************************************
961 2 : SUBROUTINE write_floquet_results(bs_env, floquet_env)
962 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
963 : TYPE(floquet_env_type), INTENT(INOUT) :: floquet_env
964 :
965 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_floquet_results'
966 :
967 : CHARACTER(LEN=default_string_length) :: fname
968 : INTEGER :: bunit, handle, i, i_E, ikp_for_file, &
969 : ispin, n_E, n_spin, nao, nkp_only_bs, &
970 : nkp_start, qunit, wunit
971 : REAL(KIND=dp) :: E_min, energy, energy_step, f_occ, kT, &
972 : mu, x
973 : REAL(KIND=dp), DIMENSION(3) :: xkp
974 :
975 2 : CALL timeset(routineN, handle)
976 :
977 : ! Collect all results on the global communicator and write them
978 2 : CALL bs_env%para_env%sum(floquet_env%all_quasi_energies)
979 2 : CALL bs_env%para_env%sum(floquet_env%all_a_k)
980 2 : CALL bs_env%para_env%sum(floquet_env%all_m0_energies)
981 2 : CALL bs_env%para_env%sum(floquet_env%all_m0_weights)
982 :
983 2 : IF (bs_env%para_env%is_source()) THEN
984 :
985 1 : nao = floquet_env%nao
986 1 : n_spin = floquet_env%n_spin
987 1 : nkp_only_bs = floquet_env%nkp_only_bs
988 1 : n_E = floquet_env%n_E
989 1 : nkp_start = bs_env%nkp_only_DOS
990 :
991 : ! All three share one energy origin, the VBM.
992 1 : mu = floquet_reference_energy(bs_env)
993 1 : energy_step = bs_env%energy_step_floquet
994 1 : E_min = mu - bs_env%energy_window_floquet
995 :
996 : ! k_B T in Hartree (only used when TEMPERATURE > 0); E[K] = E[Hartree]*kelvin.
997 1 : kT = bs_env%floquet_temperature/kelvin
998 :
999 : ! Each file is opened once (REPLACE) and written for every band-structure k-point and spin.
1000 :
1001 : ! The m=0 band structure
1002 1 : WRITE (fname, "(2A)") TRIM(bs_env%floquet_bs_file), ".bs"
1003 1 : CALL open_file(TRIM(fname), unit_number=bunit, file_status="REPLACE", file_action="WRITE")
1004 1 : WRITE (bunit, "(A)") "# Floquet m=0 (central-sector) band structure"
1005 1 : WRITE (bunit, "(A)") "# (in units of eV, relative to the VBM, not folded)"
1006 1 : WRITE (bunit, "(A)") "# w = sum_n |<n,m=0|alpha>|^2 in [0,1]: w~1 clean m=0 replica,"
1007 1 : WRITE (bunit, "(A)") "# w~0.5 hybridised with a sideband (drive near resonance)"
1008 :
1009 : ! Quasi-energies
1010 1 : WRITE (fname, "(2A)") TRIM(bs_env%floquet_qe_file), ".bs"
1011 1 : CALL open_file(TRIM(fname), unit_number=qunit, file_status="REPLACE", file_action="WRITE")
1012 1 : WRITE (qunit, "(A)") "# Quasi-energies obtained by diagonalising the Floquet Hamiltonian"
1013 1 : WRITE (qunit, "(A)") "# (in units of eV, the m=0 bands relative to the VBM, folded to"
1014 1 : WRITE (qunit, "(A)") "# the first Floquet Brillouin zone -hbar*Omega/2 < e <= hbar*Omega/2)"
1015 :
1016 : ! DOS
1017 1 : WRITE (fname, "(2A)") TRIM(bs_env%floquet_dos_file), ".out"
1018 1 : CALL open_file(TRIM(fname), unit_number=wunit, file_status="REPLACE", file_action="WRITE")
1019 1 : WRITE (wunit, "(A)") "# Floquet Density of States: D(ω,k) = -1/π*Im[Tr_KS(G^R(ω,k))]"
1020 :
1021 2 : DO ikp_for_file = 1, nkp_only_bs
1022 4 : xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, nkp_start + ikp_for_file)
1023 3 : DO ispin = 1, n_spin
1024 :
1025 : ! <floquet_bs_file>.bs the m=0 band structure (absolute energies) and its weights
1026 : WRITE (bunit, "(A,I0,T10,A,I0,A,T24,3(1X,F14.8))") &
1027 1 : "# Spin ", ispin, " Point ", ikp_for_file, ": ", xkp(1:3)
1028 1 : WRITE (bunit, "(A)") "# Floquet band Energy [eV] m=0 weight"
1029 9 : DO i = 1, nao
1030 8 : WRITE (bunit, "(I8,F21.8,F17.5)") i, &
1031 8 : floquet_env%all_m0_energies(i, ispin, ikp_for_file)*evolt, &
1032 17 : floquet_env%all_m0_weights(i, ispin, ikp_for_file)
1033 : END DO
1034 :
1035 : ! <floquet_qe_file>.bs Quasi-energies (folded to -ħΩ/2 < ε ≤ ħΩ/2)
1036 : WRITE (qunit, "(A,I0,T10,A,I0,A,T24,3(1X,F14.8))") &
1037 1 : "# Spin ", ispin, " Point ", ikp_for_file, ": ", xkp(1:3)
1038 1 : WRITE (qunit, "(A)") "# Floquet band Quasi-energy [eV]"
1039 9 : DO i = 1, nao
1040 8 : WRITE (qunit, "(I8,F21.8)") i, &
1041 17 : floquet_env%all_quasi_energies(i, ispin, ikp_for_file)*evolt
1042 : END DO
1043 :
1044 : ! <floquet_dos_file>.out the k-resolved DOS A(k,ω), and (if TEMPERATURE > 0) the
1045 : ! occupied spectral weight f(E)*A(k,ω), f the Fermi-Dirac occupation at the VBM.
1046 : WRITE (wunit, "(A,I0,T10,A,I0,A,T24,3(1X,F14.8))") &
1047 1 : "# Spin ", ispin, " Point ", ikp_for_file, ": ", xkp(1:3)
1048 2 : IF (bs_env%floquet_temperature > 0.0_dp) THEN
1049 1 : WRITE (wunit, "(A)") "#Energy-VBM (eV) A(ω,k) = DOS (1/eV) f*A = occupied DOS (1/eV)"
1050 3 : DO i_E = 1, n_E
1051 2 : energy = E_min + i_E*energy_step
1052 2 : x = (energy - mu)/kT ! (E - VBM)/k_B T, dimensionless
1053 2 : IF (x > 40.0_dp) THEN ! guard EXP overflow at low T
1054 : f_occ = 0.0_dp
1055 2 : ELSE IF (x < -40.0_dp) THEN
1056 : f_occ = 1.0_dp
1057 : ELSE
1058 2 : f_occ = 1.0_dp/(EXP(x) + 1.0_dp)
1059 : END IF
1060 2 : WRITE (wunit, "(2X,3G13.4)") (energy - mu)*evolt, &
1061 2 : floquet_env%all_a_k(i_E, ispin, ikp_for_file)/evolt, &
1062 5 : f_occ*floquet_env%all_a_k(i_E, ispin, ikp_for_file)/evolt
1063 : END DO
1064 : ELSE
1065 0 : WRITE (wunit, "(A)") "#Energy-VBM (eV) A(ω,k) = DOS (1/eV)"
1066 0 : DO i_E = 1, n_E
1067 0 : energy = E_min + i_E*energy_step
1068 0 : WRITE (wunit, "(2X,2G13.4)") (energy - mu)*evolt, &
1069 0 : floquet_env%all_a_k(i_E, ispin, ikp_for_file)/evolt
1070 : END DO
1071 : END IF
1072 :
1073 : END DO
1074 : END DO
1075 :
1076 1 : CALL close_file(bunit)
1077 1 : CALL close_file(qunit)
1078 1 : CALL close_file(wunit)
1079 : END IF
1080 :
1081 2 : CALL timestop(handle)
1082 :
1083 2 : END SUBROUTINE write_floquet_results
1084 :
1085 : END MODULE floquet_utils
|
