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

Resonant Inelastic X-ray Spectroscopy (RIXS) Module

Browse source
This commit is contained in:
BelizSertcan 2025-07-22 16:39:38 +02:00 • committed by GitHub
parent 856d28c997
commit bf603df4a0
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
25 changed files with 1571 additions and 62 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/CMakeLists.txt
View file

@ -741,6 +741,8 @@ list(
response_solver.F
restraint.F
ri_environment_methods.F
rixs_methods.F
rixs_types.F
rmsd.F
rpa_exchange.F
rpa_communication.F

9
src/common/bibliography.F
View file

@ -89,7 +89,8 @@ MODULE bibliography
Stein2022, Stein2024, Pracht2019, &
Blase2018, Blase2020, Bruneval2015, Golze2019, Gui2018, Jacquemin2017, Liu2020, &
Sander2015, Schreiber2008, vanSetten2015, Setyawan2010, Ahart2024, Knysh2024, &
Schambeck2024, Mewes2018, Sertcan2024, Drautz2019, Lysogorskiy2021, Bochkarev2024
Schambeck2024, Mewes2018, Sertcan2024, Drautz2019, Lysogorskiy2021, Bochkarev2024, &
VazdaCruz2021
CONTAINS
@ -1952,6 +1953,12 @@ CONTAINS
title="Graph Atomic Cluster Expansion for semilocal interactions beyond equivariant message passing", &
source="Phys. Rev. X", volume="14", pages="021036", &
year=2024, doi="10.1103/PhysRevX.14.021036")
CALL add_reference(key=VazdaCruz2021, &
authors=s2a("V. Vaz da Cruz", "S. Eckert", "A. Fohlisch"), &
title="TD-DFT simulations of K-edge resonant inelastic X-ray scattering within "// &
"the restricted subspace approximation", &
source="Phys. Chem. Chem. Phys.", volume="23", pages="1835-1848", &
year=2021, doi="10.1039/d0cp04726k")
END SUBROUTINE add_all_references

55
src/cp_control_types.F
View file

@ -29,6 +29,9 @@ MODULE cp_control_types
smeagol_control_type
USE xas_control, ONLY: xas_control_release,&
xas_control_type
USE xas_tdp_types, ONLY: xas_tdp_control_create,&
xas_tdp_control_release,&
xas_tdp_control_type
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -578,6 +581,16 @@ MODULE cp_control_types
TYPE(smeared_type), DIMENSION(:), POINTER :: smeared_occup => NULL()
END TYPE tddfpt2_control_type
! **************************************************************************************************
!> \brief
! **************************************************************************************************
TYPE rixs_control_type
TYPE(tddfpt2_control_type), POINTER :: tddfpt2_control => NULL()
TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control => NULL()
END TYPE rixs_control_type
! **************************************************************************************************
! \brief Control parameters for a DFT calculation
! \par History
@ -590,6 +603,7 @@ MODULE cp_control_types
TYPE(rtp_control_type), POINTER :: rtp_control => NULL()
TYPE(sccs_control_type), POINTER :: sccs_control => NULL()
TYPE(tddfpt2_control_type), POINTER :: tddfpt2_control => NULL()
TYPE(rixs_control_type), POINTER :: rixs_control => NULL()
TYPE(xas_control_type), POINTER :: xas_control => NULL()
TYPE(expot_control_type), POINTER :: expot_control => NULL()
TYPE(maxwell_control_type), POINTER :: maxwell_control => NULL()
@ -659,6 +673,7 @@ MODULE cp_control_types
qs_control_type, &
gapw_control_type, &
tddfpt2_control_type, &
rixs_control_type, &
proj_mo_type, &
efield_type, &
mulliken_restraint_type, &
@ -684,7 +699,9 @@ MODULE cp_control_types
admm_control_release, &
maxwell_control_create, &
expot_control_create, &
ddapc_control_create
ddapc_control_create, &
rixs_control_create, &
rixs_control_release
CONTAINS
@ -798,6 +815,7 @@ CONTAINS
NULLIFY (dft_control%xas_control)
NULLIFY (dft_control%qs_control)
NULLIFY (dft_control%tddfpt2_control)
NULLIFY (dft_control%rixs_control)
NULLIFY (dft_control%efield_fields)
NULLIFY (dft_control%period_efield)
NULLIFY (dft_control%admm_control)
@ -813,6 +831,7 @@ CONTAINS
dft_control%hairy_probes = .FALSE.
CALL qs_control_create(dft_control%qs_control)
CALL tddfpt2_control_create(dft_control%tddfpt2_control)
CALL rixs_control_create(dft_control%rixs_control)
CALL smeagol_control_create(dft_control%smeagol_control)
END SUBROUTINE dft_control_create
@ -830,6 +849,7 @@ CONTAINS
CALL qs_control_release(dft_control%qs_control)
CALL tddfpt2_control_release(dft_control%tddfpt2_control)
CALL rixs_control_release(dft_control%rixs_control) ! maybe check first if allocated
IF (ASSOCIATED(dft_control%xas_control)) THEN
CALL xas_control_release(dft_control%xas_control)
DEALLOCATE (dft_control%xas_control)
@ -971,6 +991,39 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE tddfpt2_control_release
! **************************************************************************************************
!> \brief Creates and initializes the rixs_control_type
!> \param rixs_control the type to initialize
! **************************************************************************************************
SUBROUTINE rixs_control_create(rixs_control)
TYPE(rixs_control_type), POINTER :: rixs_control
CPASSERT(.NOT. ASSOCIATED(rixs_control))
ALLOCATE (rixs_control)
NULLIFY (rixs_control%tddfpt2_control)
CALL tddfpt2_control_create(rixs_control%tddfpt2_control)
NULLIFY (rixs_control%xas_tdp_control)
CALL xas_tdp_control_create(rixs_control%xas_tdp_control)
END SUBROUTINE rixs_control_create
! **************************************************************************************************
!> \brief Releases the rixs_control_type
!> \param rixs_control ...
! **************************************************************************************************
SUBROUTINE rixs_control_release(rixs_control)
TYPE(rixs_control_type), POINTER :: rixs_control
IF (ASSOCIATED(rixs_control)) THEN
CALL tddfpt2_control_release(rixs_control%tddfpt2_control)
CALL xas_tdp_control_release(rixs_control%xas_tdp_control)
DEALLOCATE (rixs_control)
END IF
END SUBROUTINE rixs_control_release
! **************************************************************************************************
!> \brief ...
!> \param proj_mo_list ...

26
src/cp_control_utils.F
View file

@ -18,7 +18,8 @@ MODULE cp_control_utils
USE cp_control_types, ONLY: &
admm_control_create, admm_control_type, ddapc_control_create, ddapc_restraint_type, &
dft_control_create, dft_control_type, efield_type, expot_control_create, &
maxwell_control_create, qs_control_type, tddfpt2_control_type, xtb_control_type
maxwell_control_create, qs_control_type, rixs_control_type, tddfpt2_control_type, &
xtb_control_type
USE cp_files, ONLY: close_file,&
open_file
USE cp_log_handling, ONLY: cp_get_default_logger,&
@ -76,6 +77,7 @@ MODULE cp_control_utils
USE smeagol_control_types, ONLY: read_smeagol_control
USE string_utilities, ONLY: uppercase
USE util, ONLY: sort
USE xas_tdp_types, ONLY: read_xas_tdp_control
USE xc, ONLY: xc_uses_kinetic_energy_density,&
xc_uses_norm_drho
USE xc_input_constants, ONLY: xc_deriv_collocate
@ -89,6 +91,7 @@ MODULE cp_control_utils
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_control_utils'
PUBLIC :: read_dft_control, &
read_rixs_control, &
read_mgrid_section, &
read_qs_section, &
read_tddfpt2_control, &
@ -705,6 +708,27 @@ CONTAINS
END SUBROUTINE read_dft_control
! **************************************************************************************************
!> \brief Reads the input and stores in the rixs_control_type
!> \param rixs_control ...
!> \param rixs_section ...
!> \param qs_control ...
! **************************************************************************************************
SUBROUTINE read_rixs_control(rixs_control, rixs_section, qs_control)
TYPE(rixs_control_type), POINTER :: rixs_control
TYPE(section_vals_type), POINTER :: rixs_section
TYPE(qs_control_type), POINTER :: qs_control
TYPE(section_vals_type), POINTER :: td_section, xas_section
td_section => section_vals_get_subs_vals(rixs_section, "TDDFPT")
CALL read_tddfpt2_control(rixs_control%tddfpt2_control, td_section, qs_control)
xas_section => section_vals_get_subs_vals(rixs_section, "XAS_TDP")
CALL read_xas_tdp_control(rixs_control%xas_tdp_control, xas_section)
END SUBROUTINE read_rixs_control
! **************************************************************************************************
!> \brief ...
!> \param qs_control ...

21
src/header.F
View file

@ -23,7 +23,7 @@ MODULE header
PRIVATE
PUBLIC :: cp2k_header, cp2k_footer, vib_header
PUBLIC :: fist_header, se_header, dftb_header, qs_header, tddfpt_header
PUBLIC :: fist_header, se_header, dftb_header, qs_header, tddfpt_header, rixs_header
PUBLIC :: qmmm_header, atom_header, atom_footer, band_header, xtb_header, tblite_header
PUBLIC :: tmc_header, tmc_ana_header
PUBLIC :: sirius_header
@ -508,5 +508,24 @@ CONTAINS
"*******************************************************************************", &
"*******************************************************************************"
END SUBROUTINE tmc_ana_header
! **************************************************************************************************
!> \brief ...
!> \param iw ...
! **************************************************************************************************
SUBROUTINE rixs_header(iw)
INTEGER, INTENT(IN) :: iw
IF (iw < 0) RETURN
WRITE (iw, "(/,(T2,A79))") &
"*******************************************************************************", &
"** ######## #### ## ## ###### **", &
"** ## ## ## ## ## ## **", &
"** ######## ## ### ###### **", &
"** ## ## ## ## ## ## **", &
"** ## ## #### ## ## ###### **", &
"** Resonant Inelastic X-RAY Scattering calculations **", &
"*******************************************************************************"
END SUBROUTINE rixs_header
END MODULE header

