comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 21:55:16 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Floquet physics (#5380)

Browse source
This commit is contained in:
Shridhar Shanbhag 2026-06-14 12:46:35 +02:00 • committed by GitHub
parent 74766aa487
commit f9604532c0
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
11 changed files with 1140 additions and 3 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/CMakeLists.txt
View file

@ -199,6 +199,8 @@ list(
fist_nonbond_env_types.F
fist_nonbond_force.F
fist_pol_scf.F
floquet_main.F
floquet_utils.F
force_env_methods.F
force_env_types.F
force_env_utils.F

215
src/floquet_main.F Normal file
View file

@ -0,0 +1,215 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Floquet stuff
!> \par History
!> \author Shridhar Shanbhag (27.01.2026)
! **************************************************************************************************
MODULE floquet_main
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_cfm_diag, ONLY: cp_cfm_heevd
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_release,&
cp_cfm_set_all,&
cp_cfm_to_cfm,&
cp_cfm_type
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE floquet_utils, ONLY: build_floquet_matrix,&
calculate_floquet_spectral_function,&
check_floquet_conversion,&
write_floquet_spectral_function,&
write_quasi_energy
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi,&
z_zero
USE message_passing, ONLY: mp_para_env_type
USE physcon, ONLY: a_bohr,&
evolt
USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'floquet_main'
! Public subroutines
PUBLIC :: floquet
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
SUBROUTINE floquet(qs_env, bs_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'floquet'
INTEGER :: handle, ikp, ikp_for_file, ispin, &
max_f_index, n_E, n_f_size, n_fbands, &
n_spin, nao, nkp, nkp_start, unit_nr
INTEGER, ALLOCATABLE, DIMENSION(:) :: indices
REAL(KIND=dp) :: energy_step, energy_window
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a_k, eigenvalues
REAL(KIND=dp), DIMENSION(3) :: xkp
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_cfm_type) :: cfm_eigenvectors, floquet_copy, &
floquet_matrix
TYPE(cp_fm_struct_type), POINTER :: floquet_struct
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, &
blacs_env=blacs_env, &
para_env=para_env, &
mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
unit_nr = bs_env%unit_nr
nkp = bs_env%nkp_bs_and_DOS
nkp_start = bs_env%nkp_only_DOS
n_spin = bs_env%n_spin
max_f_index = bs_env%max_floquet_index
n_fbands = 1 + 2*max_f_index
n_f_size = nao*n_fbands
energy_step = bs_env%energy_step_floquet
energy_window = 2*bs_env%energy_window_floquet
n_E = INT(energy_window/energy_step)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A)') ' '
WRITE (unit_nr, '(T2,A)') REPEAT('-', 79)
WRITE (unit_nr, '(T2,A,A78)') '-', '-'
WRITE (unit_nr, '(T2,A,A51,A27)') '-', 'FLOQUET BANDSTRUCTURE CALCULATION', '-'
WRITE (unit_nr, '(T2,A,A78)') '-', '-'
WRITE (unit_nr, '(T2,A)') REPEAT('-', 79)
WRITE (unit_nr, '(T2,A)') ' '
WRITE (unit_nr, '(T2,A,T56,ES12.2E2)') &
"FLOQUET PARAMETERS Amplitude [V/m]:", bs_env%floquet_amplitude*evolt/a_bohr
WRITE (unit_nr, '(T2,A,T56,F12.4)') &
" Frequency [eV]:", bs_env%floquet_omega*evolt
WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
WRITE (unit_nr, '(T2,A,T56,3F8.4)') &
" Polarisation:", bs_env%floquet_polarisation(1:3)
WRITE (unit_nr, '(T2,A,T56,3F8.4)') &
" Phase offsets:", pi*bs_env%floquet_phi(1:3)
WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
WRITE (unit_nr, '(T2,A,T56,I12)') &
" Max Floquet index:", bs_env%max_floquet_index
WRITE (unit_nr, '(T2,A,T56,I12)') &
" Floquet Hamiltonian Size:", n_f_size
WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
WRITE (unit_nr, '(T2,A,T56,F12.4)') &
" Energy window [eV]:", bs_env%energy_window_floquet*evolt
WRITE (unit_nr, '(T2,A,T56,F12.4)') &
" Energy step [eV]:", bs_env%energy_step_floquet*evolt
WRITE (unit_nr, '(T2,A,T56,F12.4)') &
" Broadening [eV]:", bs_env%broadening_floquet*evolt
WRITE (unit_nr, '(T2,A)') " "//REPEAT("-", 44)
WRITE (unit_nr, '(A)') ""
WRITE (unit_nr, '(T2,A)') &
"We construct the Floquet-Bloch Hamiltonian and diagonalise to obtain the"
WRITE (unit_nr, '(T2,A)') &
"eigenvalues which give us the quasi-energies stored in QUASI_ENERGIES.bs"
WRITE (unit_nr, '(A)') ""
WRITE (unit_nr, '(T2,A)') &
"The k-resolved density of states is obtained by computing the trace of"
WRITE (unit_nr, '(T2,A)') &
"the retarded Green's function and stored in FLOQUET_DOS.out"
WRITE (unit_nr, '(T2,A)') &
"DOS(ω,k) = -1/π*Im[Tr_KS(G^R(ω,k))]"
WRITE (unit_nr, '(A)') ""
WRITE (unit_nr, '(T2,A)') "Calculations completed for:"
WRITE (unit_nr, '(2X,A,T12,A,T22,A)') "KPoint", "Spin", "Coordinate"
END IF
CALL cp_fm_struct_create(floquet_struct, &
nrow_global=n_f_size, &
ncol_global=n_f_size, &
context=blacs_env, &
para_env=para_env)
CALL cp_cfm_create(floquet_matrix, floquet_struct, set_zero=.TRUE.)
CALL cp_cfm_create(floquet_copy, floquet_struct, set_zero=.TRUE.)
CALL cp_cfm_create(cfm_eigenvectors, floquet_struct, set_zero=.TRUE.)
ALLOCATE (eigenvalues(n_f_size), indices(nao))
ALLOCATE (a_k(n_E))
DO ikp = nkp_start + 1, nkp
xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp)
ikp_for_file = ikp - nkp_start
DO ispin = 1, n_spin
! Build the Floquet-Bloch Hamiltonian matrix H_F(k), with the equilibrium
! band energies ε_{nk} on the diagonal (shifted by mħΩ per sector) and
! the light-matter coupling blocks A . p_uv on the off-diagonals
CALL build_floquet_matrix(qs_env, bs_env, xkp, ispin, floquet_matrix)
! Use a copy of the Floquet matrix for diagonalization, since
! cp_cfm_heevd overwrites it (needed later for the Green's function)
CALL cp_cfm_to_cfm(floquet_matrix, floquet_copy)
! Diagonalize the Floquet-Bloch matrix H_F(k) F_{α,k} = ε_α F_{α,k}
! to obtain the quasi-energies ε_α (eigenvalues) and the
! Floquet-Bloch coefficients F^{nu}_{α,k} (eigenvectors)
CALL cp_cfm_set_all(cfm_eigenvectors, z_zero)
CALL cp_cfm_heevd(floquet_copy, cfm_eigenvectors, eigenvalues)
! Run a conversion check to ensure that MAX_FLOQUET_INDEX,
! the size of the truncated matrix H_F(k), is sufficiently large.
CALL check_floquet_conversion(qs_env, bs_env, cfm_eigenvectors)
! Write the quasi-energies ε_α folded into the first Floquet Brillouin
! zone -ħΩ/2 < ε_α ≤ ħΩ/2 for this k-point and spin channel
CALL write_quasi_energy(qs_env, bs_env, ispin, ikp_for_file, xkp, eigenvalues)
! Compute the physical spectral function
! A_phys(k,ω) = -(1/π) Im[ Tr_KS G^R_{00}(k,ω) ]
! from the zero-zero Floquet sector of the retarded Green's function
CALL calculate_floquet_spectral_function(qs_env, bs_env, floquet_matrix, a_k)
! Write A_phys(k,ω) (or f(ω,μ)·A_phys for the occupied spectral weight)
! to output for this k-point and spin channel
CALL write_floquet_spectral_function(qs_env, bs_env, ispin, ikp_for_file, xkp, a_k)
! Write progress in the main output file
IF (unit_nr > 0) WRITE (unit_nr, '(2X,I6,T12,I4,T18,3F12.6)') ikp_for_file, ispin, xkp(1:3)
END DO
END DO
IF (unit_nr > 0) WRITE (unit_nr, '(A)') ""
DEALLOCATE (eigenvalues)
CALL cp_cfm_release(floquet_matrix)
CALL cp_cfm_release(floquet_copy)
CALL cp_cfm_release(cfm_eigenvectors)
CALL cp_fm_struct_release(floquet_struct)
CALL timestop(handle)
END SUBROUTINE floquet
END MODULE floquet_main

656
src/floquet_utils.F Normal file
View file

@ -0,0 +1,656 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Floquet stuff
!> \par History
!> \author Shridhar Shanbhag (27.01.2026)
! **************************************************************************************************
MODULE floquet_utils
USE cell_types, ONLY: cell_type,&
get_cell
USE cp_cfm_basic_linalg, ONLY: cp_cfm_lu_invert,&
cp_cfm_scale_and_add
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_get_submatrix,&
cp_cfm_release,&
cp_cfm_set_all,&
cp_cfm_set_submatrix,&
cp_cfm_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: dbcsr_p_type
USE cp_files, ONLY: close_file,&
open_file
USE kinds, ONLY: default_string_length,&
dp
USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices,&
rs_to_kp
USE kpoint_methods, ONLY: kpoint_init_cell_index
USE kpoint_types, ONLY: get_kpoint_info,&
kpoint_create,&
kpoint_release,&
kpoint_type
USE mathconstants, ONLY: gaussi,&
pi,&
z_one,&
z_zero
USE mathlib, ONLY: geeig_right,&
gemm_square
USE message_passing, ONLY: mp_para_env_type
USE physcon, ONLY: evolt
USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
USE qs_moments, ONLY: qs_moment_kpoints_deep
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE util, ONLY: sort
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'floquet_utils'
! Public subroutines
PUBLIC :: build_floquet_matrix, &
calculate_epsilon_derivative, &
build_momentum_matrix, &
check_floquet_conversion, &
calculate_floquet_spectral_function, &
write_quasi_energy, &
write_floquet_spectral_function
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param xkp ...
!> \param e_k ...
!> \param de_dk ...
! **************************************************************************************************
SUBROUTINE calculate_epsilon_derivative(qs_env, xkp, e_k, de_dk)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
OPTIONAL :: e_k
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :, :, :), OPTIONAL :: de_dk
CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_epsilon_derivative'
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: C_dH_C, C_dS_C, C_k, dH_dk_i, dS_dk_i, &
H_k, S_k
INTEGER :: handle, i_dir, ikp, ispin, n, n_img_all, &
n_spin, nao, nkp
INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
LOGICAL :: present_dedk, present_ek
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvals
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: H_rs, S_rs
REAL(KIND=dp), DIMENSION(3, 3) :: hmat
TYPE(cell_type), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_all
CALL timeset(routineN, handle)
present_ek = PRESENT(e_k) ! calculate band energies, ε_k for all kpoints xkp
present_dedk = PRESENT(de_dk) ! calculate derivative, ∇_k ε_k of band energies for all kpoints xkp
IF (.NOT. (present_ek .OR. present_dedk)) CPABORT("Subroutine needs either e_k or de_dk")
CALL get_qs_env(qs_env, &
matrix_ks_kp=matrix_ks_kp, &
matrix_s_kp=matrix_s_kp, &
sab_all=sab_all, &
cell=cell, &
kpoints=kpoints_scf, &
para_env=para_env, &
mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
CALL get_cell(cell=cell, h=hmat)
n_spin = SIZE(matrix_ks_kp, 1)
nkp = SIZE(xkp, 2)
! create kpoint environment kpoints_all which contains all neighbor cells R
! without considering any lattice symmetry
NULLIFY (kpoints_all)
CALL kpoint_create(kpoints_all)
CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_all)
CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index_all, index_to_cell=index_to_cell_all)
ALLOCATE (S_rs(1, nao, nao, n_img_all), H_rs(n_spin, nao, nao, n_img_all), source=0.0_dp)
! Convert real-space dbcsr matrices into arrays
CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, S_rs, cell_to_index_all)
CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, H_rs, cell_to_index_all)
IF (present_dedk) ALLOCATE (de_dk(n_spin, nkp, 3, nao), source=0.0_dp)
IF (present_ek) ALLOCATE (e_k(n_spin, nkp, nao), source=0.0_dp)
!$OMP PARALLEL PRIVATE(ikp, ispin, S_k, H_k, &
!$OMP eigenvals, C_k, dS_dk_i, dH_dk_i, C_dS_C, C_dH_C) &
!$OMP SHARED(nao, n_spin, de_dk, e_k, present_ek, present_dedk, &
!$OMP nkp, xkp, S_rs, H_rs, index_to_cell_all, hmat)
IF (present_dedk) ALLOCATE (dS_dk_i(nao, nao), C_dS_C(nao, nao), &
dH_dk_i(nao, nao), C_dH_C(nao, nao), source=z_zero)
ALLOCATE (C_k(nao, nao), S_k(nao, nao), H_k(nao, nao), source=z_zero)
ALLOCATE (eigenvals(nao), source=0.0_dp)
!$OMP DO COLLAPSE(2)
DO ispin = 1, n_spin
DO ikp = 1, nkp
! S^R -> S(k), H^R -> H(k)
S_k = 0
H_k = 0
CALL rs_to_kp(S_rs(1, :, :, :), S_k, index_to_cell_all, xkp(:, ikp))
CALL rs_to_kp(H_rs(ispin, :, :, :), H_k, index_to_cell_all, xkp(:, ikp))
! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
CALL geeig_right(H_k, S_k, eigenvals, C_k)
IF (present_ek) e_k(ispin, ikp, :) = eigenvals(:)
IF (present_dedk) THEN
! Evaluate the derivatives using
! ∇ ε_k = C^H(k) ∇ H_k C(k) - ε_k C^H(k) ∇ S_k C(k)
DO i_dir = 1, 3
CALL rs_to_kp(S_rs(1, :, :, :), dS_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
CALL rs_to_kp(H_rs(ispin, :, :, :), dH_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
CALL gemm_square(C_k, 'C', dS_dk_i, 'N', C_k, 'N', C_dS_C)
CALL gemm_square(C_k, 'C', dH_dk_i, 'N', C_k, 'N', C_dH_C)
DO n = 1, nao
de_dk(ispin, ikp, i_dir, n) = DBLE(C_dH_C(n, n)) - DBLE(eigenvals(n)*C_dS_C(n, n))
END DO
END DO
END IF
END DO
END DO
!$OMP END DO
IF (present_dedk) DEALLOCATE (dS_dk_i, C_dS_C, dH_dk_i, C_dH_C)
DEALLOCATE (S_k, H_k, C_k, eigenvals)
!$OMP END PARALLEL
DEALLOCATE (S_rs, H_rs)
CALL kpoint_release(kpoints_all)
CALL timestop(handle)
END SUBROUTINE calculate_epsilon_derivative
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param xkp ...
!> \param momentum ...
! **************************************************************************************************
SUBROUTINE build_momentum_matrix(qs_env, xkp, momentum)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), DIMENSION(3) :: xkp
COMPLEX(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :, :, :) :: momentum
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_momentum_matrix'
COMPLEX(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :, :, :, :) :: dipole
INTEGER :: handle, i, i_dir, ispin, j, n_spin, nao
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp1
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: e_k
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: de_dk
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
! We calculate momentum matrix elements p_nm = <ψ_n|-iħ∇_r|ψ_m>
! p_nm = <u_n|(ħk - iħ∇_r)|u_m>
! p_nm = i d_nm (ε_n - ε_m) + ∇_k ε_n δ_nm
NULLIFY (dft_control)
CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
ALLOCATE (xkp1(3, 1), source=0.0_dp)
xkp1(:, 1) = xkp(:)
CALL calculate_epsilon_derivative(qs_env, xkp1, e_k, de_dk)
n_spin = dft_control%nspins
CALL qs_moment_kpoints_deep(qs_env, xkp1, dipole)
!$OMP PARALLEL PRIVATE(ispin, i_dir, i, j) &
!$OMP SHARED(nao, n_spin, momentum, dipole, e_k, de_dk)
!$OMP DO COLLAPSE(4)
DO ispin = 1, n_spin
DO i_dir = 1, 3
DO i = 1, nao
DO j = 1, nao
IF (j == i) THEN
momentum(ispin, i_dir, i, j) = de_dk(ispin, 1, i_dir, i)
ELSE
momentum(ispin, i_dir, i, j) = gaussi*(e_k(ispin, 1, i) - e_k(ispin, 1, j))* &
dipole(ispin, 1, i_dir, i, j)
END IF
END DO
END DO
END DO
END DO
!$OMP END DO
!$OMP END PARALLEL
DEALLOCATE (xkp1, dipole, e_k, de_dk)
CALL timestop(handle)
END SUBROUTINE build_momentum_matrix
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param xkp ...
!> \param ispin ...
!> \param off_diag_m ...
! **************************************************************************************************
SUBROUTINE build_off_diagonal_matrix(qs_env, bs_env, xkp, ispin, off_diag_m)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
REAL(KIND=dp), DIMENSION(3) :: xkp
INTEGER :: ispin
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: off_diag_m
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_off_diagonal_matrix'
COMPLEX(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :, :, :) :: momentum
COMPLEX(KIND=dp), DIMENSION(3) :: efactor
INTEGER :: handle, i, i_dir, j, n_spin, nao
REAL(KIND=dp) :: amplitude, omega
REAL(KIND=dp), DIMENSION(3) :: e_vec, phi, polarisation
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
! Builds the matrix that occupies the off diagonal blocks in the Floquet Matrix H_F
polarisation(:) = bs_env%floquet_polarisation(:)
IF (SQRT(SUM(polarisation**2)) < EPSILON(0.0_dp)) THEN
CPABORT("Invalid (too small) polarisation vector specified for POLARISATION_PUMP")
END IF
amplitude = bs_env%floquet_amplitude
e_vec(:) = amplitude*polarisation
phi(:) = pi*bs_env%floquet_phi(:)
omega = bs_env%floquet_omega
n_spin = bs_env%n_spin
CALL get_qs_env(qs_env, mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
ALLOCATE (momentum(n_spin, 3, nao, nao), source=z_zero)
CALL build_momentum_matrix(qs_env, xkp, momentum)
! E_factor(α) = (i E(α) exp(i·φ(α)))/(2ω) where α = x,y,z
DO i_dir = 1, 3
efactor(i_dir) = gaussi*e_vec(i_dir)*CMPLX(DCOS(phi(i_dir)), DSIN(phi(i_dir)), KIND=dp)/(2*omega)
END DO
DO i_dir = 1, 3
DO i = 1, nao
DO j = 1, nao
off_diag_m(i, j) = off_diag_m(i, j) + &
momentum(ispin, i_dir, i, j)*efactor(i_dir)
END DO
END DO
END DO
DEALLOCATE (momentum)
CALL timestop(handle)
END SUBROUTINE build_off_diagonal_matrix
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param xkp ...
!> \param ispin ...
!> \param f_index ...
!> \param diag_e ...
! **************************************************************************************************
SUBROUTINE build_diagonal_matrix(qs_env, bs_env, xkp, ispin, f_index, diag_e)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
REAL(KIND=dp), DIMENSION(3) :: xkp
INTEGER :: ispin, f_index
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: diag_e
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_diagonal_matrix'
INTEGER :: handle, i, n_spin, nao
REAL(KIND=dp) :: omega
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp1
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: e_k
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
! Builds the matrix that occupies the off diagonal blocks in the Floquet Matrix H_F
omega = bs_env%floquet_omega
n_spin = bs_env%n_spin
CALL get_qs_env(qs_env, mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
ALLOCATE (xkp1(3, 1), source=0.0_dp)
xkp1(:, 1) = xkp(:)
CALL calculate_epsilon_derivative(qs_env, xkp1, e_k)
DO i = 1, nao
diag_e(i, i) = e_k(ispin, 1, i) + f_index*omega
END DO
DEALLOCATE (xkp1, e_k)
CALL timestop(handle)
END SUBROUTINE build_diagonal_matrix
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param xkp ...
!> \param ispin ...
!> \param floquet_matrix ...
! **************************************************************************************************
SUBROUTINE build_floquet_matrix(qs_env, bs_env, xkp, ispin, floquet_matrix)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
REAL(KIND=dp), DIMENSION(3) :: xkp
INTEGER :: ispin
TYPE(cp_cfm_type) :: floquet_matrix
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_floquet_matrix'
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: conj_off_diag_m, diag_e, off_diag_m
INTEGER :: f_index, handle, i, i_f, max_f_index, &
n_fbands, n_spin, nao
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp1
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, &
mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
max_f_index = bs_env%max_floquet_index
n_fbands = 1 + 2*max_f_index
n_spin = bs_env%n_spin
! Creates the floquet matrix H_F by placing the diagonal and off diagonal blocks
ALLOCATE (xkp1(3, 1), source=0.0_dp)
ALLOCATE (diag_e(nao, nao), source=z_zero)
ALLOCATE (off_diag_m(nao, nao), source=z_zero)
ALLOCATE (conj_off_diag_m(nao, nao), source=z_zero)
xkp1(:, 1) = xkp(:)
CALL build_off_diagonal_matrix(qs_env, bs_env, xkp, ispin, off_diag_m)
conj_off_diag_m(:, :) = CONJG(TRANSPOSE(off_diag_m(:, :)))
floquet_matrix%local_data(:, :) = z_zero
DO i = 1, n_fbands
i_f = 1 + (i - 1)*nao
f_index = i - max_f_index - 1
CALL build_diagonal_matrix(qs_env, bs_env, xkp, ispin, f_index, diag_e)
CALL cp_cfm_set_submatrix(floquet_matrix, diag_e, i_f, i_f)
IF (i > 1) THEN
CALL cp_cfm_set_submatrix(floquet_matrix, off_diag_m, i_f, i_f - nao)
CALL cp_cfm_set_submatrix(floquet_matrix, conj_off_diag_m, i_f - nao, i_f)
END IF
END DO
DEALLOCATE (diag_e, off_diag_m, conj_off_diag_m, xkp1)
CALL timestop(handle)
END SUBROUTINE build_floquet_matrix
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param cfm_eigenvectors ...
! **************************************************************************************************
SUBROUTINE check_floquet_conversion(qs_env, bs_env, cfm_eigenvectors)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
TYPE(cp_cfm_type) :: cfm_eigenvectors
CHARACTER(LEN=*), PARAMETER :: routineN = 'check_floquet_conversion'
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vector_block_0, vector_block_1
INTEGER :: handle, max_f_index, nao, start_row_0, &
start_row_1
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
max_f_index = bs_env%max_floquet_index
start_row_0 = 1 + nao*max_f_index
start_row_1 = 1 + nao*(max_f_index + 1)
ALLOCATE (vector_block_0(3*nao, nao), vector_block_1(3*nao, nao), source=z_zero)
CALL cp_cfm_get_submatrix(cfm_eigenvectors, vector_block_0, start_row_0, start_row_0)
CALL cp_cfm_get_submatrix(cfm_eigenvectors, vector_block_1, start_row_1, start_row_1)
IF (bs_env%eps_floquet > 0) THEN
IF (ABS(MAXVAL(ABS(vector_block_0)) - MAXVAL(ABS(vector_block_1))) > bs_env%eps_floquet) THEN
CPABORT("FLOQUET DIAGONALIZATION DID NOT CONVERGE. MAX_FLOQUET_INDEX TOO SMALL")
END IF
END IF
DEALLOCATE (vector_block_0, vector_block_1)
CALL timestop(handle)
END SUBROUTINE check_floquet_conversion
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param floquet_matrix ...
!> \param a_k ...
! **************************************************************************************************
SUBROUTINE calculate_floquet_spectral_function(qs_env, bs_env, floquet_matrix, a_k)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
TYPE(cp_cfm_type) :: floquet_matrix
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a_k
CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_floquet_spectral_function'
COMPLEX(KIND=dp) :: diag
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g00
INTEGER :: handle, i, i_E, max_f_index, n_E, nao, &
start_row_g00
REAL(KIND=dp) :: broad, E_min, energy, energy_step, mu
TYPE(cp_cfm_type) :: g_r
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, mos=mos)
CALL get_mo_set(mo_set=mos(1), mu=mu, nao=nao)
broad = bs_env%broadening_floquet
energy_step = bs_env%energy_step_floquet
E_min = mu - bs_env%energy_window_floquet
max_f_index = bs_env%max_floquet_index
n_E = SIZE(a_k)
start_row_g00 = 1 + nao*max_f_index
ALLOCATE (g00(nao, nao))
CALL cp_cfm_create(g_r, floquet_matrix%matrix_struct, set_zero=.TRUE.)
DO i_E = 1, n_E
energy = E_min + i_E*energy_step
diag = energy + gaussi*broad/2.0_dp
CALL cp_cfm_set_all(g_r, z_zero, diag)
CALL cp_cfm_scale_and_add(z_one, g_r, -z_one, floquet_matrix)
CALL cp_cfm_lu_invert(g_r)
g00(:, :) = 0.0_dp
CALL cp_cfm_get_submatrix(g_r, g00, start_row_g00, start_row_g00)
a_k(i_E) = -AIMAG(SUM([(g00(i, i), i=1, nao)]))/pi
END DO
DEALLOCATE (g00)
CALL cp_cfm_release(g_r)
CALL timestop(handle)
END SUBROUTINE calculate_floquet_spectral_function
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param ispin ...
!> \param ikp_for_file ...
!> \param xkp ...
!> \param eigenvalues ...
! **************************************************************************************************
SUBROUTINE write_quasi_energy(qs_env, bs_env, ispin, ikp_for_file, xkp, eigenvalues)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
INTEGER :: ispin, ikp_for_file
REAL(KIND=dp), DIMENSION(3) :: xkp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
CHARACTER(LEN=*), PARAMETER :: routineN = 'write_quasi_energy'
CHARACTER(LEN=default_string_length) :: fname
INTEGER :: handle, i, j, n, nao, nkp_start, qunit
INTEGER, ALLOCATABLE, DIMENSION(:) :: indices, work
REAL(KIND=dp) :: fbz_max, fbz_min, omega, qe
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: quasi_energies
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, mos=mos)
CALL get_mo_set(mo_set=mos(1), nao=nao)
nkp_start = bs_env%nkp_only_DOS
omega = bs_env%floquet_omega
fbz_max = omega/2.0_dp
fbz_min = -omega/2.0_dp
n = SIZE(eigenvalues)
ALLOCATE (indices(nao), work(nao), quasi_energies(nao))
indices(:) = [(i, i=1, nao)]
DO i = 1, nao
DO j = 1, n
IF (ANY(indices == j)) CYCLE
IF (ABS(eigenvalues(j)) < ABS(eigenvalues(indices(i)))) indices(i) = j
END DO
END DO
CALL sort(indices, nao, work)
quasi_energies(:) = eigenvalues(indices(:))
IF (bs_env%para_env%is_source()) THEN
WRITE (fname, "(2A)") TRIM(bs_env%floquet_qe_file), ".bs"
IF (ikp_for_file == 1) THEN
CALL open_file(TRIM(fname), unit_number=qunit, file_status="REPLACE", &
file_action="WRITE")
ELSE
CALL open_file(TRIM(fname), unit_number=qunit, file_status="OLD", &
file_action="WRITE", file_position="APPEND")
END IF
WRITE (qunit, "(A)") "# Quasi-energies obtained by diagonalising the Floquet Hamiltonian"
WRITE (qunit, "(A)") "# (in units of eV, folded to Floquet Brillouin zone)"
WRITE (qunit, "(A,I0,T10,A,I0,A,T24,3(1X,F14.8))") &
"# Spin ", ispin, " Point ", ikp_for_file, ": ", xkp(1:3)
WRITE (qunit, "(A)") "# Floquet band Quasi-energy [eV]"
DO i = 1, nao
qe = quasi_energies(i)
WRITE (qunit, "(I8,F21.8)") i, qe*evolt
END DO
CALL close_file(qunit)
END IF
DEALLOCATE (indices, work, quasi_energies)
CALL timestop(handle)
END SUBROUTINE write_quasi_energy
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param ispin ...
!> \param ikp_for_file ...
!> \param xkp ...
!> \param a_k ...
! **************************************************************************************************
SUBROUTINE write_floquet_spectral_function(qs_env, bs_env, ispin, ikp_for_file, xkp, a_k)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
INTEGER :: ispin, ikp_for_file
REAL(KIND=dp), DIMENSION(3) :: xkp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a_k
CHARACTER(LEN=*), PARAMETER :: routineN = 'write_floquet_spectral_function'
CHARACTER(LEN=default_string_length) :: fname
INTEGER :: handle, i_E, n_E, nkp_start, wunit
REAL(KIND=dp) :: E_min, energy, energy_step, mu
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, mos=mos)
CALL get_mo_set(mo_set=mos(1), mu=mu)
energy_step = bs_env%energy_step_floquet
E_min = mu - bs_env%energy_window_floquet
n_E = SIZE(a_k)
nkp_start = bs_env%nkp_only_DOS
IF (bs_env%para_env%is_source()) THEN
WRITE (fname, "(2A)") TRIM(bs_env%floquet_dos_file), ".out"
IF (ikp_for_file == 1) THEN
CALL open_file(TRIM(fname), unit_number=wunit, file_status="REPLACE", &
file_action="WRITE")
ELSE
CALL open_file(TRIM(fname), unit_number=wunit, file_status="OLD", &
file_action="WRITE", file_position="APPEND")
END IF
WRITE (wunit, "(A)") "# Floquet Density of States: D(ω,k) = -1/π*Im[Tr_KS(G^R(ω,k))]"
WRITE (wunit, "(A,I0,T10,A,I0,A,T24,3(1X,F14.8))") &
"# Spin ", ispin, " Point ", ikp_for_file, ": ", xkp(1:3)
WRITE (wunit, "(A)") "#Energy-E_F (eV) A(ω,k) = DOS (1/eV)"
DO i_E = 1, n_E
energy = E_min + i_E*energy_step
WRITE (wunit, "(2X,2G13.4)") energy*evolt, a_k(i_E)/evolt
END DO
CALL close_file(wunit)
END IF
CALL timestop(handle)
END SUBROUTINE write_floquet_spectral_function
END MODULE floquet_utils

