diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 95e4915bcd..b20da64d58 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/floquet_main.F b/src/floquet_main.F new file mode 100644 index 0000000000..72bfee0b89 --- /dev/null +++ b/src/floquet_main.F @@ -0,0 +1,215 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! 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 diff --git a/src/floquet_utils.F b/src/floquet_utils.F new file mode 100644 index 0000000000..f2b4149a35 --- /dev/null +++ b/src/floquet_utils.F @@ -0,0 +1,656 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! 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 = + ! 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 diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index 4ae3ade219..2665f3e6d9 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -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 diff --git a/src/post_scf_bandstructure_methods.F b/src/post_scf_bandstructure_methods.F index bf104630b8..f56cecc82a 100644 --- a/src/post_scf_bandstructure_methods.F +++ b/src/post_scf_bandstructure_methods.F @@ -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 diff --git a/src/post_scf_bandstructure_types.F b/src/post_scf_bandstructure_types.F index c7689e9169..2bd4a4ec97 100644 --- a/src/post_scf_bandstructure_types.F +++ b/src/post_scf_bandstructure_types.F @@ -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, & diff --git a/src/post_scf_bandstructure_utils.F b/src/post_scf_bandstructure_utils.F index 89540cc929..5f493d4d9f 100644 --- a/src/post_scf_bandstructure_utils.F +++ b/src/post_scf_bandstructure_utils.F @@ -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 diff --git a/tests/QS/regtest-floquet/Floquet_KP.inp b/tests/QS/regtest-floquet/Floquet_KP.inp new file mode 100644 index 0000000000..5970a26522 --- /dev/null +++ b/tests/QS/regtest-floquet/Floquet_KP.inp @@ -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 diff --git a/tests/QS/regtest-floquet/TEST_FILES.toml b/tests/QS/regtest-floquet/TEST_FILES.toml new file mode 100644 index 0000000000..c0c611eb5d --- /dev/null +++ b/tests/QS/regtest-floquet/TEST_FILES.toml @@ -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 diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 416ab19558..419b110645 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -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 diff --git a/tests/matchers.py b/tests/matchers.py index 11b025c82c..669dd3b686 100644 --- a/tests/matchers.py +++ b/tests/matchers.py @@ -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