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