125
src/input_cp2k_properties_dft.F
View file

@ -2374,6 +2374,10 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_floquet_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END SUBROUTINE create_bandstructure_section
! **************************************************************************************************
@ -2899,4 +2903,125 @@ CONTAINS
END SUBROUTINE create_tipscan_section
! **************************************************************************************************
!> \brief ...
!> \param section ...
!> \param section_name ...
!> \author Shridhar Shanbhag
! **************************************************************************************************
SUBROUTINE create_floquet_section(section)
TYPE(section_type), POINTER :: section
TYPE(keyword_type), POINTER :: keyword
CALL section_create(section, __LOCATION__, name="FLOQUET", &
description="Parameters controlling the calculation of the "// &
"Floquet band structure and Quasi-energies of a system driven "// &
"by a periodic monochromatic electric field with any arbitrary "// &
"polarisation state.", &
n_keywords=1, n_subsections=0, repeats=.FALSE.)
NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="AMPLITUDE", &
description="Electric field amplitude of the monochromatic light.", &
usage="AMPLITUDE [Vm-1] 1.0E9", n_var=1, type_of_var=real_t, &
unit_str="Vm-1", default_r_val=0.0_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY", &
description="Frequency of the electric field for Floquet calculations "// &
"expressed in terms of the energy of a photon in eV.", &
usage="FREQUENCY 1.5", unit_str="eV", n_var=1, type_of_var=real_t, &
default_r_val=cp_unit_to_cp2k(value=1.0_dp, unit_str="eV"))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="POLARISATION", &
description="Polarisation vector of the input light. The "// &
"amplitude in each direction is the product of the AMPLITUDE "// &
"and the POLARISATION vector component. Eg. for light polarized "// &
"in the x direction, use 1 0 0. For left-circularly polarized "// &
"light travelling in the +z direction use 1 1 0 and set "// &
"PHASE_OFFSET 0 0.5 0 for a +π/2 offset", &
usage="POLARISATION 0.0 0.0 1.0", &
default_r_vals=[0.0_dp, 0.0_dp, 1.0_dp])
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="PHASE_OFFSETS", &
description="Phase offset of the electric field in the "// &
"x, y, and z directions given in multiples of π. For "// &
"linearly polarized light, use 0 0 0 and for left-circularly "// &
"polarized light travelling in the +z direction use "// &
"PHASE_OFFSET 0.0 0.5 0 ", &
usage="PHASE_OFFSETS 0.5 0 0", &
default_r_vals=[0.0_dp, 0.0_dp, 0.0_dp])
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAX_FLOQUET_INDEX", &
description="Largest absolute Floquet index up to which "// &
"Floquet Hamiltonian is truncated. Use a larger value to "// &
"ensure convergence at higher computational cost and a "// &
"smaller value for faster computation. EPS_FLOQUET checks "// &
"to ensure that the error due to truncation is small.", &
usage="MAX_FLOQUET_INDEX 50", &
default_i_val=50)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_FLOQUET", &
description="Threshold on the error due to truncation of "// &
"the Floquet Hamiltonian. Used to verify that MAX_FLOQUET_INDEX "// &
"is large enough that truncation errors are small. "// &
"If negative, the check is removed (not recommended).", &
usage="EPS_FLOQUET 1.e-10", default_r_val=1.0e-10_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_WINDOW", &
description="Half-width of the energy range, in eV, centred "// &
"on the Fermi level, over which the Floquet spectral "// &
"function is evaluated.", &
usage="ENERGY_WINDOW 5.0", &
default_r_val=cp_unit_to_cp2k(value=10.0_dp, unit_str="eV"), &
unit_str="eV")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_STEP", &
description="Resolution of the energy E used to compute the "// &
"spectral function within the energy window.", &
usage="ENERGY_STEP 0.01", &
default_r_val=cp_unit_to_cp2k(value=0.01_dp, unit_str="eV"), &
unit_str="eV")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="BROADENING", &
description="Lorentzian broadening applied to the peaks of the "// &
"Floquet spectral function.", &
usage="BROADENING 0.02", &
default_r_val=cp_unit_to_cp2k(value=0.02_dp, unit_str="eV"), &
unit_str="eV")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FLOQUET_DOS_FILE_NAME", &
description="File name used for the Floquet spectral function.", &
usage="FLOQUET_DOS_FILE_NAME FLOQUET_DOS", &
default_lc_val="FLOQUET_DOS")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="QUASI_ENERGIES_FILE_NAME", &
description="File name used for the Floquet quasi-energies.", &
usage="QUASI_ENERGIES_FILE_NAME QUASI_ENERGIES", &
default_lc_val="QUASI_ENERGIES")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_floquet_section
END MODULE input_cp2k_properties_dft