51
src/input_cp2k_properties_dft.F
View file

@ -19,7 +19,8 @@ MODULE input_cp2k_properties_dft
Putrino2000, &
Putrino2002, &
Sebastiani2001, &
Weber2009
Weber2009, &
VazdaCruz2021
USE cp_output_handling, ONLY: add_last_numeric, &
cp_print_key_section_create, &
debug_print_level, &
@ -61,6 +62,7 @@ MODULE input_cp2k_properties_dft
lchar_t, &
logical_t, &
real_t
USE input_cp2k_xas, ONLY: create_xas_tdp_section
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
#include "./base/base_uses.f90"
@ -129,6 +131,10 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_rixs_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_bandstructure_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
@ -139,6 +145,49 @@ CONTAINS
END SUBROUTINE create_properties_section
! **************************************************************************************************
!> \brief creates the input structure used to activate
!> a resonant inelastic xray scattering (RIXS) calculation
! **************************************************************************************************
SUBROUTINE create_rixs_section(section)
TYPE(section_type), POINTER :: section
TYPE(section_type), POINTER :: subsection, print_key
TYPE(keyword_type), POINTER :: keyword
CPASSERT(.NOT. ASSOCIATED(section))
NULLIFY (keyword, subsection, print_key)
CALL section_create(section, __LOCATION__, name="RIXS", &
description="Resonant Inelastic Xray Scattering using XAS_TDP and TDDFPT.", &
n_keywords=1, n_subsections=3, repeats=.FALSE., &
citations=(/VazdaCruz2021/))
CALL create_tddfpt2_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_xas_tdp_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL section_create(subsection, __LOCATION__, "PRINT", "Controls the printing of information "// &
"during RIXS calculations", repeats=.FALSE.)
CALL cp_print_key_section_create(print_key, __LOCATION__, name="SPECTRUM", &
description="Controles the printing of the RIXS spectrum "// &
"in output files", &
print_level=low_print_level, filename="", &
common_iter_levels=3)
CALL section_add_subsection(subsection, print_key)
CALL section_release(print_key)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END SUBROUTINE create_rixs_section
! **************************************************************************************************
!> \brief creates the input structure used to activate
!> a linear response calculation

5
src/qs_energy_utils.F
View file

@ -62,6 +62,7 @@ MODULE qs_energy_utils
USE qs_tddfpt2_methods, ONLY: tddfpt
USE qs_vxc, ONLY: qs_xc_density
USE qs_vxc_atom, ONLY: calculate_vxc_atom
USE rixs_methods, ONLY: rixs
USE tip_scan_methods, ONLY: tip_scanning
USE xas_methods, ONLY: xas
USE xas_tdp_methods, ONLY: xas_tdp
@ -180,6 +181,7 @@ CONTAINS
END IF
!Properties
IF (dft_control%do_xas_calculation) THEN
CALL xas(qs_env, dft_control)
END IF
@ -197,6 +199,9 @@ CONTAINS
CALL tddfpt(qs_env, calc_forces)
END IF
! stand-alone RIXS, does not depend on previous xas_tdp and/or tddfpt2 calculations
IF (qs_env%do_rixs) CALL rixs(qs_env)
! post-SCF bandstructure calculation from higher level methods
NULLIFY (post_scf_bands_section)
post_scf_bands_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%BANDSTRUCTURE")

26
src/qs_environment.F
View file

@ -41,14 +41,10 @@ MODULE qs_environment
qs_control_type,&
semi_empirical_control_type,&
xtb_control_type
USE cp_control_utils, ONLY: read_ddapc_section,&
read_dft_control,&
read_mgrid_section,&
read_qs_section,&
read_tddfpt2_control,&
write_admm_control,&
write_dft_control,&
write_qs_control
USE cp_control_utils, ONLY: &
read_ddapc_section, read_dft_control, read_mgrid_section, read_qs_section, &
read_rixs_control, read_tddfpt2_control, write_admm_control, write_dft_control, &
write_qs_control
USE cp_ddapc_types, ONLY: cp_ddapc_ewald_create
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
@ -652,7 +648,8 @@ CONTAINS
TYPE(se_taper_type), POINTER :: se_taper
TYPE(section_vals_type), POINTER :: dft_section, et_coupling_section, et_ddapc_section, &
ewald_section, harris_section, lri_section, mp2_section, nl_section, poisson_section, &
pp_section, print_section, qs_section, se_section, tddfpt_section, xc_section
pp_section, print_section, qs_section, rixs_section, se_section, tddfpt_section, &
xc_section
TYPE(semi_empirical_control_type), POINTER :: se_control
TYPE(semi_empirical_si_type), POINTER :: se_store_int_env
TYPE(xtb_control_type), POINTER :: xtb_control
@ -684,6 +681,7 @@ CONTAINS
et_coupling_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%ET_COUPLING")
! reimplemented TDDFPT
tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
rixs_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS")
CALL qs_subsys_get(subsys, particle_set=particle_set, &
qs_kind_set=qs_kind_set, &
@ -751,6 +749,12 @@ CONTAINS
! reimplemented TDDFPT
CALL read_tddfpt2_control(dft_control%tddfpt2_control, tddfpt_section, dft_control%qs_control)
! rixs
CALL section_vals_get(rixs_section, explicit=qs_env%do_rixs)
IF (qs_env%do_rixs) THEN
CALL read_rixs_control(dft_control%rixs_control, rixs_section, dft_control%qs_control)
END IF
! Create relativistic control section
BLOCK
TYPE(rel_control_type), POINTER :: rel_control
@ -1058,7 +1062,7 @@ CONTAINS
END IF
IF (dft_control%do_xas_tdp_calculation) THEN
IF (dft_control%do_xas_tdp_calculation .OR. qs_env%do_rixs) THEN
! Check if RI_XAS basis is given, auto-generate if not
CALL get_qs_env(qs_env, nkind=nkind)
DO ikind = 1, nkind
@ -1090,7 +1094,7 @@ CONTAINS
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_HXC")
maxlgto = MAX(maxlgto, maxlgto_lri)
END IF
IF (dft_control%do_xas_tdp_calculation) THEN
IF (dft_control%do_xas_tdp_calculation .OR. qs_env%do_rixs) THEN
!done as a precaution
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_XAS")
maxlgto = MAX(maxlgto, maxlgto_lri)

53
src/qs_environment_types.F
View file

@ -229,6 +229,7 @@ MODULE qs_environment_types
LOGICAL :: given_embed_pot = .FALSE.
LOGICAL :: energy_correction = .FALSE.
LOGICAL :: harris_method = .FALSE.
LOGICAL :: do_rixs = .FALSE.
REAL(KIND=dp) :: sim_time = -1.0_dp
REAL(KIND=dp) :: start_time = -1.0_dp, target_time = -1.0_dp
REAL(KIND=dp), DIMENSION(:, :), POINTER :: image_matrix => NULL()
@ -483,6 +484,7 @@ CONTAINS
!> \param mos_last_converged ... [SGh]
!> \param eeq ...
!> \param rhs ...
!> \param do_rixs ...
!> \param tb_tblite ...
!> \date 23.01.2002
!> \author MK
@ -514,7 +516,8 @@ CONTAINS
WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, &
s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, &
gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, &
eeq, rhs, tb_tblite)
eeq, rhs, do_rixs, tb_tblite)
TYPE(qs_environment_type), INTENT(IN) :: qs_env
TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, &
POINTER :: atomic_kind_set
@ -648,6 +651,7 @@ CONTAINS
TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env
TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos_last_converged
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eeq, rhs
LOGICAL, OPTIONAL :: do_rixs
TYPE(tblite_type), OPTIONAL, POINTER :: tb_tblite
TYPE(rho0_mpole_type), POINTER :: rho0_m
@ -707,6 +711,7 @@ CONTAINS
IF (PRESENT(mscfg_env)) mscfg_env => qs_env%molecular_scf_guess_env
IF (PRESENT(active_space)) active_space => qs_env%active_space
IF (PRESENT(admm_env)) admm_env => qs_env%admm_env
IF (PRESENT(do_rixs)) do_rixs = qs_env%do_rixs
! Embedding potential
IF (PRESENT(embed_pot)) embed_pot => qs_env%embed_pot
@ -939,6 +944,7 @@ CONTAINS
qs_env%calc_image_preconditioner = .TRUE.
qs_env%do_transport = .FALSE.
qs_env%given_embed_pot = .FALSE.
qs_env%do_rixs = .FALSE.
IF (PRESENT(globenv)) THEN
qs_env%target_time = globenv%cp2k_target_time
qs_env%start_time = globenv%cp2k_start_time
@ -1041,6 +1047,7 @@ CONTAINS
!> \param mos_last_converged ... [SGh]
!> \param eeq ...
!> \param rhs ...
!> \param do_rixs ...
!> \param tb_tblite ...
!> \date 23.01.2002
!> \author MK
@ -1061,7 +1068,7 @@ CONTAINS
do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, &
harris_env, gcp_env, mp2_env, bs_env, kg_env, force, &
kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, &
spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, tb_tblite)
spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
TYPE(qs_environment_type), INTENT(INOUT) :: qs_env
TYPE(cell_type), OPTIONAL, POINTER :: super_cell
@ -1138,6 +1145,7 @@ CONTAINS
TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env
TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos_last_converged
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eeq, rhs
LOGICAL, OPTIONAL :: do_rixs
TYPE(tblite_type), OPTIONAL, POINTER :: tb_tblite
TYPE(qs_subsys_type), POINTER :: subsys
@ -1174,13 +1182,14 @@ CONTAINS
IF (PRESENT(rtp)) qs_env%rtp => rtp
IF (PRESENT(efield)) qs_env%efield => efield
IF (PRESENT(active_space)) qs_env%active_space => active_space
IF (PRESENT(do_rixs)) do_rixs = qs_env%do_rixs
IF (PRESENT(ewald_env)) THEN ! accept also null pointers?
IF (ASSOCIATED(qs_env%ewald_env)) THEN
IF (.NOT. ASSOCIATED(qs_env%ewald_env, ewald_env)) THEN
CALL ewald_env_release(qs_env%ewald_env)
DEALLOCATE (qs_env%ewald_env)
END IF
IF (.NOT. ASSOCIATED(qs_env%ewald_env, ewald_env)) THEN
CALL ewald_env_release(qs_env%ewald_env)
DEALLOCATE (qs_env%ewald_env)
END IF
END IF
qs_env%ewald_env => ewald_env
END IF
@ -1213,20 +1222,20 @@ CONTAINS
END IF
IF (PRESENT(linres_control)) THEN ! accept also null pointers?
IF (ASSOCIATED(qs_env%linres_control)) THEN
IF (.NOT. ASSOCIATED(qs_env%linres_control, linres_control)) THEN
CALL linres_control_release(qs_env%linres_control)
DEALLOCATE (qs_env%linres_control)
END IF
IF (.NOT. ASSOCIATED(qs_env%linres_control, linres_control)) THEN
CALL linres_control_release(qs_env%linres_control)
DEALLOCATE (qs_env%linres_control)
END IF
END IF
qs_env%linres_control => linres_control
END IF
! ZMP associating variables
IF (PRESENT(rho_external)) THEN
IF (ASSOCIATED(qs_env%rho_external)) THEN
IF (.NOT. ASSOCIATED(qs_env%rho_external, rho_external)) THEN
CALL qs_rho_release(qs_env%rho_external)
DEALLOCATE (qs_env%rho_external)
END IF
IF (.NOT. ASSOCIATED(qs_env%rho_external, rho_external)) THEN
CALL qs_rho_release(qs_env%rho_external)
DEALLOCATE (qs_env%rho_external)
END IF
END IF
qs_env%rho_external => rho_external
END IF
@ -1273,10 +1282,10 @@ CONTAINS
END IF
IF (PRESENT(scf_env)) THEN ! accept also null pointers ?
IF (ASSOCIATED(qs_env%scf_env)) THEN
IF (.NOT. ASSOCIATED(qs_env%scf_env, scf_env)) THEN
CALL scf_env_release(qs_env%scf_env)
DEALLOCATE (qs_env%scf_env)
END IF
IF (.NOT. ASSOCIATED(qs_env%scf_env, scf_env)) THEN
CALL scf_env_release(qs_env%scf_env)
DEALLOCATE (qs_env%scf_env)
END IF
END IF
qs_env%scf_env => scf_env
END IF
@ -1344,10 +1353,10 @@ CONTAINS
END IF
IF (PRESENT(se_nonbond_env)) THEN
IF (ASSOCIATED(qs_env%se_nonbond_env)) THEN
IF (.NOT. ASSOCIATED(qs_env%se_nonbond_env, se_nonbond_env)) THEN
CALL fist_nonbond_env_release(qs_env%se_nonbond_env)
DEALLOCATE (qs_env%se_nonbond_env)
END IF
IF (.NOT. ASSOCIATED(qs_env%se_nonbond_env, se_nonbond_env)) THEN
CALL fist_nonbond_env_release(qs_env%se_nonbond_env)
DEALLOCATE (qs_env%se_nonbond_env)
END IF
END IF
qs_env%se_nonbond_env => se_nonbond_env
END IF

