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-Bloch band-structure calculations: builds and diagonalises the
10 : !> truncated Floquet-Bloch Hamiltonian H_F(k) for each DOS k-point and writes the m=0
11 : !> band structure, the quasi-energies and the Floquet density of states.
12 : !> \par History
13 : !> \author Shridhar Shanbhag (27.01.2026)
14 : ! **************************************************************************************************
15 : MODULE floquet_main
16 : USE cp_blacs_env, ONLY: cp_blacs_env_release,&
17 : cp_blacs_env_type
18 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
19 : USE cp_cfm_types, ONLY: cp_cfm_create,&
20 : cp_cfm_release,&
21 : cp_cfm_set_all,&
22 : cp_cfm_type
23 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 : cp_fm_struct_release,&
25 : cp_fm_struct_type
26 : USE floquet_types, ONLY: floquet_env_create,&
27 : floquet_env_release,&
28 : floquet_env_type
29 : USE floquet_utils, ONLY: &
30 : build_floquet_matrix, calculate_floquet_observables, check_floquet_convergence, &
31 : compute_e_k_p_k, distribute_floquet_kp_data, floquet_sector_weights, &
32 : make_floquet_subgroups, write_floquet_header, write_floquet_results
33 : USE kinds, ONLY: dp
34 : USE mathconstants, ONLY: z_zero
35 : USE message_passing, ONLY: mp_para_env_release,&
36 : mp_para_env_type
37 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
38 : USE qs_environment_types, ONLY: get_qs_env,&
39 : qs_environment_type
40 : #include "./base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 :
44 : PRIVATE
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'floquet_main'
47 :
48 : PUBLIC :: floquet
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Computes the Floquet-Bloch band structure, density of states and quasi-energies.
54 : !> \param qs_env ...
55 : !> \param bs_env ...
56 : ! **************************************************************************************************
57 2 : SUBROUTINE floquet(qs_env, bs_env)
58 : TYPE(qs_environment_type), POINTER :: qs_env
59 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
60 :
61 : CHARACTER(LEN=*), PARAMETER :: routineN = 'floquet'
62 :
63 : COMPLEX(KIND=dp), ALLOCATABLE, &
64 2 : DIMENSION(:, :, :, :, :) :: p_k
65 : INTEGER :: handle, ikp, ispin, my_group, &
66 : n_est_done, n_source_done, ngroups
67 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: group_distribution
68 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: e_k
69 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
70 : TYPE(cp_cfm_type) :: floquet_eigenvectors, floquet_matrix
71 : TYPE(cp_fm_struct_type), POINTER :: floquet_struct
72 2 : TYPE(floquet_env_type) :: floquet_env
73 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
74 :
75 2 : CALL timeset(routineN, handle)
76 :
77 2 : CALL get_qs_env(qs_env, para_env=para_env)
78 :
79 : ! Floquet-Bloch theory: under a monochromatic electric field drive the Hamiltonian is periodic,
80 : ! H(k,t) = H(k,t+T) with T = 2π/Ω, where Ω is the angular frequency of the drive.
81 : ! By Floquet's theorem the solutions are:
82 : ! ψ_αk(t) = e^(-i λ_αk t) Σ_m e^(-i m Ω t) |F^m_αk>, with quasi-energies λ_αk and
83 : ! Fourier index m.
84 : ! In order to obtain the solutions λ_αk, we start with the Floquet Hamiltonian defined as,
85 : ! H_F(k,t) = H(k,t) - iħ d/dt.
86 : ! Taking the Fourier transform, the solutions are obtained by diagonalizing the matrix H_F(k)
87 : ! H_F(k) F_αk = λ_αk F_αk,
88 : ! where H_F(k) is block-tridiagonal in the sector index m = -M..M (M = MAX_FLOQUET_INDEX)
89 : ! and each block is n_mo×n_mo. We build and diagonalise H_F(k) per DOS k-point k and spin.
90 : ! The matrix H_F(k) has size (n_mo*(2*MAX_FLOQUET_INDEX+1) x n_mo*(2*MAX_FLOQUET_INDEX+1)).
91 :
92 2 : CALL write_floquet_header(bs_env)
93 :
94 2 : CALL floquet_env_create(floquet_env, bs_env)
95 :
96 : ! Equilibrium band data at every band-structure k-point k from H(k)C(k) = S(k)C(k)ε(k):
97 : ! band energies: ε_nk
98 : ! momentum matrix elements:
99 : ! p_nm = <ψ_n|-iħ∇_r|ψ_m> = <u_n|(ħk - iħ∇_r)|u_m>
100 : ! Results are distributed across ranks (k-point ikp stored on rank MOD(ikp-1,num_pe)).
101 2 : CALL compute_e_k_p_k(qs_env, bs_env, e_k, p_k)
102 :
103 : ! Split the global communicator into Floquet subgroups (Hamiltonian distributed within each
104 : ! subgroup, k-points distributed across subgroups) and create the subgroup BLACS context.
105 : CALL make_floquet_subgroups(qs_env, bs_env, para_env_sub, blacs_env_sub, &
106 2 : group_distribution, ngroups)
107 :
108 : !get inside make_floquet_subroutine
109 2 : my_group = group_distribution(para_env%mepos)
110 :
111 : CALL cp_fm_struct_create(floquet_struct, context=blacs_env_sub, para_env=para_env_sub, &
112 2 : nrow_global=floquet_env%n_f_size, ncol_global=floquet_env%n_f_size)
113 2 : CALL cp_cfm_create(floquet_matrix, floquet_struct, set_zero=.TRUE.)
114 2 : CALL cp_cfm_create(floquet_eigenvectors, floquet_struct, set_zero=.TRUE.)
115 :
116 2 : n_source_done = 0
117 :
118 4 : DO ikp = 1, floquet_env%nkp_only_bs
119 : ! Only the subgroup that owns this k-point's data works on it
120 2 : IF (group_distribution(MOD(ikp - 1, para_env%num_pe)) /= my_group) CYCLE
121 :
122 4 : DO ispin = 1, floquet_env%n_spin
123 : ! Broadcast this k-point/spin's ε_nk, ∇_k ε_nk and d_nm(k) to every rank of the
124 : ! owning subgroup (its owner rank holds them from the distributed precompute above).
125 : CALL distribute_floquet_kp_data(para_env, para_env_sub, ispin, ikp, e_k, &
126 2 : p_k, floquet_env)
127 :
128 : ! Assemble H_F(k) in Sambe space (blocks m,m' = -M..M, each n_mo×n_mo)
129 : ! diagonal (H_F)^{m,m}_{nn'} = (ε_nk + m ħΩ) δ_{nn'}
130 : ! off diagonal (H_F)^{m,m±1}_{nn'} = Σ_a p^a_{nn'}(k) (i E_a e^{±iφ_a})/(2Ω), a=x,y,z
131 : ! with momentum p^a_{nn'}(k) = i d^a_{nn'}(k)(ε_nk - ε_n'k) + ∇_{k_a} ε_nk δ_{nn'}
132 : ! and field E_a = |E| pol_a.
133 2 : CALL build_floquet_matrix(bs_env, floquet_env, floquet_matrix)
134 :
135 : ! Diagonalise H_F(k) F_αk = λ_αk F_αk -> quasi-energies λ_αk and eigenvector components
136 : ! F^{n,m}_αk = <n,m|αk> (band index n, sector m).
137 2 : CALL cp_cfm_set_all(floquet_eigenvectors, z_zero)
138 2 : CALL cp_cfm_heevd(floquet_matrix, floquet_eigenvectors, floquet_env%eigenvalues)
139 :
140 : ! Calculate sector weights,
141 : ! w_αk^(m) = Σ_{n=1}^{n_mo} |<n,m|αk>|^2 (with Σ_m w_αk^(m) = 1).
142 : ! Keep only the central weight w0_αk = w_αk^(0) and the outermost-rung weight
143 : ! wE_αk = w_αk^(-M) + w_αk^(+M)
144 2 : CALL floquet_sector_weights(floquet_env, floquet_eigenvectors)
145 :
146 : ! Truncation check: the m=0-dominated states (w0_αk > 1/2) must leak negligibly onto the
147 : ! outermost rungs, max_{α: w0_αk > 1/2} wE_αk ≤ EPS_FLOQUET, else M is too small.
148 2 : CALL check_floquet_convergence(bs_env, floquet_env)
149 :
150 : ! From λ_αk and w0_αk obtain, and store on the subgroup source for (ispin,ikp):
151 : ! DOS:
152 : ! A(k,E) = -(1/π) Im Tr_KS G^R_00(k,E) = -(1/π) Σ_α w0_αk Im[1/(E+iη-λ_αk)]
153 : ! with G^R_00 = [(E+iη) - H_F(k)]^{-1}_{m=0,m=0}, η = broadening/2.
154 : ! The central (m=0) block of the retarded Floquet Green's function gives the
155 : ! physical, time-averaged spectral function / density of states; see
156 : ! T. Oka and S. Kitamura, Annu. Rev. Condens. Matter Phys. 10, 387 (2019), and
157 : ! N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B 78, 235124 (2008).
158 : ! m=0 bands:
159 : ! ε^{0}_ik = λ_{α(i)}k, with band i placed at the eigenvalue where the cumulative
160 : ! m=0 weight Σ_{λ_αk<E} w0_αk first crosses i - 1/2.
161 : ! We know that: Σ_α w0_αk = n_mo. So we choose the matrix eigenvalue such that:
162 : ! N(ε) = Σ_{α : λ_αk ≤ ε} w0_αk first crosses i - 1/2 to find the m=0 band.
163 : ! Quasi-energy:
164 : ! ε̃_ik = ε^{0}_ik - mħΩ (choose m such that ε̃_ik is folded to (-ħΩ/2, ħΩ/2]).
165 4 : CALL calculate_floquet_observables(bs_env, para_env_sub, ispin, ikp, floquet_env)
166 : END DO
167 :
168 : ! Rough live K-Point progress
169 4 : IF (bs_env%unit_nr > 0) THEN
170 1 : n_source_done = n_source_done + 1
171 1 : n_est_done = MIN(n_source_done*ngroups, floquet_env%nkp_only_bs)
172 : WRITE (bs_env%unit_nr, '(T2,A,I6,A,F6.1,A)') &
173 1 : "FLOQUET PROGRESS | Estimated no. of k-points done: ", n_est_done, " (", &
174 2 : 100.0_dp*REAL(n_est_done, dp)/REAL(floquet_env%nkp_only_bs, dp), " %)"
175 : END IF
176 : END DO
177 :
178 : ! Write the m=0 bands ε^{0}_ik, the quasi-energies ε̃_ik and the DOS A(k,E) in their
179 : ! respective files
180 2 : CALL write_floquet_results(bs_env, floquet_env)
181 :
182 2 : DEALLOCATE (e_k, p_k)
183 2 : DEALLOCATE (group_distribution)
184 2 : CALL floquet_env_release(floquet_env)
185 2 : CALL cp_cfm_release(floquet_matrix)
186 2 : CALL cp_cfm_release(floquet_eigenvectors)
187 2 : CALL cp_fm_struct_release(floquet_struct)
188 2 : CALL cp_blacs_env_release(blacs_env_sub)
189 2 : CALL mp_para_env_release(para_env_sub)
190 :
191 2 : CALL timestop(handle)
192 :
193 4 : END SUBROUTINE floquet
194 :
195 : END MODULE floquet_main
|