6
src/post_scf_bandstructure_methods.F
View file

@ -6,6 +6,7 @@
!--------------------------------------------------------------------------------------------------!
MODULE post_scf_bandstructure_methods
USE floquet_main, ONLY: floquet
USE gw_main, ONLY: gw
USE input_section_types, ONLY: section_vals_type
USE post_scf_bandstructure_utils, ONLY: create_and_init_bs_env,&
@ -58,6 +59,11 @@ CONTAINS
! density of states (DOS), projected DOS, local DOS for DFT, DFT+SOC, G0W0, G0W0+SOC
CALL eval_bandstructure_properties(qs_env, qs_env%bs_env)
! Floquet Hamiltonian Diagonalization and calculation of Floquet Density of States
IF (qs_env%bs_env%do_floquet) THEN
CALL floquet(qs_env, qs_env%bs_env)
END IF
CALL timestop(handle)
END SUBROUTINE post_scf_bandstructure

17
src/post_scf_bandstructure_types.F
View file

@ -25,6 +25,7 @@ MODULE post_scf_bandstructure_types
USE input_constants, ONLY: rtp_method_bse,&
small_cell_full_kp
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE kpoint_types, ONLY: kpoint_release,&
kpoint_type
@ -88,7 +89,8 @@ MODULE post_scf_bandstructure_types
do_soc = .FALSE., &
do_ldos = .FALSE., &
do_gw_ri_rs = .FALSE., &
do_dos_pdos = .FALSE.
do_dos_pdos = .FALSE., &
do_floquet = .FALSE.
! various eigenvalues computed in GW code, some depend on k-points
! and have therefore three dimensions (band index, k-point, spin)
@ -310,6 +312,19 @@ MODULE post_scf_bandstructure_types
INTEGER :: n_bins_max_for_printing = -1
REAL(KIND=dp) :: unit_ldos_int_z_inv_Ang2_eV = -1.0_dp
! parameters for Floquet band structure calculations
INTEGER :: max_floquet_index = -1
REAL(KIND=dp), DIMENSION(:), POINTER :: floquet_polarisation => NULL(), &
floquet_phi => NULL()
REAL(KIND=dp) :: floquet_omega = -1.0_dp, &
floquet_amplitude = -1.0_dp, &
eps_floquet = -1.0_dp, &
broadening_floquet = -1.0_dp, &
energy_step_floquet = -1.0_dp, &
energy_window_floquet = -1.0_dp
CHARACTER(LEN=default_string_length) :: floquet_dos_file = "", &
floquet_qe_file = ""
! quantities only needed for small cells and k-point sampling in DFT (small_cell_full_kp)
INTEGER :: nkp_scf_desymm = -1, &
nimages_3c = -1, &