2
src/qs_scf_initialization.F
View file

@ -660,6 +660,8 @@ CONTAINS
END IF
IF (dft_control%do_xas_calculation) &
CPABORT("No XAS implemented with kpoints")
IF (qs_env%do_rixs) &
CPABORT("RIXS not implemented with kpoints")
DO ik = 1, SIZE(kpoints%kp_env)
CALL mpools_get(kpoints%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
mos_k => kpoints%kp_env(ik)%kpoint_env%mos

55
src/qs_tddfpt2_methods.F
View file

@ -17,6 +17,7 @@ MODULE qs_tddfpt2_methods
USE cell_types, ONLY: cell_type
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type,&
rixs_control_type,&
tddfpt2_control_type
USE cp_dbcsr_api, ONLY: dbcsr_p_type
USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
@ -108,6 +109,8 @@ MODULE qs_tddfpt2_methods
tddfpt_init_mos,&
tddfpt_oecorr,&
tddfpt_release_ground_state_mos
USE rixs_types, ONLY: rixs_env_type,&
tddfpt2_valence_type
USE string_utilities, ONLY: integer_to_string
USE xc_write_output, ONLY: xc_write
#include "./base/base_uses.f90"
@ -133,15 +136,17 @@ CONTAINS
!> \brief Perform TDDFPT calculation.
!> \param qs_env Quickstep environment
!> \param calc_forces ...
!> \param rixs_env ...
!> \par History
!> * 05.2016 created [Sergey Chulkov]
!> * 06.2016 refactored to be used with Davidson eigensolver [Sergey Chulkov]
!> * 03.2017 cleaned and refactored [Sergey Chulkov]
!> \note Based on the subroutines tddfpt_env_init(), and tddfpt_env_deallocate().
! **************************************************************************************************
SUBROUTINE tddfpt(qs_env, calc_forces)
SUBROUTINE tddfpt(qs_env, calc_forces, rixs_env)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN) :: calc_forces
TYPE(rixs_env_type), OPTIONAL, POINTER :: rixs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt'
@ -150,7 +155,7 @@ CONTAINS
nstate_max, nstates, nvirt, old_state
INTEGER, DIMENSION(maxspins) :: nactive
LOGICAL :: do_admm, do_exck, do_hfx, do_hfxlr, &
do_hfxsr, do_soc, lmult_tmp, &
do_hfxsr, do_rixs, do_soc, lmult_tmp, &
state_change
REAL(kind=dp) :: gsmin, gsval, xsval
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength
@ -173,12 +178,14 @@ CONTAINS
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_scf_env_type), POINTER :: scf_env
TYPE(rixs_control_type), POINTER :: rixs_control
TYPE(section_vals_type), POINTER :: hfxsr_section, kernel_section, &
lri_section, soc_section, &
tddfpt_print_section, tddfpt_section, &
xc_section
TYPE(stda_env_type), TARGET :: stda_kernel
TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
TYPE(tddfpt2_valence_type), POINTER :: valence_state
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
POINTER :: gs_mos
TYPE(tddfpt_subgroup_env_type) :: sub_env
@ -190,8 +197,7 @@ CONTAINS
logger => cp_get_default_logger()
! input section print/xc
NULLIFY (tddfpt_section)
tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
NULLIFY (tddfpt_section, tddfpt_control)
CALL tddfpt_input(qs_env, do_hfx, do_admm, do_exck, &
do_hfxsr, do_hfxlr, xc_section, tddfpt_print_section, &
@ -207,8 +213,20 @@ CONTAINS
matrix_ks=matrix_ks, &
matrix_s=matrix_s, &
mos=mos, &
scf_env=scf_env)
tddfpt_control => dft_control%tddfpt2_control
scf_env=scf_env, &
do_rixs=do_rixs)
IF (do_rixs) THEN
tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS%TDDFPT")
NULLIFY (rixs_control, valence_state)
rixs_control => dft_control%rixs_control
tddfpt_control => rixs_control%tddfpt2_control
valence_state => rixs_env%valence_state
ELSE
tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
tddfpt_control => dft_control%tddfpt2_control
END IF
tddfpt_control%do_hfx = do_hfx
tddfpt_control%do_admm = do_admm
tddfpt_control%do_hfxsr = do_hfxsr
@ -495,6 +513,31 @@ CONTAINS
END IF
END IF
! share evals, evects and mo_coefs with rixs
IF (do_rixs) THEN
! copy evals
valence_state%nstates = nstates
ALLOCATE (valence_state%evals(SIZE(evals)))
valence_state%evals(:) = evals(:)
ALLOCATE (valence_state%evects(nspins, nstates))
ALLOCATE (valence_state%mos_occ(nspins))
DO ispin = 1, nspins
! copy evects
DO istate = 1, nstates
CALL cp_fm_get_info(matrix=evects(ispin, istate), &
matrix_struct=matrix_struct)
CALL cp_fm_create(valence_state%evects(ispin, istate), matrix_struct)
CALL cp_fm_to_fm(evects(ispin, istate), valence_state%evects(ispin, istate))
END DO
! copy mos_occ
CALL cp_fm_get_info(matrix=gs_mos(ispin)%mos_occ, &
matrix_struct=matrix_struct)
CALL cp_fm_create(valence_state%mos_occ(ispin), matrix_struct)
CALL cp_fm_to_fm(gs_mos(ispin)%mos_occ, valence_state%mos_occ(ispin))
END DO
END IF
! clean up
CALL cp_fm_release(evects)
CALL cp_fm_release(S_evects)

470
src/rixs_methods.F Normal file
View file

