From bf603df4a084d7f4ea8f31e645a7621a95ed892c Mon Sep 17 00:00:00 2001 From: BelizSertcan <48097107+BelizSertcan@users.noreply.github.com> Date: Tue, 22 Jul 2025 16:39:38 +0200 Subject: [PATCH] Resonant Inelastic X-ray Spectroscopy (RIXS) Module --- src/CMakeLists.txt | 2 + src/common/bibliography.F | 9 +- src/cp_control_types.F | 55 +- src/cp_control_utils.F | 26 +- src/header.F | 21 +- src/input_cp2k_properties_dft.F | 51 +- src/qs_energy_utils.F | 5 + src/qs_environment.F | 26 +- src/qs_environment_types.F | 53 +- src/qs_scf_initialization.F | 2 + src/qs_tddfpt2_methods.F | 55 +- src/rixs_methods.F | 470 ++++++++++++++++++ src/rixs_types.F | 135 +++++ src/xas_tdp_atom.F | 6 +- src/xas_tdp_methods.F | 77 ++- src/xas_tdp_types.F | 10 + src/xas_tdp_utils.F | 1 + tests/QS/regtest-rixs/CO-PBE0.inp | 82 +++ tests/QS/regtest-rixs/H2O-PBE0.inp | 84 ++++ tests/QS/regtest-rixs/He-LDA_RIbas.inp | 63 +++ tests/QS/regtest-rixs/Ne-LDA-e_range.inp | 64 +++ .../QS/regtest-rixs/SiH4-PBE0-admm-pseudo.inp | 98 ++++ tests/QS/regtest-rixs/TEST_FILES.toml | 7 + tests/TEST_DIRS | 1 + tools/rixs_postprocess.py | 230 +++++++++ 25 files changed, 1571 insertions(+), 62 deletions(-) create mode 100644 src/rixs_methods.F create mode 100644 src/rixs_types.F create mode 100644 tests/QS/regtest-rixs/CO-PBE0.inp create mode 100644 tests/QS/regtest-rixs/H2O-PBE0.inp create mode 100644 tests/QS/regtest-rixs/He-LDA_RIbas.inp create mode 100644 tests/QS/regtest-rixs/Ne-LDA-e_range.inp create mode 100644 tests/QS/regtest-rixs/SiH4-PBE0-admm-pseudo.inp create mode 100644 tests/QS/regtest-rixs/TEST_FILES.toml create mode 100755 tools/rixs_postprocess.py diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f3f2c38342..1334bb6da5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/common/bibliography.F b/src/common/bibliography.F index 602aadba57..19b82df8cc 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -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 diff --git a/src/cp_control_types.F b/src/cp_control_types.F index 8b2a13cf1b..c7eb4f8719 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -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 ... diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 6e6c319f97..f9adc88a7d 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -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 ... diff --git a/src/header.F b/src/header.F index 75d83ce0a6..47d7ed507e 100644 --- a/src/header.F +++ b/src/header.F @@ -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 diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index d91a13a09b..8a97a2df9e 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -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 diff --git a/src/qs_energy_utils.F b/src/qs_energy_utils.F index 8fd6da0822..448acc8800 100644 --- a/src/qs_energy_utils.F +++ b/src/qs_energy_utils.F @@ -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") diff --git a/src/qs_environment.F b/src/qs_environment.F index 277347c683..707d514698 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -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) diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index 15ca1bde2a..4f5caab124 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -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 diff --git a/src/qs_scf_initialization.F b/src/qs_scf_initialization.F index fe51cf35f9..03f0db452a 100644 --- a/src/qs_scf_initialization.F +++ b/src/qs_scf_initialization.F @@ -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 diff --git a/src/qs_tddfpt2_methods.F b/src/qs_tddfpt2_methods.F index f29c00b6b2..0d3fbbed54 100644 --- a/src/qs_tddfpt2_methods.F +++ b/src/qs_tddfpt2_methods.F @@ -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) diff --git a/src/rixs_methods.F b/src/rixs_methods.F new file mode 100644 index 0000000000..14df612844 --- /dev/null +++ b/src/rixs_methods.F @@ -0,0 +1,470 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! 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 diff --git a/src/rixs_types.F b/src/rixs_types.F new file mode 100644 index 0000000000..b8107e899d --- /dev/null +++ b/src/rixs_types.F @@ -0,0 +1,135 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! 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 diff --git a/src/xas_tdp_atom.F b/src/xas_tdp_atom.F index 8970b5ad80..116bcd1603 100644 --- a/src/xas_tdp_atom.F +++ b/src/xas_tdp_atom.F @@ -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 diff --git a/src/xas_tdp_methods.F b/src/xas_tdp_methods.F index eff58370a6..61998c496c 100644 --- a/src/xas_tdp_methods.F +++ b/src/xas_tdp_methods.F @@ -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), & diff --git a/src/xas_tdp_types.F b/src/xas_tdp_types.F index 4d32bcbb74..827fbf7ef8 100644 --- a/src/xas_tdp_types.F +++ b/src/xas_tdp_types.F @@ -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 diff --git a/src/xas_tdp_utils.F b/src/xas_tdp_utils.F index 4e90b9aa96..535323c16c 100644 --- a/src/xas_tdp_utils.F +++ b/src/xas_tdp_utils.F @@ -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) diff --git a/tests/QS/regtest-rixs/CO-PBE0.inp b/tests/QS/regtest-rixs/CO-PBE0.inp new file mode 100644 index 0000000000..bddba1eb7f --- /dev/null +++ b/tests/QS/regtest-rixs/CO-PBE0.inp @@ -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 diff --git a/tests/QS/regtest-rixs/H2O-PBE0.inp b/tests/QS/regtest-rixs/H2O-PBE0.inp new file mode 100644 index 0000000000..8ff82ec0ba --- /dev/null +++ b/tests/QS/regtest-rixs/H2O-PBE0.inp @@ -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 diff --git a/tests/QS/regtest-rixs/He-LDA_RIbas.inp b/tests/QS/regtest-rixs/He-LDA_RIbas.inp new file mode 100644 index 0000000000..e6deeda346 --- /dev/null +++ b/tests/QS/regtest-rixs/He-LDA_RIbas.inp @@ -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 diff --git a/tests/QS/regtest-rixs/Ne-LDA-e_range.inp b/tests/QS/regtest-rixs/Ne-LDA-e_range.inp new file mode 100644 index 0000000000..f1f7df3f5e --- /dev/null +++ b/tests/QS/regtest-rixs/Ne-LDA-e_range.inp @@ -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 diff --git a/tests/QS/regtest-rixs/SiH4-PBE0-admm-pseudo.inp b/tests/QS/regtest-rixs/SiH4-PBE0-admm-pseudo.inp new file mode 100644 index 0000000000..4a20cdf31b --- /dev/null +++ b/tests/QS/regtest-rixs/SiH4-PBE0-admm-pseudo.inp @@ -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 diff --git a/tests/QS/regtest-rixs/TEST_FILES.toml b/tests/QS/regtest-rixs/TEST_FILES.toml new file mode 100644 index 0000000000..acc8e927d9 --- /dev/null +++ b/tests/QS/regtest-rixs/TEST_FILES.toml @@ -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 diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index d8a3a51d63..39bcc9b294 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -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 diff --git a/tools/rixs_postprocess.py b/tools/rixs_postprocess.py new file mode 100755 index 0000000000..1ef5d89e45 --- /dev/null +++ b/tools/rixs_postprocess.py @@ -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")