19
src/post_scf_bandstructure_utils.F
View file

@ -197,8 +197,8 @@ CONTAINS
REAL(KIND=dp), DIMENSION(3) :: kpptr
REAL(KIND=dp), DIMENSION(3, 3) :: cart_hmat
TYPE(cell_type), POINTER :: cell
TYPE(section_vals_type), POINTER :: dos_pdos_sec, gw_ri_rs_sec, gw_sec, &
kp_bs_sec, ldos_sec, soc_sec
TYPE(section_vals_type), POINTER :: dos_pdos_sec, floquet_sec, gw_ri_rs_sec, &
gw_sec, kp_bs_sec, ldos_sec, soc_sec
CALL timeset(routineN, handle)
NULLIFY (cell)
@ -247,6 +247,21 @@ CONTAINS
CALL uppercase(ustr)
CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", n_rep_val=bs_env%input_kp_bs_n_sp_pts)
NULLIFY (floquet_sec)
floquet_sec => section_vals_get_subs_vals(bs_sec, "FLOQUET")
CALL section_vals_get(floquet_sec, explicit=bs_env%do_floquet)
CALL section_vals_val_get(floquet_sec, "AMPLITUDE", r_val=bs_env%floquet_amplitude)
CALL section_vals_val_get(floquet_sec, "FREQUENCY", r_val=bs_env%floquet_omega)
CALL section_vals_val_get(floquet_sec, "POLARISATION", r_vals=bs_env%floquet_polarisation)
CALL section_vals_val_get(floquet_sec, "PHASE_OFFSETS", r_vals=bs_env%floquet_phi)
CALL section_vals_val_get(floquet_sec, "MAX_FLOQUET_INDEX", i_val=bs_env%max_floquet_index)
CALL section_vals_val_get(floquet_sec, "EPS_FLOQUET", r_val=bs_env%eps_floquet)
CALL section_vals_val_get(floquet_sec, "ENERGY_WINDOW", r_val=bs_env%energy_window_floquet)
CALL section_vals_val_get(floquet_sec, "ENERGY_STEP", r_val=bs_env%energy_step_floquet)
CALL section_vals_val_get(floquet_sec, "BROADENING", r_val=bs_env%broadening_floquet)
CALL section_vals_val_get(floquet_sec, "FLOQUET_DOS_FILE_NAME", c_val=bs_env%floquet_dos_file)
CALL section_vals_val_get(floquet_sec, "QUASI_ENERGIES_FILE_NAME", c_val=bs_env%floquet_qe_file)
! read special points for band structure
ALLOCATE (bs_env%xkp_special(3, bs_env%input_kp_bs_n_sp_pts))
DO ikp = 1, bs_env%input_kp_bs_n_sp_pts