@ -0,0 +1,470 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Methods for Resonant Inelastic XRAY Scattering (RIXS) calculations
!> \author BSG (02.2025)
! **************************************************************************************************
MODULE rixs_methods
USE bibliography, ONLY: VazdaCruz2021,&
cite_reference
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type,&
rixs_control_create,&
rixs_control_release,&
rixs_control_type
USE cp_control_utils, ONLY: read_rixs_control
USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
dbcsr_type
USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_get_info,&
cp_fm_get_submatrix,&
cp_fm_release,&
cp_fm_to_fm,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE header, ONLY: rixs_header
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_para_env_type
USE parallel_gemm_api, ONLY: parallel_gemm
USE physcon, ONLY: evolt
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_tddfpt2_methods, ONLY: tddfpt
USE rixs_types, ONLY: rixs_env_create,&
rixs_env_release,&
rixs_env_type,&
tddfpt2_valence_type
USE xas_tdp_methods, ONLY: xas_tdp
USE xas_tdp_types, ONLY: donor_state_type,&
xas_tdp_env_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rixs_methods'
PUBLIC :: rixs, rixs_core
CONTAINS
! **************************************************************************************************
!> \brief Driver for RIXS calculations.
!> \param qs_env the inherited qs_environment
!> \author BSG
! **************************************************************************************************
SUBROUTINE rixs(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'rixs'
INTEGER :: handle, output_unit
TYPE(dft_control_type), POINTER :: dft_control
TYPE(section_vals_type), POINTER :: rixs_section, tddfp2_section, &
xas_tdp_section
CALL timeset(routineN, handle)
NULLIFY (rixs_section)
rixs_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS")
output_unit = cp_logger_get_default_io_unit()
qs_env%do_rixs = .TRUE.
CALL cite_reference(VazdaCruz2021)
CALL get_qs_env(qs_env, dft_control=dft_control)
IF (dft_control%uks .OR. dft_control%roks) CPABORT("RIXS not implemented for LSD/ROKS")
xas_tdp_section => section_vals_get_subs_vals(rixs_section, "XAS_TDP")
tddfp2_section => section_vals_get_subs_vals(rixs_section, "TDDFPT")
IF (.NOT. ASSOCIATED(xas_tdp_section)) THEN
CPABORT("XAS_TDP calculation missing")
END IF
IF (.NOT. ASSOCIATED(tddfp2_section)) THEN
CPABORT("TDDFPT calculation missing")
END IF
CALL rixs_core(rixs_section, qs_env)
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,(T2,A79))") &
"*******************************************************************************", &
"! Normal termination of Resonant Inelastic X-RAY Scattering calculation !", &
"*******************************************************************************"
END IF
CALL timestop(handle)
END SUBROUTINE rixs
! **************************************************************************************************
!> \brief Perform RIXS calculation.
!> \param rixs_section ...
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE rixs_core(rixs_section, qs_env)
TYPE(section_vals_type), POINTER :: rixs_section
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'rixs_core'
INTEGER :: ax, current_state_index, fstate, handle, &
iatom, istate, nao, nex_atoms, nocc, &
nstates, nvirt, output_unit, td_state
REAL(dp) :: osc_xyz
REAL(dp), ALLOCATABLE, DIMENSION(:) :: w_i0, w_if
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dip_block, mu_i0
REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: mu_if
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_struct_type), POINTER :: dip_0_struct, dip_f_struct, &
i_dip_0_struct, i_dip_f_struct
TYPE(cp_fm_type) :: dip_0, dip_f, i_dip_0, i_dip_f
TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: valence_evects
TYPE(cp_fm_type), POINTER :: core_evects, local_gs_coeffs, mo_coeffs
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(donor_state_type), POINTER :: current_state
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(rixs_control_type), POINTER :: rixs_control
TYPE(rixs_env_type), POINTER :: rixs_env
TYPE(tddfpt2_valence_type), POINTER :: valence_state
TYPE(xas_tdp_env_type), POINTER :: core_state
NULLIFY (rixs_control, dft_control, rixs_env)
NULLIFY (valence_state, core_state)
NULLIFY (para_env, blacs_env)
NULLIFY (local_gs_coeffs, mo_coeffs, valence_evects)
NULLIFY (dipmat, dip_0_struct, i_dip_0_struct, dip_f_struct, i_dip_f_struct)
output_unit = cp_logger_get_default_io_unit()
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
matrix_s=matrix_s, &
para_env=para_env, &
blacs_env=blacs_env)
CALL rixs_control_create(rixs_control)
CALL read_rixs_control(rixs_control, rixs_section, dft_control%qs_control)
! create rixs_env
CALL rixs_env_create(rixs_env)
! first, xas_tdp calculation
CALL xas_tdp(qs_env, rixs_env)
IF (rixs_control%xas_tdp_control%check_only) THEN
CPWARN("CHECK_ONLY run for XAS_TDP requested, RIXS will not be performed.")
ELSE
! then, tddfpt2 calculation
CALL tddfpt(qs_env, calc_forces=.FALSE., rixs_env=rixs_env)
IF (output_unit > 0) THEN
CALL rixs_header(output_unit)
END IF
! timings for rixs only, excluding xas_tdp and tddft calls
CALL timeset(routineN, handle)
core_state => rixs_env%core_state
valence_state => rixs_env%valence_state
! gs coefficients from tddfpt
mo_coeffs => valence_state%mos_occ(1)
! localised gs coefficients from xas_tdp
local_gs_coeffs => core_state%mo_coeff(1) ! TODO (1)=ispin
valence_evects => valence_state%evects
IF (rixs_control%xas_tdp_control%do_loc) THEN
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(T2,A)") &
"RIXS| Found localised XAS_TDP orbitals"
WRITE (UNIT=output_unit, FMT="(T2,A)") &
"RIXS| Rotating TDDFPT vectors..."
END IF
CALL rotate_vectors(valence_evects, local_gs_coeffs, mo_coeffs, matrix_s(1)%matrix, output_unit)
END IF
CALL cp_fm_get_info(matrix=valence_evects(1, 1), nrow_global=nao, ncol_global=nocc) ! TODO evects
nex_atoms = core_state%nex_atoms
nstates = valence_state%nstates
dipmat => core_state%dipmat
nvirt = core_state%nvirt
ALLOCATE (dip_block(1, 1))
ALLOCATE (mu_i0(4, nvirt))
mu_i0 = 0.0_dp
ALLOCATE (mu_if(4, nvirt, nstates)) ! mu per (donor state -> nstate) per (x,y,z) ! experimental
mu_if = 0.0_dp
ALLOCATE (w_i0(nvirt), w_if(nstates))
w_if(:) = valence_state%evals(:)*evolt
! initialise matrices for i->0
CALL cp_fm_struct_create(dip_0_struct, para_env=para_env, context=blacs_env, &
nrow_global=nao, ncol_global=1)
CALL cp_fm_create(dip_0, dip_0_struct)
CALL cp_fm_struct_create(i_dip_0_struct, para_env=para_env, context=blacs_env, &
nrow_global=nvirt, ncol_global=1)
CALL cp_fm_create(i_dip_0, i_dip_0_struct)
! initialise matrices for i->f
CALL cp_fm_struct_create(dip_f_struct, para_env=para_env, context=blacs_env, &
nrow_global=nao, ncol_global=nocc)
CALL cp_fm_create(dip_f, dip_f_struct)
CALL cp_fm_struct_create(i_dip_f_struct, para_env=para_env, context=blacs_env, &
nrow_global=nvirt, ncol_global=nocc)
CALL cp_fm_create(i_dip_f, i_dip_f_struct)
! looping over ex_atoms and ex_kinds is enough as excited atoms have to be unique
current_state_index = 1
DO iatom = 1, nex_atoms
current_state => core_state%donor_states(current_state_index)
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(T2,A,A,A,A)") &
"RIXS| Calculating dipole moment from core-excited state ", &
core_state%state_type_char(current_state%state_type), " of ", TRIM(current_state%at_symbol)
END IF
core_evects => current_state%sg_coeffs
w_i0(:) = current_state%sg_evals(:)*evolt
! 0 -> i
DO ax = 1, 3
! R*0
CALL cp_dbcsr_sm_fm_multiply(dipmat(ax)%matrix, current_state%gs_coeffs, dip_0, ncol=1)
! i*R*0
CALL parallel_gemm('T', 'N', nvirt, 1, nao, 1.0_dp, core_evects, dip_0, 0.0_dp, i_dip_0)
DO istate = 1, nvirt
CALL cp_fm_get_submatrix(fm=i_dip_0, target_m=dip_block, start_row=istate, &
start_col=1, n_rows=1, n_cols=1)
mu_i0(ax, istate) = dip_block(1, 1)
osc_xyz = mu_i0(ax, istate)**2
mu_i0(4, istate) = mu_i0(4, istate) + osc_xyz
END DO ! istate
END DO ! ax
! i -> f
DO td_state = 1, nstates
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(T9,A,I3,A,F10.4)") &
"to valence-excited state ", td_state, " with energy ", w_if(td_state)
END IF
DO ax = 1, 3
! core_evects x dipmat x valence_evects
CALL cp_dbcsr_sm_fm_multiply(dipmat(ax)%matrix, valence_evects(1, td_state), dip_f, ncol=nocc)
CALL parallel_gemm('T', 'N', nvirt, nocc, nao, 1.0_dp, core_evects, dip_f, 0.0_dp, i_dip_f)
DO istate = 1, nvirt
DO fstate = 1, nocc ! 5
CALL cp_fm_get_submatrix(fm=i_dip_f, target_m=dip_block, start_row=istate, &
start_col=fstate, n_rows=1, n_cols=1)
mu_if(ax, istate, td_state) = mu_if(ax, istate, td_state) + dip_block(1, 1)
END DO ! fstate (tddft)
osc_xyz = mu_if(ax, istate, td_state)**2
mu_if(4, istate, td_state) = mu_if(4, istate, td_state) + osc_xyz
END DO ! istate (core)
END DO ! ax
END DO ! td_state
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T2,A,/)") "RIXS| Printing spectrum to file"
END IF
CALL print_rixs_to_file(current_state, mu_i0, mu_if, w_i0, w_if, rixs_env, rixs_section)
current_state_index = current_state_index + 1
END DO ! iatom
NULLIFY (current_state)
! cleanup
CALL cp_fm_struct_release(i_dip_0_struct)
CALL cp_fm_struct_release(dip_0_struct)
CALL cp_fm_release(dip_0)
CALL cp_fm_release(i_dip_0)
CALL cp_fm_struct_release(i_dip_f_struct)
CALL cp_fm_struct_release(dip_f_struct)
CALL cp_fm_release(dip_f)
CALL cp_fm_release(i_dip_f)
END IF
! nullify rixs_control, rixs_env
CALL rixs_control_release(rixs_control)
CALL rixs_env_release(rixs_env)
NULLIFY (valence_state, core_state)
CALL timestop(handle)
END SUBROUTINE rixs_core
! **************************************************************************************************
!> \brief Rotate vectors. Returns rotated mo_occ and evects.
!> \param evects ...
!> \param mo_ref ...
!> \param mo_occ ...
!> \param overlap_matrix ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE rotate_vectors(evects, mo_ref, mo_occ, overlap_matrix, unit_nr)
TYPE(cp_fm_type), DIMENSION(:, :) :: evects
TYPE(cp_fm_type) :: mo_ref, mo_occ
TYPE(dbcsr_type), POINTER :: overlap_matrix
INTEGER :: unit_nr
INTEGER :: istate, ncol, nrow, nstates
REAL(kind=dp) :: diff
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_struct_type), POINTER :: emat_struct
TYPE(cp_fm_type) :: emat, rotated_mo_coeffs, smo
TYPE(cp_fm_type), POINTER :: current_evect
TYPE(mp_para_env_type), POINTER :: para_env
NULLIFY (emat_struct, para_env, blacs_env, current_evect)
CALL cp_fm_get_info(matrix=mo_occ, nrow_global=nrow, ncol_global=ncol, &
para_env=para_env, context=blacs_env)
CALL cp_fm_create(smo, mo_occ%matrix_struct)
! rotate mo_occ
! smo = matrix_s x mo_occ
CALL cp_dbcsr_sm_fm_multiply(overlap_matrix, mo_occ, smo, ncol, alpha=1.0_dp, beta=0.0_dp)
CALL cp_fm_struct_create(emat_struct, nrow_global=ncol, ncol_global=ncol, &
para_env=para_env, context=blacs_env)
CALL cp_fm_create(emat, emat_struct)
! emat = mo_ref^T x smo
CALL parallel_gemm('T', 'N', ncol, ncol, nrow, 1.0_dp, mo_ref, smo, 0.0_dp, emat)
CALL cp_fm_create(rotated_mo_coeffs, mo_occ%matrix_struct)
! rotated_mo_coeffs = cpmos x emat
CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, mo_occ, emat, 0.0_dp, rotated_mo_coeffs)
diff = MAXVAL(ABS(rotated_mo_coeffs%local_data - mo_occ%local_data))
IF (unit_nr > 0) THEN
WRITE (unit_nr, FMT="(T9,A,F10.6,/)") "Max difference between orbitals = ", diff
END IF
CALL cp_fm_to_fm(rotated_mo_coeffs, mo_occ)
nstates = SIZE(evects, DIM=2)
DO istate = 1, nstates
ASSOCIATE (current_evect => evects(1, istate))
CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, current_evect, emat, 0.0_dp, smo)
diff = MAXVAL(ABS(smo%local_data - current_evect%local_data))
CALL cp_fm_to_fm(smo, current_evect)
END ASSOCIATE
END DO
CALL cp_fm_struct_release(emat_struct)
CALL cp_fm_release(smo)
CALL cp_fm_release(emat)
CALL cp_fm_release(rotated_mo_coeffs)
END SUBROUTINE rotate_vectors
!**************************************************************************************************
!> \brief Print RIXS spectrum.
!> \param donor_state ...
!> \param mu_i0 ...
!> \param mu_if ...
!> \param w_i0 ...
!> \param w_if ...
!> \param rixs_env ...
!> \param rixs_section ...
! **************************************************************************************************
SUBROUTINE print_rixs_to_file(donor_state, mu_i0, mu_if, w_i0, w_if, &
rixs_env, rixs_section)
TYPE(donor_state_type), POINTER :: donor_state
REAL(dp), DIMENSION(:, :) :: mu_i0
REAL(dp), DIMENSION(:, :, :) :: mu_if
REAL(dp), DIMENSION(:) :: w_i0, w_if
TYPE(rixs_env_type), POINTER :: rixs_env
TYPE(section_vals_type), POINTER :: rixs_section
INTEGER :: f, i, output_unit, rixs_unit
TYPE(cp_logger_type), POINTER :: logger
NULLIFY (logger)
logger => cp_get_default_logger()
rixs_unit = cp_print_key_unit_nr(logger, rixs_section, "PRINT%SPECTRUM", &
extension=".rixs", file_position="APPEND", &
file_action="WRITE", file_form="FORMATTED")
output_unit = cp_logger_get_default_io_unit()
IF (rixs_unit > 0) THEN
WRITE (rixs_unit, FMT="(A,/,T2,A,A,A,A,A,/,A)") &
"====================================================================================", &
"Excitation from ground-state (", &
rixs_env%core_state%state_type_char(donor_state%state_type), " of kind ", &
TRIM(donor_state%at_symbol), ") to core-excited state i ", &
"===================================================================================="
WRITE (rixs_unit, FMT="(T3,A)") &
"w_0i (eV) mu^x_0i (a.u.) mu^y_0i (a.u.) mu^z_0i (a.u.) mu^2_0i (a.u.)"
DO i = 1, SIZE(mu_i0, DIM=2)
WRITE (rixs_unit, FMT="(T2,F10.4,T26,E12.5,T42,E12.5,T58,E12.5,T74,E12.5)") &
w_i0(i), mu_i0(1, i), mu_i0(2, i), mu_i0(3, i), mu_i0(4, i)
END DO
WRITE (rixs_unit, FMT="(A,/,T2,A,/,A)") &
"====================================================================================", &
"Emission from core-excited state i to valence-excited state f ", &
"===================================================================================="
WRITE (rixs_unit, FMT="(T3,A)") &
"w_0i (eV) w_if (eV) mu^x_if (a.u.) mu^y_if (a.u.) mu^z_if (a.u.) mu^2_if (a.u.)"
DO i = 1, SIZE(mu_if, DIM=2)
DO f = 1, SIZE(mu_if, DIM=3)
WRITE (rixs_unit, FMT="(T2,F10.4,T14,F8.4,T26,E12.5,T42,E12.5,T58,E12.5,T74,E12.5)") &
w_i0(i), w_if(f), mu_if(1, i, f), mu_if(2, i, f), mu_if(3, i, f), mu_if(4, i, f)
END DO
END DO
END IF
CALL cp_print_key_finished_output(rixs_unit, logger, rixs_section, "PRINT%SPECTRUM")
END SUBROUTINE print_rixs_to_file
END MODULE rixs_methods

