Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Floquet stuff
10 : !> \par History
11 : !> \author Shridhar Shanbhag (27.01.2026)
12 : ! **************************************************************************************************
13 : MODULE floquet_main
14 : USE cp_blacs_env, ONLY: cp_blacs_env_type
15 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
16 : USE cp_cfm_types, ONLY: cp_cfm_create,&
17 : cp_cfm_release,&
18 : cp_cfm_set_all,&
19 : cp_cfm_to_cfm,&
20 : cp_cfm_type
21 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22 : cp_fm_struct_release,&
23 : cp_fm_struct_type
24 : USE floquet_utils, ONLY: build_floquet_matrix,&
25 : calculate_floquet_spectral_function,&
26 : check_floquet_conversion,&
27 : write_floquet_spectral_function,&
28 : write_quasi_energy
29 : USE kinds, ONLY: dp
30 : USE mathconstants, ONLY: pi,&
31 : z_zero
32 : USE message_passing, ONLY: mp_para_env_type
33 : USE physcon, ONLY: a_bohr,&
34 : evolt
35 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type
38 : USE qs_mo_types, ONLY: get_mo_set,&
39 : mo_set_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 subroutines
49 : PUBLIC :: floquet
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief ...
55 : !> \param qs_env ...
56 : !> \param bs_env ...
57 : ! **************************************************************************************************
58 2 : SUBROUTINE floquet(qs_env, bs_env)
59 : TYPE(qs_environment_type), POINTER :: qs_env
60 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
61 :
62 : CHARACTER(LEN=*), PARAMETER :: routineN = 'floquet'
63 :
64 : INTEGER :: handle, ikp, ikp_for_file, ispin, &
65 : max_f_index, n_E, n_f_size, n_fbands, &
66 : n_spin, nao, nkp, nkp_start, unit_nr
67 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indices
68 : REAL(KIND=dp) :: energy_step, energy_window
69 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a_k, eigenvalues
70 : REAL(KIND=dp), DIMENSION(3) :: xkp
71 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
72 : TYPE(cp_cfm_type) :: cfm_eigenvectors, floquet_copy, &
73 : floquet_matrix
74 : TYPE(cp_fm_struct_type), POINTER :: floquet_struct
75 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
76 : TYPE(mp_para_env_type), POINTER :: para_env
77 :
78 2 : CALL timeset(routineN, handle)
79 :
80 : CALL get_qs_env(qs_env, &
81 : blacs_env=blacs_env, &
82 : para_env=para_env, &
83 2 : mos=mos)
84 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
85 :
86 2 : unit_nr = bs_env%unit_nr
87 2 : nkp = bs_env%nkp_bs_and_DOS
88 2 : nkp_start = bs_env%nkp_only_DOS
89 2 : n_spin = bs_env%n_spin
90 2 : max_f_index = bs_env%max_floquet_index
91 2 : n_fbands = 1 + 2*max_f_index
92 2 : n_f_size = nao*n_fbands
93 :
94 2 : energy_step = bs_env%energy_step_floquet
95 2 : energy_window = 2*bs_env%energy_window_floquet
96 :
97 2 : n_E = INT(energy_window/energy_step)
98 :
99 2 : IF (unit_nr > 0) THEN
100 :
101 1 : WRITE (unit_nr, '(T2,A)') ' '
102 1 : WRITE (unit_nr, '(T2,A)') REPEAT('-', 79)
103 1 : WRITE (unit_nr, '(T2,A,A78)') '-', '-'
104 1 : WRITE (unit_nr, '(T2,A,A51,A27)') '-', 'FLOQUET BANDSTRUCTURE CALCULATION', '-'
105 1 : WRITE (unit_nr, '(T2,A,A78)') '-', '-'
106 1 : WRITE (unit_nr, '(T2,A)') REPEAT('-', 79)
107 1 : WRITE (unit_nr, '(T2,A)') ' '
108 :
109 : WRITE (unit_nr, '(T2,A,T56,ES12.2E2)') &
110 1 : "FLOQUET PARAMETERS Amplitude [V/m]:", bs_env%floquet_amplitude*evolt/a_bohr
111 : WRITE (unit_nr, '(T2,A,T56,F12.4)') &
112 1 : " Frequency [eV]:", bs_env%floquet_omega*evolt
113 1 : WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
114 : WRITE (unit_nr, '(T2,A,T56,3F8.4)') &
115 4 : " Polarisation:", bs_env%floquet_polarisation(1:3)
116 : WRITE (unit_nr, '(T2,A,T56,3F8.4)') &
117 4 : " Phase offsets:", pi*bs_env%floquet_phi(1:3)
118 1 : WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
119 : WRITE (unit_nr, '(T2,A,T56,I12)') &
120 1 : " Max Floquet index:", bs_env%max_floquet_index
121 : WRITE (unit_nr, '(T2,A,T56,I12)') &
122 1 : " Floquet Hamiltonian Size:", n_f_size
123 1 : WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
124 : WRITE (unit_nr, '(T2,A,T56,F12.4)') &
125 1 : " Energy window [eV]:", bs_env%energy_window_floquet*evolt
126 : WRITE (unit_nr, '(T2,A,T56,F12.4)') &
127 1 : " Energy step [eV]:", bs_env%energy_step_floquet*evolt
128 : WRITE (unit_nr, '(T2,A,T56,F12.4)') &
129 1 : " Broadening [eV]:", bs_env%broadening_floquet*evolt
130 1 : WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
131 1 : WRITE (unit_nr, '(A)') ""
132 :
133 : WRITE (unit_nr, '(T2,A)') &
134 1 : "We construct the Floquet-Bloch Hamiltonian and diagonalise to obtain the"
135 : WRITE (unit_nr, '(T2,A)') &
136 1 : "eigenvalues which give us the quasi-energies stored in QUASI_ENERGIES.bs"
137 1 : WRITE (unit_nr, '(A)') ""
138 : WRITE (unit_nr, '(T2,A)') &
139 1 : "The k-resolved density of states is obtained by computing the trace of"
140 : WRITE (unit_nr, '(T2,A)') &
141 1 : "the retarded Green's function and stored in FLOQUET_DOS.out"
142 : WRITE (unit_nr, '(T2,A)') &
143 1 : "DOS(ω,k) = -1/π*Im[Tr_KS(G^R(ω,k))]"
144 1 : WRITE (unit_nr, '(A)') ""
145 1 : WRITE (unit_nr, '(T2,A)') "Calculations completed for:"
146 1 : WRITE (unit_nr, '(2X,A,T12,A,T22,A)') "KPoint", "Spin", "Coordinate"
147 : END IF
148 :
149 : CALL cp_fm_struct_create(floquet_struct, &
150 : nrow_global=n_f_size, &
151 : ncol_global=n_f_size, &
152 : context=blacs_env, &
153 2 : para_env=para_env)
154 :
155 2 : CALL cp_cfm_create(floquet_matrix, floquet_struct, set_zero=.TRUE.)
156 2 : CALL cp_cfm_create(floquet_copy, floquet_struct, set_zero=.TRUE.)
157 2 : CALL cp_cfm_create(cfm_eigenvectors, floquet_struct, set_zero=.TRUE.)
158 10 : ALLOCATE (eigenvalues(n_f_size), indices(nao))
159 6 : ALLOCATE (a_k(n_E))
160 4 : DO ikp = nkp_start + 1, nkp
161 :
162 8 : xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp)
163 2 : ikp_for_file = ikp - nkp_start
164 :
165 6 : DO ispin = 1, n_spin
166 : ! Build the Floquet-Bloch Hamiltonian matrix H_F(k), with the equilibrium
167 : ! band energies ε_{nk} on the diagonal (shifted by mħΩ per sector) and
168 : ! the light-matter coupling blocks A . p_uv on the off-diagonals
169 2 : CALL build_floquet_matrix(qs_env, bs_env, xkp, ispin, floquet_matrix)
170 :
171 : ! Use a copy of the Floquet matrix for diagonalization, since
172 : ! cp_cfm_heevd overwrites it (needed later for the Green's function)
173 2 : CALL cp_cfm_to_cfm(floquet_matrix, floquet_copy)
174 :
175 : ! Diagonalize the Floquet-Bloch matrix H_F(k) F_{α,k} = ε_α F_{α,k}
176 : ! to obtain the quasi-energies ε_α (eigenvalues) and the
177 : ! Floquet-Bloch coefficients F^{nu}_{α,k} (eigenvectors)
178 2 : CALL cp_cfm_set_all(cfm_eigenvectors, z_zero)
179 2 : CALL cp_cfm_heevd(floquet_copy, cfm_eigenvectors, eigenvalues)
180 :
181 : ! Run a conversion check to ensure that MAX_FLOQUET_INDEX,
182 : ! the size of the truncated matrix H_F(k), is sufficiently large.
183 2 : CALL check_floquet_conversion(qs_env, bs_env, cfm_eigenvectors)
184 :
185 : ! Write the quasi-energies ε_α folded into the first Floquet Brillouin
186 : ! zone -ħΩ/2 < ε_α ≤ ħΩ/2 for this k-point and spin channel
187 2 : CALL write_quasi_energy(qs_env, bs_env, ispin, ikp_for_file, xkp, eigenvalues)
188 :
189 : ! Compute the physical spectral function
190 : ! A_phys(k,ω) = -(1/π) Im[ Tr_KS G^R_{00}(k,ω) ]
191 : ! from the zero-zero Floquet sector of the retarded Green's function
192 2 : CALL calculate_floquet_spectral_function(qs_env, bs_env, floquet_matrix, a_k)
193 :
194 : ! Write A_phys(k,ω) (or f(ω,μ)·A_phys for the occupied spectral weight)
195 : ! to output for this k-point and spin channel
196 2 : CALL write_floquet_spectral_function(qs_env, bs_env, ispin, ikp_for_file, xkp, a_k)
197 :
198 : ! Write progress in the main output file
199 4 : IF (unit_nr > 0) WRITE (unit_nr, '(2X,I6,T12,I4,T18,3F12.6)') ikp_for_file, ispin, xkp(1:3)
200 :
201 : END DO
202 : END DO
203 :
204 2 : IF (unit_nr > 0) WRITE (unit_nr, '(A)') ""
205 2 : DEALLOCATE (eigenvalues)
206 2 : CALL cp_cfm_release(floquet_matrix)
207 2 : CALL cp_cfm_release(floquet_copy)
208 2 : CALL cp_cfm_release(cfm_eigenvectors)
209 2 : CALL cp_fm_struct_release(floquet_struct)
210 :
211 2 : CALL timestop(handle)
212 :
213 6 : END SUBROUTINE floquet
214 :
215 : END MODULE floquet_main
|