87
tests/QS/regtest-floquet/Floquet_KP.inp Normal file
View file

@ -0,0 +1,87 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT Floquet_KP
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
&KPOINTS
FULL_GRID T
SCHEME MONKHORST-PACK 4 4 1
&END KPOINTS
&POISSON
PERIODIC XYZ
&END POISSON
&QS
EXTRAPOLATION USE_GUESS ! required for K-Point sampling
!EPS_DEFAULT 1.0E-30
&END QS
&SCF
ADDED_MOS -1
EPS_SCF 1.0E-6
MAX_SCF 10000
SCF_GUESS ATOMIC
&DIAGONALIZATION
ALGORITHM STANDARD
EPS_ADAPT 0.01
&END DIAGONALIZATION
&MIXING
ALPHA 0.2
BETA 1.5
METHOD BROYDEN_MIXING
NBROYDEN 8
&END MIXING
&SMEAR ON
ELECTRONIC_TEMPERATURE [K] 300
METHOD FERMI_DIRAC
&END SMEAR
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&BANDSTRUCTURE
&BANDSTRUCTURE_PATH
NPOINTS 1
SPECIAL_POINT K 1/3 2/3 0.0
&END BANDSTRUCTURE_PATH
&FLOQUET
AMPLITUDE 1.06E9
BROADENING 0.02
ENERGY_STEP 1.69
ENERGY_WINDOW 3.38
FREQUENCY 1.5
MAX_FLOQUET_INDEX 50
PHASE_OFFSETS 0.0 0.5 0.0
POLARISATION 1 1 0
&END FLOQUET
&END BANDSTRUCTURE
&END PROPERTIES
&SUBSYS
&CELL
ABC [angstrom] 2.4612 2.4612 16
ALPHA_BETA_GAMMA 90. 90. 60.
PERIODIC XY
&END CELL
&COORD
SCALED
C 1./3. 1./3. 0.
C 2./3. 2./3. 0.
&END COORD
&KIND C
BASIS_SET SZV-MOLOPT-GTH
ELEMENT C
POTENTIAL GTH-PBE
&END KIND
&TOPOLOGY
&CENTER_COORDINATES T
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

