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_de_dk_dipole, 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(:, :, :, :, :) :: dipole
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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: de_dk
70 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
71 : TYPE(cp_cfm_type) :: cfm_eigenvectors, floquet_matrix
72 : TYPE(cp_fm_struct_type), POINTER :: floquet_struct
73 2 : TYPE(floquet_env_type) :: floquet_env
74 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
75 :
76 2 : CALL timeset(routineN, handle)
77 :
78 2 : CALL get_qs_env(qs_env, para_env=para_env)
79 :
80 : ! Floquet-Bloch theory: under a monochromatic electric field drive the Hamiltonian is periodic,
81 : ! H(k,t) = H(k,t+T) with T = 2π/Ω, where Ω is the angular frequency of the drive.
82 : ! By Floquet's theorem the solutions are:
83 : ! ψ_αk(t) = e^(-i λ_αk t) Σ_m e^(-i m Ω t) |F^m_αk>, with quasi-energies λ_αk and
84 : ! Fourier index m.
85 : ! In order to obtain the solutions λ_αk, we start with the Floquet Hamiltonian defined as,
86 : ! H_F(k,t) = H(k,t) - iħ d/dt.
87 : ! Taking the Fourier transform, the solutions are obtained by diagonalizing the matrix H_F(k)
88 : ! H_F(k) F_αk = λ_αk F_αk,
89 : ! where H_F(k) is block-tridiagonal in the sector index m = -M..M (M = MAX_FLOQUET_INDEX)
90 : ! and each block is nao×nao. We build and diagonalise H_F(k) per DOS k-point k and spin.
91 : ! The matrix H_F(k) has size (nao*(2*MAX_FLOQUET_INDEX+1) x nao*(2*MAX_FLOQUET_INDEX+1)).
92 :
93 2 : CALL write_floquet_header(bs_env)
94 :
95 2 : CALL floquet_env_create(floquet_env, bs_env)
96 :
97 : ! Equilibrium band data at every band-structure k-point k from H(k)C(k) = S(k)C(k)ε(k):
98 : ! band energies ε_nk
99 : ! k-derivatives ∇_k ε_nk = C_n(k)^† ∇_k H(k) C_n(k) - ε_nk C_n(k)^† ∇_k S(k) C_n(k)
100 : ! dipoles d_nm(k) = <ψ_nk| r |ψ_mk>
101 : ! Results are distributed across ranks (k-point ikp stored on rank MOD(ikp-1,num_pe)).
102 2 : CALL compute_e_k_de_dk_dipole(qs_env, bs_env, e_k, de_dk, dipole)
103 :
104 : ! Split the global communicator into Floquet subgroups (Hamiltonian distributed within each
105 : ! subgroup, k-points distributed across subgroups) and create the subgroup BLACS context.
106 : CALL make_floquet_subgroups(qs_env, bs_env, para_env_sub, blacs_env_sub, &
107 2 : group_distribution, ngroups)
108 :
109 : !get inside make_floquet_subroutine
110 2 : my_group = group_distribution(para_env%mepos)
111 :
112 : CALL cp_fm_struct_create(floquet_struct, context=blacs_env_sub, para_env=para_env_sub, &
113 2 : nrow_global=floquet_env%n_f_size, ncol_global=floquet_env%n_f_size)
114 2 : CALL cp_cfm_create(floquet_matrix, floquet_struct, set_zero=.TRUE.)
115 2 : CALL cp_cfm_create(cfm_eigenvectors, floquet_struct, set_zero=.TRUE.)
116 :
117 2 : n_source_done = 0
118 :
119 4 : DO ikp = 1, floquet_env%nkp_only_bs
120 : ! Only the subgroup that owns this k-point's data works on it
121 2 : IF (group_distribution(MOD(ikp - 1, para_env%num_pe)) /= my_group) CYCLE
122 :
123 4 : DO ispin = 1, floquet_env%n_spin
124 : ! Broadcast this k-point/spin's ε_nk, ∇_k ε_nk and d_nm(k) to every rank of the
125 : ! owning subgroup (its owner rank holds them from the distributed precompute above).
126 : CALL distribute_floquet_kp_data(para_env, para_env_sub, ispin, ikp, e_k, de_dk, &
127 2 : dipole, floquet_env)
128 :
129 : ! Assemble H_F(k) in Sambe space (blocks m,m' = -M..M, each nao×nao)
130 : ! diagonal (H_F)^{m,m}_{nn'} = (ε_nk + m ħΩ) δ_{nn'}
131 : ! off diagonal (H_F)^{m,m±1}_{nn'} = Σ_a p^a_{nn'}(k) (i E_a e^{±iφ_a})/(2Ω), a=x,y,z
132 : ! with momentum p^a_{nn'}(k) = i d^a_{nn'}(k)(ε_nk - ε_n'k) + ∇_{k_a} ε_nk δ_{nn'}
133 : ! and field E_a = |E| pol_a.
134 2 : CALL build_floquet_matrix(bs_env, floquet_env, floquet_matrix)
135 :
136 : ! Diagonalise H_F(k) F_αk = λ_αk F_αk -> quasi-energies λ_αk and eigenvector components
137 : ! F^{n,m}_αk = <n,m|αk> (band index n, sector m).
138 2 : CALL cp_cfm_set_all(cfm_eigenvectors, z_zero)
139 2 : CALL cp_cfm_heevd(floquet_matrix, cfm_eigenvectors, floquet_env%eigenvalues)
140 :
141 : ! Calculate sector weights,
142 : ! w_αk^(m) = Σ_{n=1}^{nao} |<n,m|αk>|^2 (with Σ_m w_αk^(m) = 1).
143 : ! Keep only the central weight w0_αk = w_αk^(0) and the outermost-rung weight
144 : ! wE_αk = w_αk^(-M) + w_αk^(+M)
145 2 : CALL floquet_sector_weights(floquet_env, cfm_eigenvectors)
146 :
147 : ! Truncation check: the m=0-dominated states (w0_αk > 1/2) must leak negligibly onto the
148 : ! outermost rungs, max_{α: w0_αk > 1/2} wE_αk ≤ EPS_FLOQUET, else M is too small.
149 2 : CALL check_floquet_convergence(bs_env, floquet_env)
150 :
151 : ! From λ_αk and w0_αk obtain, and store on the subgroup source for (ispin,ikp):
152 : ! DOS:
153 : ! A(k,E) = -(1/π) Im Tr_KS G^R_00(k,E) = -(1/π) Σ_α w0_αk Im[1/(E+iη-λ_αk)]
154 : ! with G^R_00 = [(E+iη) - H_F(k)]^{-1}_{m=0,m=0}, η = broadening/2.
155 : ! The central (m=0) block of the retarded Floquet Green's function gives the
156 : ! physical, time-averaged spectral function / density of states; see
157 : ! T. Oka and S. Kitamura, Annu. Rev. Condens. Matter Phys. 10, 387 (2019), and
158 : ! N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B 78, 235124 (2008).
159 : ! m=0 bands:
160 : ! ε^{0}_ik = λ_{α(i)}k, with band i placed at the eigenvalue where the cumulative
161 : ! m=0 weight Σ_{λ_αk<E} w0_αk first crosses i - 1/2.
162 : ! We know that: Σ_α w0_αk = nao. So we choose the matrix eigenvalue such that:
163 : ! N(ε) = Σ_{α : λ_αk ≤ ε} w0_αk first crosses i - 1/2 to find the m=0 band.
164 : ! Quasi-energy:
165 : ! ε̃_ik = ε^{0}_ik - mħΩ (choose m such that ε̃_ik is folded to (-ħΩ/2, ħΩ/2]).
166 4 : CALL calculate_floquet_observables(bs_env, para_env_sub, ispin, ikp, floquet_env)
167 : END DO
168 :
169 : ! Rough live K-Point progress
170 4 : IF (bs_env%unit_nr > 0) THEN
171 1 : n_source_done = n_source_done + 1
172 1 : n_est_done = MIN(n_source_done*ngroups, floquet_env%nkp_only_bs)
173 : WRITE (bs_env%unit_nr, '(T2,A,I6,A,F6.1,A)') &
174 1 : "FLOQUET PROGRESS | Estimated no. of k-points done: ", n_est_done, " (", &
175 2 : 100.0_dp*REAL(n_est_done, dp)/REAL(floquet_env%nkp_only_bs, dp), " %)"
176 : END IF
177 : END DO
178 :
179 : ! Write the m=0 bands ε^{0}_ik, the quasi-energies ε̃_ik and the DOS A(k,E) in their
180 : ! respective files
181 2 : CALL write_floquet_results(bs_env, floquet_env)
182 :
183 2 : DEALLOCATE (e_k, de_dk, dipole)
184 2 : DEALLOCATE (group_distribution)
185 2 : CALL floquet_env_release(floquet_env)
186 2 : CALL cp_cfm_release(floquet_matrix)
187 2 : CALL cp_cfm_release(cfm_eigenvectors)
188 2 : CALL cp_fm_struct_release(floquet_struct)
189 2 : CALL cp_blacs_env_release(blacs_env_sub)
190 2 : CALL mp_para_env_release(para_env_sub)
191 :
192 2 : CALL timestop(handle)
193 :
194 4 : END SUBROUTINE floquet
195 :
196 : END MODULE floquet_main
|