135
src/rixs_types.F Normal file
View file

@ -0,0 +1,135 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
!> *************************************************************************************************
!> \brief Define Resonant Inelastic XRAY Scattering (RIXS) control type and associated create,
!> release, etc subroutines
!> \author BSG (02.2025)
!> *************************************************************************************************
MODULE rixs_types
USE cp_fm_types, ONLY: cp_fm_release,&
cp_fm_type
USE kinds, ONLY: dp
USE xas_tdp_types, ONLY: xas_tdp_env_create,&
xas_tdp_env_release,&
xas_tdp_env_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
! **************************************************************************************************
!> \brief
! **************************************************************************************************
TYPE rixs_env_type
TYPE(xas_tdp_env_type), POINTER :: core_state => NULL()
TYPE(tddfpt2_valence_type), &
POINTER :: valence_state => NULL()
END TYPE rixs_env_type
! **************************************************************************************************
!> \brief Valence state coming from the qs_tddfpt2 routine
! **************************************************************************************************
TYPE tddfpt2_valence_type
INTEGER :: nstates = 0
TYPE(cp_fm_type), DIMENSION(:, :), &
POINTER :: evects => NULL() ! eigenvectors
REAL(dp), DIMENSION(:), ALLOCATABLE :: evals ! energies
! entities below are coming from tddfpt_ground_state_mos type
TYPE(cp_fm_type), DIMENSION(:), &
POINTER :: mos_occ => NULL()
END TYPE tddfpt2_valence_type
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rixs_types'
PUBLIC :: rixs_env_type, tddfpt2_valence_type, tddfpt2_valence_state_create
PUBLIC :: rixs_env_create, rixs_env_release
CONTAINS
! **************************************************************************************************
!> \brief Creates a rixs environment type
!> \param rixs_env the type to create
! **************************************************************************************************
SUBROUTINE rixs_env_create(rixs_env)
TYPE(rixs_env_type), POINTER :: rixs_env
ALLOCATE (rixs_env)
NULLIFY (rixs_env%core_state)
CALL xas_tdp_env_create(rixs_env%core_state)
NULLIFY (rixs_env%valence_state)
CALL tddfpt2_valence_state_create(rixs_env%valence_state)
END SUBROUTINE rixs_env_create
! **************************************************************************************************
!> \brief Releases the rixs environment type
!> \param rixs_env the type to release
! **************************************************************************************************
SUBROUTINE rixs_env_release(rixs_env)
TYPE(rixs_env_type), POINTER :: rixs_env
IF (ASSOCIATED(rixs_env)) THEN
IF (ASSOCIATED(rixs_env%core_state)) THEN
CALL xas_tdp_env_release(rixs_env%core_state)
END IF
IF (ASSOCIATED(rixs_env%valence_state)) THEN
CALL tddfpt2_valence_state_release(rixs_env%valence_state)
END IF
END IF
DEALLOCATE (rixs_env)
END SUBROUTINE rixs_env_release
! **************************************************************************************************
!> \brief Creates the valence state type
!> \param valence_state ...
! **************************************************************************************************
SUBROUTINE tddfpt2_valence_state_create(valence_state)
TYPE(tddfpt2_valence_type), POINTER :: valence_state
ALLOCATE (valence_state)
NULLIFY (valence_state%evects)
! entities below come from tddfpt_ground_state_mos type
NULLIFY (valence_state%mos_occ)
END SUBROUTINE tddfpt2_valence_state_create
! **************************************************************************************************
!> \brief Releases the valence state type
!> \param valence_state ...
! **************************************************************************************************
SUBROUTINE tddfpt2_valence_state_release(valence_state)
TYPE(tddfpt2_valence_type), POINTER :: valence_state
IF (ASSOCIATED(valence_state)) THEN
IF (ASSOCIATED(valence_state%evects)) THEN
CALL cp_fm_release(valence_state%evects)
END IF
IF (ALLOCATED(valence_state%evals)) THEN
DEALLOCATE (valence_state%evals)
END IF
IF (ASSOCIATED(valence_state%mos_occ)) THEN
CALL cp_fm_release(valence_state%mos_occ)
END IF
END IF
DEALLOCATE (valence_state)
END SUBROUTINE tddfpt2_valence_state_release
END MODULE rixs_types