11
tests/QS/regtest-floquet/TEST_FILES.toml Normal file
View file

@ -0,0 +1,11 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
#
# compute ground state for graphene with PBE and SZV
# compute the floquet DOS and QUASI_ENERGIES at K-point
"Floquet_KP.inp" = [{matcher="Quasienergy", file="QUASI_ENERGIES.bs", tol=1.0E-06, ref=0.00534452},
{matcher="Floquet_DOS", file="FLOQUET_DOS.out", tol=1e-8, ref=24.91}]
#EOF

1
tests/TEST_DIRS
View file

@ -3,6 +3,7 @@
# Directories have been reordered according the execution time needed for a gfortran pdbg run using 2 MPI tasks
# in case a new directory is added just add it at the top of the list..
# the order will be regularly checked and modified...
QS/regtest-floquet
QS/regtest-dft-d4-auto-ref libdftd4
QS/regtest-dft-d2-auto-ref
QS/regtest-dft-d3-auto-ref libint

4
tests/matchers.py
View file

@ -385,4 +385,8 @@ registry["RIRS_CUTOFF"] = GenericMatcher(
)
registry["E_RIRS_HOMO"] = GenericMatcher(r"G0W0 valence band maximum", col=6)
registry["E_RIRS_LUMO"] = GenericMatcher(r"G0W0 conduction band minimum", col=6)
# Floquet Calculations
registry["Quasienergy"] = GenericMatcher(r" 4", col=2)
registry["Floquet_DOS"] = GenericMatcher(r"-1.690", col=2)
# EOF