LCOV - code coverage report
Current view: top level - src - rixs_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:8ebf9ad) Lines: 99.1 % 220 218
Test Date: 2026-01-22 06:43:13 Functions: 100.0 % 4 4

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Methods for Resonant Inelastic XRAY Scattering (RIXS) calculations
      10              : !> \author BSG (02.2025)
      11              : ! **************************************************************************************************
      12              : MODULE rixs_methods
      13              :    USE bibliography,                    ONLY: VazdaCruz2021,&
      14              :                                               cite_reference
      15              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      16              :    USE cp_control_types,                ONLY: dft_control_type,&
      17              :                                               rixs_control_create,&
      18              :                                               rixs_control_release,&
      19              :                                               rixs_control_type
      20              :    USE cp_control_utils,                ONLY: read_rixs_control
      21              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type,&
      22              :                                               dbcsr_type
      23              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      24              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      25              :                                               cp_fm_struct_release,&
      26              :                                               cp_fm_struct_type
      27              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      28              :                                               cp_fm_get_info,&
      29              :                                               cp_fm_get_submatrix,&
      30              :                                               cp_fm_release,&
      31              :                                               cp_fm_to_fm,&
      32              :                                               cp_fm_to_fm_submat,&
      33              :                                               cp_fm_type
      34              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      35              :                                               cp_logger_get_default_io_unit,&
      36              :                                               cp_logger_type
      37              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      38              :                                               cp_print_key_unit_nr
      39              :    USE header,                          ONLY: rixs_header
      40              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      41              :                                               section_vals_type
      42              :    USE kinds,                           ONLY: dp
      43              :    USE message_passing,                 ONLY: mp_para_env_type
      44              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      45              :    USE physcon,                         ONLY: evolt
      46              :    USE qs_environment_types,            ONLY: get_qs_env,&
      47              :                                               qs_environment_type
      48              :    USE qs_tddfpt2_methods,              ONLY: tddfpt
      49              :    USE rixs_types,                      ONLY: rixs_env_create,&
      50              :                                               rixs_env_release,&
      51              :                                               rixs_env_type,&
      52              :                                               tddfpt2_valence_type
      53              :    USE xas_tdp_methods,                 ONLY: xas_tdp
      54              :    USE xas_tdp_types,                   ONLY: donor_state_type,&
      55              :                                               xas_tdp_env_type
      56              : #include "./base/base_uses.f90"
      57              : 
      58              :    IMPLICIT NONE
      59              :    PRIVATE
      60              : 
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rixs_methods'
      62              : 
      63              :    PUBLIC :: rixs, rixs_core
      64              : 
      65              : CONTAINS
      66              : 
      67              : ! **************************************************************************************************
      68              : !> \brief Driver for RIXS calculations.
      69              : !> \param qs_env the inherited qs_environment
      70              : !> \author BSG
      71              : ! **************************************************************************************************
      72              : 
      73           14 :    SUBROUTINE rixs(qs_env)
      74              : 
      75              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      76              : 
      77              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'rixs'
      78              : 
      79              :       INTEGER                                            :: handle, output_unit
      80              :       TYPE(dft_control_type), POINTER                    :: dft_control
      81              :       TYPE(section_vals_type), POINTER                   :: rixs_section, tddfp2_section, &
      82              :                                                             xas_tdp_section
      83              : 
      84           14 :       CALL timeset(routineN, handle)
      85              : 
      86           14 :       NULLIFY (rixs_section)
      87           14 :       rixs_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS")
      88           14 :       output_unit = cp_logger_get_default_io_unit()
      89              : 
      90           14 :       qs_env%do_rixs = .TRUE.
      91              : 
      92           14 :       CALL cite_reference(VazdaCruz2021)
      93              : 
      94           14 :       CALL get_qs_env(qs_env, dft_control=dft_control)
      95              : 
      96           14 :       xas_tdp_section => section_vals_get_subs_vals(rixs_section, "XAS_TDP")
      97           14 :       tddfp2_section => section_vals_get_subs_vals(rixs_section, "TDDFPT")
      98              : 
      99           14 :       CALL rixs_core(rixs_section, qs_env)
     100              : 
     101           14 :       IF (output_unit > 0) THEN
     102              :          WRITE (UNIT=output_unit, FMT="(/,(T2,A79))") &
     103            7 :             "*******************************************************************************", &
     104            7 :             "!    Normal termination of Resonant Inelastic X-RAY Scattering calculation    !", &
     105           14 :             "*******************************************************************************"
     106              :       END IF
     107              : 
     108           14 :       CALL timestop(handle)
     109              : 
     110           14 :    END SUBROUTINE rixs
     111              : 
     112              : ! **************************************************************************************************
     113              : !> \brief Perform RIXS calculation.
     114              : !> \param rixs_section ...
     115              : !> \param qs_env ...
     116              : ! **************************************************************************************************
     117           14 :    SUBROUTINE rixs_core(rixs_section, qs_env)
     118              : 
     119              :       TYPE(section_vals_type), POINTER                   :: rixs_section
     120              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     121              : 
     122              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'rixs_core'
     123              : 
     124              :       INTEGER :: ax, current_state_index, fstate, handle, iatom, ispin, istate, nactive_max, nao, &
     125              :          ncol, nex_atoms, nocc, nspins, nstates, nvirt, output_unit, td_state
     126           14 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nactive
     127              :       LOGICAL                                            :: do_sc, do_sg, roks, uks
     128           14 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: w_i0, w_if
     129           14 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: dip_block, mu_i0
     130           14 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :, :)          :: mu_if
     131              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     132              :       TYPE(cp_fm_struct_type), POINTER                   :: core_evect_struct, dip_0_struct, &
     133              :                                                             dip_f_struct, gs_coeff_struct, &
     134              :                                                             i_dip_0_struct, i_dip_f_struct
     135              :       TYPE(cp_fm_type)                                   :: dip_0
     136           14 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: core_evects, dip_f, i_dip_0, i_dip_f, &
     137           14 :                                                             state_gs_coeffs
     138           14 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: local_gs_coeffs, mo_coeffs
     139           14 :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: valence_evects
     140              :       TYPE(cp_fm_type), POINTER                          :: target_ex_coeffs
     141           14 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: dipmat, matrix_s
     142              :       TYPE(dft_control_type), POINTER                    :: dft_control
     143              :       TYPE(donor_state_type), POINTER                    :: current_state
     144              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     145              :       TYPE(rixs_control_type), POINTER                   :: rixs_control
     146              :       TYPE(rixs_env_type), POINTER                       :: rixs_env
     147              :       TYPE(tddfpt2_valence_type), POINTER                :: valence_state
     148              :       TYPE(xas_tdp_env_type), POINTER                    :: core_state
     149              : 
     150           14 :       NULLIFY (rixs_control, dft_control, rixs_env)
     151           14 :       NULLIFY (valence_state, core_state)
     152           14 :       NULLIFY (para_env, blacs_env)
     153           14 :       NULLIFY (local_gs_coeffs, mo_coeffs, valence_evects)
     154           14 :       NULLIFY (dipmat, dip_0_struct, i_dip_0_struct, dip_f_struct, i_dip_f_struct, &
     155           14 :                core_evect_struct, gs_coeff_struct)
     156              : 
     157           28 :       output_unit = cp_logger_get_default_io_unit()
     158              : 
     159              :       CALL get_qs_env(qs_env, &
     160              :                       dft_control=dft_control, &
     161              :                       matrix_s=matrix_s, &
     162              :                       para_env=para_env, &
     163           14 :                       blacs_env=blacs_env)
     164           14 :       CALL rixs_control_create(rixs_control)
     165           14 :       CALL read_rixs_control(rixs_control, rixs_section, dft_control%qs_control)
     166              : 
     167              :       ! create rixs_env
     168           14 :       CALL rixs_env_create(rixs_env)
     169              : 
     170              :       ! first, xas_tdp calculation
     171           14 :       CALL xas_tdp(qs_env, rixs_env)
     172              : 
     173           14 :       do_sg = rixs_control%xas_tdp_control%do_singlet
     174           14 :       do_sc = rixs_control%xas_tdp_control%do_spin_cons
     175              : 
     176           14 :       IF (rixs_control%xas_tdp_control%check_only) THEN
     177            0 :          CPWARN("CHECK_ONLY run for XAS_TDP requested, RIXS and TDDFPT will not be performed.")
     178              :       ELSE
     179              : 
     180              :          ! then, tddfpt calculation
     181           14 :          CALL tddfpt(qs_env, calc_forces=.FALSE., rixs_env=rixs_env)
     182              : 
     183           14 :          IF (output_unit > 0) THEN
     184            7 :             CALL rixs_header(output_unit)
     185              :          END IF
     186              : 
     187              :          ! timings for rixs only, excluding xas_tdp and tddft calls
     188           14 :          CALL timeset(routineN, handle)
     189              : 
     190           14 :          IF (do_sg) THEN ! singlet
     191              :             nspins = 1
     192            4 :          ELSE IF (do_sc) THEN ! spin-conserving
     193              :             nspins = 2
     194              :          ELSE
     195            0 :             CPABORT("RIXS only implemented for singlet and spin-conserving excitations")
     196              :          END IF
     197              : 
     198           14 :          IF (output_unit > 0) THEN
     199            7 :             IF (dft_control%uks) THEN
     200            1 :                uks = .TRUE.
     201            1 :                WRITE (UNIT=output_unit, FMT="(T2,A)") "RIXS| Unrestricted Open-Shell Kohn-Sham"
     202            6 :             ELSE IF (dft_control%roks) THEN
     203            1 :                roks = .TRUE.
     204            1 :                WRITE (UNIT=output_unit, FMT="(T2,A)") "RIXS| Restricted Open-Shell Kohn-Sham"
     205              :             END IF
     206              :          END IF
     207              : 
     208           14 :          core_state => rixs_env%core_state
     209           14 :          valence_state => rixs_env%valence_state
     210              : 
     211              :          ! gs coefficients from tddfpt
     212           14 :          mo_coeffs => valence_state%mos_active
     213              :          ! localised gs coefficients from xas_tdp
     214           14 :          local_gs_coeffs => core_state%mo_coeff
     215           14 :          valence_evects => valence_state%evects
     216              : 
     217              :          ! res tddft
     218              :          IF (.NOT. roks) THEN
     219              :             CALL cp_fm_get_info(matrix=local_gs_coeffs(1), ncol_global=nocc)
     220              :             CALL cp_fm_get_info(matrix=mo_coeffs(1), ncol_global=ncol)
     221              :             IF (ncol /= nocc) THEN
     222              :                CPABORT("RIXS with restricted space excitations NYI")
     223              :             END IF
     224              :          END IF
     225              : 
     226           14 :          IF (rixs_control%xas_tdp_control%do_loc) THEN
     227            2 :             IF (output_unit > 0) THEN
     228              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
     229            1 :                   "RIXS| Found localised XAS_TDP orbitals"
     230              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
     231            1 :                   "RIXS| Rotating TDDFPT vectors..."
     232              :             END IF
     233            2 :             CALL rotate_vectors(valence_state%evects, local_gs_coeffs, mo_coeffs, matrix_s(1)%matrix, output_unit)
     234              :          END IF
     235              : 
     236              :          ! find max nactive for open-shell cases
     237           28 :          ALLOCATE (nactive(nspins))
     238           32 :          DO ispin = 1, nspins
     239           32 :             CALL cp_fm_get_info(matrix=valence_state%mos_active(ispin), nrow_global=nao, ncol_global=nactive(ispin))
     240              :          END DO
     241           32 :          nactive_max = MAXVAL(nactive)
     242              : 
     243           14 :          nex_atoms = core_state%nex_atoms
     244           14 :          nstates = valence_state%nstates
     245           14 :          nvirt = core_state%nvirt
     246              : 
     247           14 :          IF (rixs_control%core_states > 0) THEN
     248            4 :             rixs_control%core_states = MIN(rixs_control%core_states, nvirt)
     249              :          ELSE
     250           10 :             rixs_control%core_states = nvirt
     251              :          END IF
     252              : 
     253           14 :          IF (rixs_control%valence_states > 0) THEN
     254            2 :             rixs_control%valence_states = MIN(rixs_control%valence_states, nstates)
     255              :          ELSE
     256           12 :             rixs_control%valence_states = nstates
     257              :          END IF
     258              : 
     259           14 :          IF (output_unit > 0) THEN
     260              :             WRITE (UNIT=output_unit, FMT="(T2,A,I5,A,I5)") &
     261            7 :                "RIXS| Using ", rixs_control%core_states, " core states out of ", core_state%nvirt
     262              :             WRITE (UNIT=output_unit, FMT="(T2,A,I5,A,I5,/)") &
     263            7 :                "RIXS| Using ", rixs_control%valence_states, " valence states out of ", valence_state%nstates
     264              :          END IF
     265              : 
     266           14 :          dipmat => core_state%dipmat
     267              : 
     268           78 :          ALLOCATE (core_evects(nspins), state_gs_coeffs(nspins))
     269          154 :          ALLOCATE (dip_block(1, nspins), mu_i0(4, nvirt), mu_if(4, nvirt, nstates), w_i0(nvirt), w_if(nstates))
     270           62 :          w_if(:) = valence_state%evals(:)*evolt
     271           46 :          ALLOCATE (i_dip_0(nspins))
     272           78 :          ALLOCATE (dip_f(nspins), i_dip_f(nspins))
     273              : 
     274              :          CALL cp_fm_struct_create(core_evect_struct, para_env=para_env, context=blacs_env, &
     275           14 :                                   nrow_global=nao, ncol_global=nvirt)
     276              :          CALL cp_fm_struct_create(gs_coeff_struct, para_env=para_env, context=blacs_env, &
     277           14 :                                   nrow_global=nao, ncol_global=1)
     278              : 
     279              :          ! looping over ex_atoms and ex_kinds is enough as excited atoms have to be unique
     280           14 :          current_state_index = 1
     281           30 :          DO iatom = 1, nex_atoms
     282           16 :             current_state => core_state%donor_states(current_state_index)
     283           16 :             IF (output_unit > 0) THEN
     284              :                WRITE (UNIT=output_unit, FMT="(T2,A,A,A,I3,A,A)") &
     285            8 :                   "RIXS| Calculating dipole moment from core-excited state ", &
     286            8 :                   core_state%state_type_char(current_state%state_type), " for atom ", &
     287           16 :                   current_state%at_index, " of kind ", TRIM(current_state%at_symbol)
     288              :             END IF
     289              : 
     290         1096 :             mu_i0 = 0.0_dp
     291         3660 :             mu_if = 0.0_dp
     292              : 
     293           16 :             IF (do_sg) THEN ! singlet
     294           12 :                target_ex_coeffs => current_state%sg_coeffs
     295          180 :                w_i0(:) = current_state%sg_evals(:)*evolt
     296            4 :             ELSE IF (do_sc) THEN ! spin-conserving
     297            4 :                target_ex_coeffs => current_state%sc_coeffs
     298           52 :                w_i0(:) = current_state%sc_evals(:)*evolt
     299              :             END IF
     300              : 
     301              :             ! reshape sc and sg coeffs (separate spins to columns)
     302           36 :             DO ispin = 1, nspins
     303           20 :                CALL cp_fm_create(core_evects(ispin), core_evect_struct)
     304              :                CALL cp_fm_to_fm_submat(msource=target_ex_coeffs, mtarget=core_evects(ispin), s_firstrow=1, &
     305           36 :                                        s_firstcol=(nvirt*(ispin - 1) + 1), t_firstrow=1, t_firstcol=1, nrow=nao, ncol=nvirt)
     306              :             END DO
     307           36 :             DO ispin = 1, nspins
     308           20 :                CALL cp_fm_create(state_gs_coeffs(ispin), gs_coeff_struct)
     309           16 :                IF (roks) THEN
     310              :                   ! store same coeffs for both spins, easier later on
     311              :                   CALL cp_fm_to_fm_submat(msource=current_state%gs_coeffs, mtarget=state_gs_coeffs(ispin), s_firstrow=1, &
     312           20 :                                           s_firstcol=1, t_firstrow=1, t_firstcol=1, nrow=nao, ncol=1)
     313              :                ELSE
     314              :                   CALL cp_fm_to_fm_submat(msource=current_state%gs_coeffs, mtarget=state_gs_coeffs(ispin), s_firstrow=1, &
     315              :                                           s_firstcol=ispin, t_firstrow=1, t_firstcol=1, nrow=nao, ncol=1)
     316              :                END IF
     317              :             END DO
     318              : 
     319              :             ! initialise matrices for i->0
     320              :             CALL cp_fm_struct_create(dip_0_struct, para_env=para_env, context=blacs_env, &
     321           16 :                                      nrow_global=nao, ncol_global=1)
     322           16 :             CALL cp_fm_create(dip_0, dip_0_struct)
     323              :             CALL cp_fm_struct_create(i_dip_0_struct, para_env=para_env, context=blacs_env, &
     324           16 :                                      nrow_global=nvirt, ncol_global=1)
     325           36 :             DO ispin = 1, nspins
     326           36 :                CALL cp_fm_create(i_dip_0(ispin), i_dip_0_struct)
     327              :             END DO
     328              : 
     329              :             ! initialise matrices for i->f
     330           36 :             DO ispin = 1, nspins
     331              :                CALL cp_fm_struct_create(dip_f_struct, para_env=para_env, context=blacs_env, &
     332           20 :                                         nrow_global=nao, ncol_global=nactive(ispin))
     333              :                CALL cp_fm_struct_create(i_dip_f_struct, para_env=para_env, context=blacs_env, &
     334           20 :                                         nrow_global=nvirt, ncol_global=nactive(ispin))
     335           20 :                CALL cp_fm_create(dip_f(ispin), dip_f_struct)
     336           20 :                CALL cp_fm_create(i_dip_f(ispin), i_dip_f_struct)
     337           20 :                CALL cp_fm_struct_release(i_dip_f_struct)
     338           36 :                CALL cp_fm_struct_release(dip_f_struct)
     339              :             END DO
     340              : 
     341              :             ! 0 -> i
     342           64 :             DO ax = 1, 3
     343              : 
     344              :                ! i*R*0
     345          108 :                DO ispin = 1, nspins
     346           60 :                   CALL cp_dbcsr_sm_fm_multiply(dipmat(ax)%matrix, state_gs_coeffs(ispin), dip_0, ncol=1)
     347          108 :                   CALL parallel_gemm('T', 'N', nvirt, 1, nao, 1.0_dp, core_evects(ispin), dip_0, 0.0_dp, i_dip_0(ispin))
     348              :                END DO
     349              : 
     350          562 :                DO istate = 1, rixs_control%core_states
     351         1782 :                   dip_block = 0.0_dp
     352         1140 :                   DO ispin = 1, nspins
     353              :                      CALL cp_fm_get_submatrix(fm=i_dip_0(ispin), target_m=dip_block, start_row=istate, &
     354          642 :                                               start_col=1, n_rows=1, n_cols=1)
     355         1140 :                      mu_i0(ax, istate) = mu_i0(ax, istate) + dip_block(1, 1)
     356              :                   END DO ! ispin
     357          546 :                   mu_i0(4, istate) = mu_i0(4, istate) + mu_i0(ax, istate)**2
     358              :                END DO ! istate
     359              : 
     360              :             END DO ! ax
     361              : 
     362              :             ! i -> f
     363           66 :             DO td_state = 1, rixs_control%valence_states
     364              : 
     365           50 :                IF (output_unit > 0) THEN
     366              :                   WRITE (UNIT=output_unit, FMT="(T9,A,I3,A,F10.4)") &
     367           25 :                      "to valence-excited state ", td_state, " with energy ", w_if(td_state)
     368              :                END IF
     369              : 
     370          216 :                DO ax = 1, 3
     371              : 
     372              :                   ! core_evects x dipmat x valence_evects (per spin)
     373          360 :                   DO ispin = 1, nspins
     374              :                      CALL cp_dbcsr_sm_fm_multiply(dipmat(ax)%matrix, valence_evects(ispin, td_state), dip_f(ispin), &
     375          210 :                                                   ncol=nactive(ispin))
     376              :                      CALL parallel_gemm('T', 'N', nvirt, nactive(ispin), nao, 1.0_dp, core_evects(ispin), &
     377          360 :                                         dip_f(ispin), 0.0_dp, i_dip_f(ispin))
     378              :                   END DO
     379              : 
     380         1832 :                   DO istate = 1, rixs_control%core_states
     381         9696 :                      DO fstate = 1, nactive_max
     382        31392 :                         dip_block = 0.0_dp
     383        21360 :                         DO ispin = 1, nspins
     384        19728 :                            IF (fstate <= nactive(ispin)) THEN
     385              :                               CALL cp_fm_get_submatrix(fm=i_dip_f(ispin), target_m=dip_block, start_row=istate, &
     386        10944 :                                                        start_col=fstate, n_rows=1, n_cols=1)
     387        10944 :                               mu_if(ax, istate, td_state) = mu_if(ax, istate, td_state) + dip_block(1, 1)
     388              :                            END IF
     389              :                         END DO ! ispin
     390              :                      END DO ! fstate (tddft)
     391         1782 :                      mu_if(4, istate, td_state) = mu_if(4, istate, td_state) + mu_if(ax, istate, td_state)**2
     392              :                   END DO ! istate (core)
     393              : 
     394              :                END DO ! ax
     395              : 
     396              :             END DO ! td_state
     397              : 
     398           16 :             IF (output_unit > 0) THEN
     399            8 :                WRITE (UNIT=output_unit, FMT="(/,T2,A,/)") "RIXS| Printing spectrum to file"
     400              :             END IF
     401           16 :             CALL print_rixs_to_file(current_state, mu_i0, mu_if, w_i0, w_if, rixs_env, rixs_section, rixs_control)
     402              : 
     403           16 :             current_state_index = current_state_index + 1
     404              : 
     405              :             ! cleanup
     406           36 :             DO ispin = 1, nspins
     407           20 :                CALL cp_fm_release(core_evects(ispin))
     408           20 :                CALL cp_fm_release(state_gs_coeffs(ispin))
     409           20 :                CALL cp_fm_release(i_dip_0(ispin))
     410           20 :                CALL cp_fm_release(i_dip_f(ispin))
     411           36 :                CALL cp_fm_release(dip_f(ispin))
     412              :             END DO
     413           16 :             CALL cp_fm_struct_release(i_dip_0_struct)
     414           16 :             CALL cp_fm_struct_release(dip_0_struct)
     415           46 :             CALL cp_fm_release(dip_0)
     416              : 
     417              :          END DO ! iatom
     418              : 
     419              :          NULLIFY (current_state)
     420              : 
     421              :          ! cleanup
     422           14 :          CALL cp_fm_struct_release(core_evect_struct)
     423           28 :          CALL cp_fm_struct_release(gs_coeff_struct)
     424              : 
     425              :       END IF
     426              : 
     427              :       ! more cleanup
     428           14 :       CALL rixs_control_release(rixs_control)
     429           14 :       CALL rixs_env_release(rixs_env)
     430           14 :       NULLIFY (valence_state, core_state)
     431              : 
     432           14 :       CALL timestop(handle)
     433              : 
     434           28 :    END SUBROUTINE rixs_core
     435              : 
     436              : ! **************************************************************************************************
     437              : !> \brief Rotate vectors. Returns rotated mo_occ and evects.
     438              : !> \param evects ...
     439              : !> \param mo_ref ...
     440              : !> \param mo_occ ...
     441              : !> \param overlap_matrix ...
     442              : !> \param unit_nr ...
     443              : ! **************************************************************************************************
     444              : 
     445            2 :    SUBROUTINE rotate_vectors(evects, mo_ref, mo_occ, overlap_matrix, unit_nr)
     446              :       TYPE(cp_fm_type), DIMENSION(:, :)                  :: evects
     447              :       TYPE(cp_fm_type), DIMENSION(:)                     :: mo_ref, mo_occ
     448              :       TYPE(dbcsr_type), POINTER                          :: overlap_matrix
     449              :       INTEGER                                            :: unit_nr
     450              : 
     451              :       INTEGER                                            :: ispin, istate, ncol, nrow, nspins, &
     452              :                                                             nstates
     453              :       REAL(kind=dp)                                      :: diff
     454              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     455              :       TYPE(cp_fm_struct_type), POINTER                   :: emat_struct
     456              :       TYPE(cp_fm_type)                                   :: emat, rotated_mo_coeffs, smo
     457              :       TYPE(cp_fm_type), POINTER                          :: current_evect
     458              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     459              : 
     460            2 :       NULLIFY (emat_struct, para_env, blacs_env, current_evect)
     461              : 
     462            2 :       nspins = SIZE(evects, DIM=1)
     463            4 :       DO ispin = 1, nspins
     464              : 
     465              :          CALL cp_fm_get_info(matrix=mo_occ(ispin), nrow_global=nrow, ncol_global=ncol, &
     466            2 :                              para_env=para_env, context=blacs_env)
     467            2 :          CALL cp_fm_create(smo, mo_occ(ispin)%matrix_struct)
     468              : 
     469              :          ! rotate mo_occ
     470              :          ! smo = matrix_s x mo_occ
     471            2 :          CALL cp_dbcsr_sm_fm_multiply(overlap_matrix, mo_occ(ispin), smo, ncol, alpha=1.0_dp, beta=0.0_dp)
     472              :          CALL cp_fm_struct_create(emat_struct, nrow_global=ncol, ncol_global=ncol, &
     473            2 :                                   para_env=para_env, context=blacs_env)
     474            2 :          CALL cp_fm_create(emat, emat_struct)
     475              :          ! emat = mo_ref^T x smo
     476            2 :          CALL parallel_gemm('T', 'N', ncol, ncol, nrow, 1.0_dp, mo_ref(ispin), smo, 0.0_dp, emat)
     477            2 :          CALL cp_fm_create(rotated_mo_coeffs, mo_occ(ispin)%matrix_struct)
     478              :          ! rotated_mo_coeffs = cpmos x emat
     479            2 :          CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, mo_occ(ispin), emat, 0.0_dp, rotated_mo_coeffs)
     480              : 
     481           77 :          diff = MAXVAL(ABS(rotated_mo_coeffs%local_data - mo_occ(ispin)%local_data))
     482            2 :          IF (unit_nr > 0) THEN
     483            1 :             WRITE (unit_nr, FMT="(T9,A,I2,A,F10.6,/)") "For spin ", ispin, ": Max difference between orbitals = ", diff
     484              :          END IF
     485              : 
     486            2 :          CALL cp_fm_to_fm(rotated_mo_coeffs, mo_occ(ispin))
     487              : 
     488            2 :          nstates = SIZE(evects, DIM=2)
     489            8 :          DO istate = 1, nstates
     490            2 :             ASSOCIATE (current_evect => evects(ispin, istate))
     491            6 :                CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, current_evect, emat, 0.0_dp, smo)
     492            6 :                CALL cp_fm_to_fm(smo, current_evect)
     493              :             END ASSOCIATE
     494              :          END DO
     495              : 
     496            2 :          CALL cp_fm_struct_release(emat_struct)
     497            2 :          CALL cp_fm_release(smo)
     498            2 :          CALL cp_fm_release(emat)
     499           10 :          CALL cp_fm_release(rotated_mo_coeffs)
     500              : 
     501              :       END DO ! ispin
     502              : 
     503            2 :    END SUBROUTINE rotate_vectors
     504              : 
     505              : !**************************************************************************************************
     506              : !> \brief Print RIXS spectrum.
     507              : !> \param donor_state ...
     508              : !> \param mu_i0 ...
     509              : !> \param mu_if ...
     510              : !> \param w_i0 ...
     511              : !> \param w_if ...
     512              : !> \param rixs_env ...
     513              : !> \param rixs_section ...
     514              : !> \param rixs_control ...
     515              : ! **************************************************************************************************
     516           16 :    SUBROUTINE print_rixs_to_file(donor_state, mu_i0, mu_if, w_i0, w_if, &
     517              :                                  rixs_env, rixs_section, rixs_control)
     518              : 
     519              :       TYPE(donor_state_type), POINTER                    :: donor_state
     520              :       REAL(dp), DIMENSION(:, :)                          :: mu_i0
     521              :       REAL(dp), DIMENSION(:, :, :)                       :: mu_if
     522              :       REAL(dp), DIMENSION(:)                             :: w_i0, w_if
     523              :       TYPE(rixs_env_type), POINTER                       :: rixs_env
     524              :       TYPE(section_vals_type), POINTER                   :: rixs_section
     525              :       TYPE(rixs_control_type), POINTER                   :: rixs_control
     526              : 
     527              :       INTEGER                                            :: f, i, output_unit, rixs_unit
     528              :       TYPE(cp_logger_type), POINTER                      :: logger
     529              : 
     530           16 :       NULLIFY (logger)
     531           16 :       logger => cp_get_default_logger()
     532              : 
     533              :       rixs_unit = cp_print_key_unit_nr(logger, rixs_section, "PRINT%SPECTRUM", &
     534              :                                        extension=".rixs", file_position="APPEND", &
     535           16 :                                        file_action="WRITE", file_form="FORMATTED")
     536              : 
     537           16 :       output_unit = cp_logger_get_default_io_unit()
     538              : 
     539           16 :       IF (rixs_unit > 0) THEN
     540              : 
     541              :          WRITE (rixs_unit, FMT="(A,/,T2,A,A,A,I3,A,A,A/,A)") &
     542            8 :             "=====================================================================================", &
     543            8 :             "Excitation from ground-state (", &
     544            8 :             rixs_env%core_state%state_type_char(donor_state%state_type), " for atom ", &
     545            8 :             donor_state%at_index, " of kind ", TRIM(donor_state%at_symbol), &
     546            8 :             ") to core-excited state i ", &
     547           16 :             "====================================================================================="
     548              : 
     549              :          WRITE (rixs_unit, FMT="(T3,A)") &
     550            8 :             "w_0i (eV)            mu^x_0i (a.u.)  mu^y_0i (a.u.)  mu^z_0i (a.u.)  mu^2_0i (a.u.)"
     551           91 :          DO i = 1, rixs_control%core_states
     552              :             WRITE (rixs_unit, FMT="(T2,F10.4,T26,E12.5,T42,E12.5,T58,E12.5,T74,E12.5)") &
     553           91 :                w_i0(i), mu_i0(1, i), mu_i0(2, i), mu_i0(3, i), mu_i0(4, i)
     554              :          END DO
     555              : 
     556              :          WRITE (rixs_unit, FMT="(A,/,T2,A,/,A)") &
     557            8 :             "=====================================================================================", &
     558            8 :             "Emission from core-excited state i to valence-excited state f ", &
     559           16 :             "====================================================================================="
     560              : 
     561              :          WRITE (rixs_unit, FMT="(T3,A)") &
     562            8 :             "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.)"
     563              : 
     564           91 :          DO i = 1, rixs_control%core_states
     565          363 :             DO f = 1, rixs_control%valence_states
     566              :                WRITE (rixs_unit, FMT="(T2,F10.4,T14,F8.4,T26,E12.5,T42,E12.5,T58,E12.5,T74,E12.5)") &
     567          355 :                   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)
     568              :             END DO
     569              :          END DO
     570              : 
     571              :       END IF
     572              : 
     573           16 :       CALL cp_print_key_finished_output(rixs_unit, logger, rixs_section, "PRINT%SPECTRUM")
     574              : 
     575           16 :    END SUBROUTINE print_rixs_to_file
     576              : 
     577              : END MODULE rixs_methods
        

Generated by: LCOV version 2.0-1