6
src/xas_tdp_atom.F
View file

@ -1919,7 +1919,11 @@ CONTAINS
do_sf = xas_tdp_control%do_spin_flip
! Get some info on the functionals
xc_functionals => section_vals_get_subs_vals(input, "DFT%XAS_TDP%KERNEL%XC_FUNCTIONAL")
IF (qs_env%do_rixs) THEN
xc_functionals => section_vals_get_subs_vals(input, "PROPERTIES%RIXS%XAS_TDP%KERNEL%XC_FUNCTIONAL")
ELSE
xc_functionals => section_vals_get_subs_vals(input, "DFT%XAS_TDP%KERNEL%XC_FUNCTIONAL")
END IF
! ask for lsd in any case
needs = xc_functionals_get_needs(xc_functionals, lsd=.TRUE., calc_potential=.TRUE.)
do_gga = needs%drho_spin !because either LDA or GGA, and the former does not need gradient

77
src/xas_tdp_methods.F
View file

@ -121,6 +121,7 @@ MODULE xas_tdp_methods
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_scf_types, ONLY: ot_method_nr
USE rixs_types, ONLY: rixs_env_type
USE util, ONLY: get_limit,&
locate,&
sort_unique
@ -158,25 +159,34 @@ CONTAINS
! **************************************************************************************************
!> \brief Driver for XAS TDDFT calculations.
!> \param qs_env the inherited qs_environment
!> \param rixs_env ...
!> \author AB
!> \note Empty for now...
! **************************************************************************************************
SUBROUTINE xas_tdp(qs_env)
SUBROUTINE xas_tdp(qs_env, rixs_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(rixs_env_type), OPTIONAL, POINTER :: rixs_env
CHARACTER(len=*), PARAMETER :: routineN = 'xas_tdp'
CHARACTER(default_string_length) :: rst_filename
INTEGER :: handle, n_rep, output_unit
LOGICAL :: do_restart
LOGICAL :: do_restart, do_rixs
TYPE(section_vals_type), POINTER :: xas_tdp_section
CALL timeset(routineN, handle)
! Logger initialization and XAS TDP banner printing
NULLIFY (xas_tdp_section)
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
! check if subroutine is called as part of rixs calculation
CALL get_qs_env(qs_env, do_rixs=do_rixs)
IF (do_rixs) THEN
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS%XAS_TDP")
ELSE
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
END IF
output_unit = cp_logger_get_default_io_unit()
IF (output_unit > 0) THEN
@ -211,7 +221,11 @@ CONTAINS
! or run the core XAS_TDP routine if not
ELSE
CALL xas_tdp_core(xas_tdp_section, qs_env)
IF (PRESENT(rixs_env)) THEN
CALL xas_tdp_core(xas_tdp_section, qs_env, rixs_env)
ELSE
CALL xas_tdp_core(xas_tdp_section, qs_env)
END IF
END IF
IF (output_unit > 0) THEN
@ -229,18 +243,20 @@ CONTAINS
!> \brief The core workflow of the XAS_TDP method
!> \param xas_tdp_section the input values for XAS_TDP
!> \param qs_env ...
!> \param rixs_env ...
! **************************************************************************************************
SUBROUTINE xas_tdp_core(xas_tdp_section, qs_env)
SUBROUTINE xas_tdp_core(xas_tdp_section, qs_env, rixs_env)
TYPE(section_vals_type), POINTER :: xas_tdp_section
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(rixs_env_type), OPTIONAL, POINTER :: rixs_env
CHARACTER(LEN=default_string_length) :: kind_name
INTEGER :: batch_size, bo(2), current_state_index, iat, iatom, ibatch, ikind, ispin, istate, &
nbatch, nex_atom, output_unit, tmp_index
INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_atoms, ex_atoms_of_kind
INTEGER, DIMENSION(:), POINTER :: atoms_of_kind
LOGICAL :: do_os, end_of_batch, unique
LOGICAL :: do_os, do_rixs, end_of_batch, unique
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
@ -262,8 +278,12 @@ CONTAINS
WRITE (UNIT=output_unit, FMT="(/,T3,A)") &
"# Create and initialize the XAS_TDP environment"
END IF
CALL get_qs_env(qs_env, dft_control=dft_control)
CALL xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
CALL get_qs_env(qs_env, dft_control=dft_control, do_rixs=do_rixs)
IF (PRESENT(rixs_env)) THEN
CALL xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env, rixs_env)
ELSE
CALL xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
END IF
CALL print_info(output_unit, xas_tdp_control, qs_env)
IF (output_unit > 0) THEN
@ -508,8 +528,7 @@ CONTAINS
IF (xas_tdp_control%do_gw2x) CALL print_xps(current_state, xas_tdp_env, xas_tdp_control, qs_env)
! Free some unneeded attributes of current_state
CALL free_ds_memory(current_state)
IF (.NOT. do_rixs) CALL free_ds_memory(current_state) ! donor-state will be cleaned in rixs
current_state_index = current_state_index + 1
NULLIFY (current_state)
@ -542,7 +561,7 @@ CONTAINS
END IF
! Clean-up
CALL xas_tdp_env_release(xas_tdp_env)
IF (.NOT. do_rixs) CALL xas_tdp_env_release(xas_tdp_env) ! is released at the end of rixs
CALL xas_tdp_control_release(xas_tdp_control)
END SUBROUTINE xas_tdp_core
@ -552,12 +571,14 @@ CONTAINS
!> \param xas_tdp_env the environment type to initialize
!> \param xas_tdp_control the control type to initialize
!> \param qs_env the inherited qs environement type
!> \param rixs_env ...
! **************************************************************************************************
SUBROUTINE xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
SUBROUTINE xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env, rixs_env)
TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(rixs_env_type), OPTIONAL, POINTER :: rixs_env
CHARACTER(LEN=default_string_length) :: kind_name
INTEGER :: at_ind, i, ispin, j, k, kind_ind, &
@ -566,7 +587,7 @@ CONTAINS
nex_kinds, nmatch, nspins
INTEGER, DIMENSION(2) :: homo, n_mo, n_moloc
INTEGER, DIMENSION(:), POINTER :: ind_of_kind
LOGICAL :: do_os, do_uks, unique
LOGICAL :: do_os, do_rixs, do_uks, unique
REAL(dp) :: fact
REAL(dp), DIMENSION(:), POINTER :: mo_evals
TYPE(admm_type), POINTER :: admm_env
@ -590,9 +611,14 @@ CONTAINS
NULLIFY (mo_coeff, matrix_ks, admm_env, dummy_section)
! XAS TDP control type initialization
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
CALL get_qs_env(qs_env, dft_control=dft_control, do_rixs=do_rixs)
IF (do_rixs) THEN
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS%XAS_TDP")
ELSE
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
END IF
CALL get_qs_env(qs_env, dft_control=dft_control)
CALL xas_tdp_control_create(xas_tdp_control)
CALL read_xas_tdp_control(xas_tdp_control, xas_tdp_section)
@ -603,7 +629,11 @@ CONTAINS
do_os = do_uks .OR. xas_tdp_control%do_roks
! XAS TDP environment type initialization
CALL xas_tdp_env_create(xas_tdp_env)
IF (PRESENT(rixs_env)) THEN
xas_tdp_env => rixs_env%core_state
ELSE
CALL xas_tdp_env_create(xas_tdp_env)
END IF
! Retrieving the excited atoms indices and correspondig state types
IF (xas_tdp_control%define_excited == xas_tdp_by_index) THEN
@ -1483,7 +1513,7 @@ CONTAINS
REAL(dp), DIMENSION(:), POINTER :: mo_evals, zeta
REAL(dp), DIMENSION(:, :), POINTER :: overlap_matrix, tmp_coeff
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_struct_type), POINTER :: eval_mat_struct, gs_struct
TYPE(cp_fm_struct_type), POINTER :: eval_mat_struct, gs_struct, matrix_struct
TYPE(cp_fm_type) :: eval_mat, work_mat
TYPE(cp_fm_type), POINTER :: gs_coeffs, mo_coeff
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
@ -1650,8 +1680,21 @@ CONTAINS
ALLOCATE (donor_state%gs_coeffs)
CALL cp_fm_create(donor_state%gs_coeffs, gs_struct)
IF (.NOT. ASSOCIATED(xas_tdp_env%mo_coeff)) THEN
ALLOCATE (xas_tdp_env%mo_coeff(nspins))
END IF
DO ispin = 1, nspins
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
! check if mo_coeff is copied before for another donor_state
IF (.NOT. ASSOCIATED(xas_tdp_env%mo_coeff(ispin)%local_data)) THEN
! copy mo_coeff
CALL cp_fm_get_info(matrix=mo_coeff, &
matrix_struct=matrix_struct)
CALL cp_fm_create(xas_tdp_env%mo_coeff(ispin), matrix_struct)
CALL cp_fm_to_fm(mo_coeff, xas_tdp_env%mo_coeff(ispin))
END IF
DO i = 1, n_states
CALL cp_fm_to_fm_submat(msource=mo_coeff, mtarget=donor_state%gs_coeffs, nrow=nao, &
ncol=1, s_firstrow=1, s_firstcol=my_mos(i, ispin), &

10
src/xas_tdp_types.F
View file

@ -203,6 +203,8 @@ MODULE xas_tdp_types
TYPE(donor_state_type), DIMENSION(:), &
POINTER :: donor_states => NULL()
INTEGER, DIMENSION(:, :, :), POINTER :: mos_of_ex_atoms => NULL()
TYPE(cp_fm_type), DIMENSION(:), &
POINTER :: mo_coeff => NULL()
TYPE(qs_loc_env_type), POINTER :: qs_loc_env => NULL()
REAL(dp), DIMENSION(:, :), POINTER :: ri_inv_coul => NULL()
REAL(dp), DIMENSION(:, :), POINTER :: ri_inv_ex => NULL()
@ -230,6 +232,7 @@ MODULE xas_tdp_types
TYPE(dbcsr_p_type), DIMENSION(:), &
POINTER :: fock_matrix => NULL()
TYPE(cp_fm_type), POINTER :: lumo_coeffs => NULL()
INTEGER :: nvirt = 0
END TYPE xas_tdp_env_type
!> *************************************************************************************************
@ -735,6 +738,7 @@ CONTAINS
NULLIFY (xas_tdp_env%donor_states)
NULLIFY (xas_tdp_env%qs_loc_env)
NULLIFY (xas_tdp_env%mos_of_ex_atoms)
NULLIFY (xas_tdp_env%mo_coeff)
NULLIFY (xas_tdp_env%ri_inv_coul)
NULLIFY (xas_tdp_env%ri_inv_ex)
NULLIFY (xas_tdp_env%opt_dist2d_coul)
@ -791,6 +795,12 @@ CONTAINS
IF (ASSOCIATED(xas_tdp_env%mos_of_ex_atoms)) THEN
DEALLOCATE (xas_tdp_env%mos_of_ex_atoms)
END IF
IF (ASSOCIATED(xas_tdp_env%mo_coeff)) THEN
DO i = 1, SIZE(xas_tdp_env%mo_coeff)
CALL cp_fm_release(xas_tdp_env%mo_coeff(i))
END DO
DEALLOCATE (xas_tdp_env%mo_coeff)
END IF
IF (ASSOCIATED(xas_tdp_env%ri_inv_coul)) THEN
DEALLOCATE (xas_tdp_env%ri_inv_coul)
END IF

1
src/xas_tdp_utils.F
View file

@ -670,6 +670,7 @@ CONTAINS
END IF
xas_tdp_env%nvirt = nevals
ALLOCATE (lr_evals(nevals))
lr_evals(:) = tmp_evals(first_ex:first_ex + nevals - 1)

82
tests/QS/regtest-rixs/CO-PBE0.inp Normal file
View file

@ -0,0 +1,82 @@
&GLOBAL
PRINT_LEVEL low
PROJECT CO-PBE0
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
&DFT
AUTO_BASIS RI_XAS SMALL
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GAPW
&END QS
&XC
&HF
FRACTION 0.25
&END HF
&XC_FUNCTIONAL
&GGA_C_PBE
&END GGA_C_PBE
&GGA_X_PBE
SCALE 0.75
&END GGA_X_PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&RIXS
&TDDFPT
KERNEL FULL
MAX_ITER 5
NSTATES 3
&END TDDFPT
&XAS_TDP
GRID O 60 120
GRID C 60 120
TDA T
&DONOR_STATES
ATOM_LIST 1 2
DEFINE_EXCITED BY_INDEX
N_SEARCH 2
STATE_TYPES 1s 1s
&END DONOR_STATES
&KERNEL
&EXACT_EXCHANGE
FRACTION 0.25
&END EXACT_EXCHANGE
&XC_FUNCTIONAL
&GGA_C_PBE
&END GGA_C_PBE
&GGA_X_PBE
SCALE 0.75
&END GGA_X_PBE
&END XC_FUNCTIONAL
&END KERNEL
&END XAS_TDP
&END RIXS
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
C 0.00000 0.00000 1.12832
O 0.00000 0.00000 0.00000
&END COORD
&KIND C
BASIS_SET 3-21G*
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET 3-21G*
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL

84
tests/QS/regtest-rixs/H2O-PBE0.inp Normal file
View file

@ -0,0 +1,84 @@
&GLOBAL
PRINT_LEVEL low
PROJECT H2O-PBE0
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
&DFT
AUTO_BASIS RI_XAS MEDIUM
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GAPW
&END QS
&XC
&HF
FRACTION 0.25
&END HF
&XC_FUNCTIONAL
&GGA_C_PBE
&END GGA_C_PBE
&GGA_X_PBE
SCALE 0.75
&END GGA_X_PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&RIXS
&TDDFPT
KERNEL FULL
MAX_ITER 5
NSTATES 3
&END TDDFPT
&XAS_TDP
GRID O 60 120
TDA T
&DONOR_STATES
ATOM_LIST 1
DEFINE_EXCITED BY_INDEX
LOCALIZE
N_SEARCH 1
STATE_TYPES 1s
&END DONOR_STATES
&KERNEL
RI_REGION 2
&EXACT_EXCHANGE
FRACTION 0.25
&END EXACT_EXCHANGE
&XC_FUNCTIONAL
&GGA_C_PBE
&END GGA_C_PBE
&GGA_X_PBE
SCALE 0.75
&END GGA_X_PBE
&END XC_FUNCTIONAL
&END KERNEL
&END XAS_TDP
&END RIXS
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
O 0.00000 0.00000 0.11779
H 0.00000 0.75545 -0.47116
H 0.00000 -0.75545 -0.47116
&END COORD
&KIND H
BASIS_SET 3-21G*
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET 3-21G*
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL

63
tests/QS/regtest-rixs/He-LDA_RIbas.inp Normal file
View file

@ -0,0 +1,63 @@
&GLOBAL
PRINT_LEVEL low
PROJECT He-LDA_RIbas
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
&DFT
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
BASIS_SET_FILE_NAME BASIS_def2_QZVP_RI_ALL
POTENTIAL_FILE_NAME POTENTIAL
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GAPW
&END QS
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&RIXS
&TDDFPT
KERNEL FULL
MAX_ITER 5
NSTATES 3
&END TDDFPT
&XAS_TDP
GRID He 150 250
TDA T
&DONOR_STATES
ATOM_LIST 1
DEFINE_EXCITED BY_INDEX
N_SEARCH 1
STATE_TYPES 1s
&END DONOR_STATES
&KERNEL
&XC_FUNCTIONAL
&LDA_XC_TETER93
&END LDA_XC_TETER93
&END XC_FUNCTIONAL
&END KERNEL
&END XAS_TDP
&END RIXS
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
He 2.5 2.5 2.5
&END COORD
&KIND He
BASIS_SET def2-QZVP
BASIS_SET RI_XAS RI-5Z
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL

64
tests/QS/regtest-rixs/Ne-LDA-e_range.inp Normal file
View file

@ -0,0 +1,64 @@
&GLOBAL
PRINT_LEVEL low
PROJECT Ne-LDA-e_range
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
&DFT
AUTO_BASIS RI_XAS MEDIUM
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GAPW
&END QS
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&RIXS
&TDDFPT
KERNEL FULL
MAX_ITER 5
NSTATES 3
&END TDDFPT
&XAS_TDP
ENERGY_RANGE 40.0
GRID Ne 150 250
TDA T
&DONOR_STATES
ATOM_LIST 1
DEFINE_EXCITED BY_INDEX
N_SEARCH 1
STATE_TYPES 1s
&END DONOR_STATES
&KERNEL
RI_REGION 2
&XC_FUNCTIONAL
&LDA_XC_TETER93
&END LDA_XC_TETER93
&END XC_FUNCTIONAL
&END KERNEL
&END XAS_TDP
&END RIXS
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
Ne 2.5 2.5 2.5
&END COORD
&KIND Ne
BASIS_SET aug-cc-pVDZ
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL

98
tests/QS/regtest-rixs/SiH4-PBE0-admm-pseudo.inp Normal file
View file

@ -0,0 +1,98 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT SiH4-PBE0-admm-pseudo
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
&DFT
AUTO_BASIS RI_XAS MEDIUM
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
BASIS_SET_FILE_NAME BASIS_MOLOPT
BASIS_SET_FILE_NAME BASIS_ADMM
POTENTIAL_FILE_NAME POTENTIAL
&AUXILIARY_DENSITY_MATRIX_METHOD
ADMM_PURIFICATION_METHOD MO_DIAG
METHOD BASIS_PROJECTION
&END AUXILIARY_DENSITY_MATRIX_METHOD
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GAPW
&END QS
&SCF
&OT
ENERGY_GAP 0.01
PRECONDITIONER FULL_ALL
&END OT
&END SCF
&XC
&HF
FRACTION 0.25
&END HF
&XC_FUNCTIONAL
&GGA_C_PBE
&END GGA_C_PBE
&GGA_X_PBE
SCALE 0.75
&END GGA_X_PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&RIXS
&TDDFPT
KERNEL FULL
MAX_ITER 5
NSTATES 2
&END TDDFPT
&XAS_TDP
GRID Si 150 150
TDA T
&DONOR_STATES
DEFINE_EXCITED BY_KIND
KIND_LIST Si
N_SEARCH 5
STATE_TYPES 1s
&END DONOR_STATES
&KERNEL
&EXACT_EXCHANGE
FRACTION 0.25
&END EXACT_EXCHANGE
&XC_FUNCTIONAL
&GGA_C_PBE
&END GGA_C_PBE
&GGA_X_PBE
SCALE 0.75
&END GGA_X_PBE
&END XC_FUNCTIONAL
&END KERNEL
&END XAS_TDP
&END RIXS
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
Si 0.000000 0.000000 0.000000
H 0.853686 0.853686 0.853686
H -0.853686 -0.853686 0.853686
H -0.853686 0.853686 -0.853686
H 0.853686 -0.853686 -0.853686
&END COORD
&KIND Si
BASIS_SET Ahlrichs-def2-SVP
BASIS_SET AUX_FIT Ahlrichs-def2-SVP
POTENTIAL ALL
&END KIND
&KIND H
BASIS_SET SZV-MOLOPT-SR-GTH
BASIS_SET AUX_FIT FIT3
POTENTIAL GTH-PBE
&END KIND
&END SUBSYS
&END FORCE_EVAL

7
tests/QS/regtest-rixs/TEST_FILES.toml Normal file
View file

@ -0,0 +1,7 @@
#Testing the RIXS method and the keyword combinations it involves
"CO-PBE0.inp" = [{matcher="M088", tol=1e-07, ref=519.422017}]
"H2O-PBE0.inp" = [{matcher="M088", tol=1e-07, ref=519.371583}]
"He-LDA_RIbas.inp" = [{matcher="M088", tol=1e-07, ref=25.875010}]
"Ne-LDA-e_range.inp" = [{matcher="M088", tol=1e-08, ref=829.640319}]
"SiH4-PBE0-admm-pseudo.inp" = [{matcher="M088", tol=1e-08, ref=1802.938610}]
#EOF

1
tests/TEST_DIRS
View file

@ -55,6 +55,7 @@ SIRIUS/regtest-1 sirius
Fist/regtest-7-2
QS/regtest-gapw_xc libint
QS/regtest-xastdp libint libxc
QS/regtest-rixs libint libxc
QS/regtest-pao-1
QS/regtest-dcdft-force libint libxc
QS/regtest-admm-gapw libint

230
tools/rixs_postprocess.py Executable file
View file

@ -0,0 +1,230 @@
#!/usr/bin/env python3
# author: BSG
import numpy as np
import matplotlib.pyplot as plt
import argparse
def parse_args():
parser = argparse.ArgumentParser(description="Run RIXS convolution script.")
parser.add_argument(
"--filename", type=str, help="Input data file from cp2k (*.rixs)"
)
parser.add_argument("--v_states", type=int, help="Number of valence excited states")
parser.add_argument(
"--gamma", type=float, default=0.075, help="Broadening width in eV"
)
parser.add_argument(
"--kappa", type=float, default=60.0, help="Scattering angle in degrees"
)
return parser.parse_args()
# hbar = 6.6260689 * 10**-34 / (2 * np.pi)
hbar = 1 # a.u.
dirs = ["x", "y", "z"]
ndirs = len(dirs)
def get_sigma_tensor(f_state, w_in, w_out, w_0i, w_if, mu_0i, mu_if, Gamma):
"""
Calculate absorption cross-section tensor of final valence state f:
sigma_{alpha beta gamma delta}^f(w)
Parameters:
- f_index: int, index of the final valence state
- w_in: float, incident photon frequency
- w_out: float, outgoing photon frequency
- w_0i: array (c_states,)
- w_if: array (v_states,)
- mu_0i: array (3) dipole moments (i->0)
- mu_if: array (c_states, v_states, 3) dipole moments (i->f)
- Gamma: float, broadening parameter
Returns:
- sigma: array (3,3,3,3) tensor
"""
sigma = np.zeros((ndirs, ndirs, ndirs, ndirs), dtype=np.float64)
c_states = len(w_0i)
for i in range(c_states):
denom = (w_in - w_0i[i]) ** 2 + Gamma**2
prefac = (w_if[f] ** 2 * w_0i[i] ** 2) / denom
for alpha in range(ndirs):
for beta in range(ndirs):
for gamma in range(ndirs):
for delta in range(ndirs):
sigma[alpha, beta, gamma, delta] += (
mu_if[i, f_state, alpha]
* mu_if[i, f_state, beta]
* mu_0i[i, gamma]
* mu_0i[i, delta]
)
sigma *= prefac
return sigma
def parse_spectrum_file(filename, v_states):
"""
Parameters:
- filename: str, path to the spectrum file
- v_states: int, number of valence excited states
"""
w_0i = []
mu_0i = []
w_if = []
mu_if = []
w_f0 = []
with open(filename, "r") as f:
lines = f.readlines()
in_core_block = False
in_valence_block = False
count = 0
for line in lines:
line = line.strip()
if "Excitation from ground-state" in line:
in_core_block = True
in_valence_block = False
continue
elif "Emission from core-excited state" in line:
in_core_block = False
in_valence_block = True
continue
elif line.startswith("=") or line.startswith("w_0i") or not line:
continue
tokens = line.split()
if in_core_block and len(tokens) >= 5:
e0i, mu_x, mu_y, mu_z = map(float, tokens[0:4])
w_0i.append(e0i / hbar)
mu_0i.append([mu_x, mu_y, mu_z])
elif in_valence_block and len(tokens) >= 6:
count += 1
ei, ef = map(float, tokens[0:2])
mu_x, mu_y, mu_z = map(float, tokens[2:5])
mu_if.append((mu_x, mu_y, mu_z))
w_if.append((ef - ei) / hbar)
if count <= v_states:
w_f0.append(ef / hbar)
xas_w_0i = np.array(w_0i)
xas_mu_0i = np.array(mu_0i)
rixs_w_if = np.array(w_if)
tddft_w_f0 = np.array(w_f0)
c_states = len(mu_0i)
# sanity check
expected_total = c_states * v_states
if len(mu_if) != expected_total:
raise ValueError(
f"Expected {expected_total} spectrum lines, but found {len(mu_if)}"
)
rixs_mu_if = np.zeros((c_states, v_states, 3))
for idx, (mx, my, mz) in enumerate(mu_if):
i = idx // v_states
f = idx % v_states
if i < c_states and f < v_states:
rixs_mu_if[i, f] = [mx, my, mz]
return xas_w_0i, xas_mu_0i, rixs_w_if, rixs_mu_if, tddft_w_f0
def convolve_sigma_averaged_tensor(
sigma_tensor, v_states, w_f0, w_0i, w_in_window, kappa_rad, Gamma
):
"""
Compute the RIXS spectrum with custom broadening.
Parameters:
- sigma_tensor: shape (v_states, 3, 3, 3, 3)
- v_states: int, number of valence excited states
- w_f0: (v_states,) energy loss values (valence excitation energies)
- w_0i: (N,) array of incident photon energies (XAS axis)
- kappa_rad: float, scattering angle in radians
- Gamma: float, broadening width
Returns:
- 2D array representing the RIXS spectrum
"""
ndirs = 3
A = 3 + np.cos(kappa_rad) ** 2
B = 0.5 * (1 - 3 * np.cos(kappa_rad) ** 2)
w_out_window = np.linspace(0, np.max(w_f0), len(w_in_window))
N_incident = len(w_in_window)
N_loss = len(w_out_window)
spectrum = np.zeros((N_incident, N_loss))
W_IN, W_LOSS = np.meshgrid(w_in_window, w_out_window, indexing="ij")
W_OUT = W_IN - W_LOSS
for f in range(v_states):
denom = (W_IN - w_0i) ** 2 + Gamma**2
prefac = (w_if[f] ** 2 * w_0i**2) / denom
Omega = W_OUT + w_f0[f] - W_IN
Delta = Gamma / np.pi / (Omega**2 + Gamma**2)
val = 0.0
for alpha in range(ndirs):
for beta in range(ndirs):
val += A * sigma_tensor[f, alpha, alpha, beta, beta]
val += B * (
sigma_tensor[f, alpha, beta, alpha, beta]
+ sigma_tensor[f, alpha, beta, beta, alpha]
)
spectrum += prefac * val * Delta / 30.0
return spectrum / np.max(spectrum)
if __name__ == "__main__":
args = parse_args()
filename = args.filename
v_states = args.v_states
kappa_deg = args.kappa
kappa_rad = np.radians(kappa_deg)
Gamma = args.gamma
molecule = filename.split(".")[0]
w_0i, mu_0i, w_if, mu_if, w_f0 = parse_spectrum_file(filename, v_states)
print(np.shape(w_if), np.shape(mu_if), np.shape(w_f0))
w_in = w_0i[0]
w_out = w_if[0]
sigma_tensor = np.zeros((v_states, 3, 3, 3, 3), dtype=np.float64)
for f in range(v_states):
sigma_tensor[f] = get_sigma_tensor(
f, w_in, w_out, w_0i, w_if, mu_0i, mu_if, Gamma
)
w_in = w_0i[0] # incident photon energy = first xas energy
e_window = 0.2
n_points = 50
w_in_window = np.linspace(w_in - e_window, w_in + e_window, n_points)
spectrum_2d = convolve_sigma_averaged_tensor(
sigma_tensor, v_states, w_f0, w_0i[0], w_in_window, kappa_rad, Gamma
)
project_name = filename.split(".")[0]
np.savetxt(f"{project_name}_rixs.dat", spectrum_2d, fmt="%.6f")