LCOV - code coverage report
Current view: top level - src - xas_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 708 779 90.9 %
Date: 2024-03-28 07:31:50 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief driver for the xas calculation and xas_scf for the tp method
      10             : !> \par History
      11             : !>      created 05.2005
      12             : !>      replace overlap integral routine [07.2014,JGH]
      13             : !> \author MI (05.2005)
      14             : ! **************************************************************************************************
      15             : MODULE xas_methods
      16             : 
      17             :    USE ai_contraction,                  ONLY: block_add,&
      18             :                                               contraction
      19             :    USE ai_overlap,                      ONLY: overlap_ab
      20             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      21             :                                               get_atomic_kind
      22             :    USE basis_set_types,                 ONLY: &
      23             :         allocate_sto_basis_set, create_gto_from_sto_basis, deallocate_sto_basis_set, &
      24             :         get_gto_basis_set, gto_basis_set_type, init_orb_basis_set, set_sto_basis_set, srules, &
      25             :         sto_basis_set_type
      26             :    USE cell_types,                      ONLY: cell_type,&
      27             :                                               pbc
      28             :    USE cp_array_utils,                  ONLY: cp_2d_r_p_type
      29             :    USE cp_control_types,                ONLY: dft_control_type
      30             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      31             :    USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr,&
      32             :                                               cp_dbcsr_sm_fm_multiply,&
      33             :                                               dbcsr_allocate_matrix_set
      34             :    USE cp_external_control,             ONLY: external_control
      35             :    USE cp_fm_pool_types,                ONLY: fm_pool_create_fm
      36             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      37             :                                               cp_fm_struct_release,&
      38             :                                               cp_fm_struct_type
      39             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      40             :                                               cp_fm_get_element,&
      41             :                                               cp_fm_get_submatrix,&
      42             :                                               cp_fm_release,&
      43             :                                               cp_fm_set_all,&
      44             :                                               cp_fm_set_submatrix,&
      45             :                                               cp_fm_to_fm,&
      46             :                                               cp_fm_type
      47             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      48             :                                               cp_logger_get_default_io_unit,&
      49             :                                               cp_logger_type,&
      50             :                                               cp_to_string
      51             :    USE cp_output_handling,              ONLY: cp_p_file,&
      52             :                                               cp_print_key_finished_output,&
      53             :                                               cp_print_key_should_output,&
      54             :                                               cp_print_key_unit_nr
      55             :    USE dbcsr_api,                       ONLY: dbcsr_convert_offsets_to_sizes,&
      56             :                                               dbcsr_copy,&
      57             :                                               dbcsr_create,&
      58             :                                               dbcsr_distribution_type,&
      59             :                                               dbcsr_p_type,&
      60             :                                               dbcsr_set,&
      61             :                                               dbcsr_type,&
      62             :                                               dbcsr_type_antisymmetric
      63             :    USE input_constants,                 ONLY: &
      64             :         do_loc_none, state_loc_list, state_loc_range, xas_1s_type, xas_2p_type, xas_2s_type, &
      65             :         xas_3d_type, xas_3p_type, xas_3s_type, xas_4d_type, xas_4f_type, xas_4p_type, xas_4s_type, &
      66             :         xas_dip_len, xas_dip_vel, xas_dscf, xas_tp_fh, xas_tp_flex, xas_tp_hh, xas_tp_xfh, &
      67             :         xas_tp_xhh, xes_tp_val
      68             :    USE input_section_types,             ONLY: section_get_lval,&
      69             :                                               section_vals_get_subs_vals,&
      70             :                                               section_vals_type,&
      71             :                                               section_vals_val_get
      72             :    USE kinds,                           ONLY: default_string_length,&
      73             :                                               dp
      74             :    USE memory_utilities,                ONLY: reallocate
      75             :    USE message_passing,                 ONLY: mp_para_env_type
      76             :    USE orbital_pointers,                ONLY: ncoset
      77             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      78             :    USE particle_methods,                ONLY: get_particle_set
      79             :    USE particle_types,                  ONLY: particle_type
      80             :    USE periodic_table,                  ONLY: ptable
      81             :    USE physcon,                         ONLY: evolt
      82             :    USE qs_diis,                         ONLY: qs_diis_b_clear,&
      83             :                                               qs_diis_b_create
      84             :    USE qs_environment_types,            ONLY: get_qs_env,&
      85             :                                               qs_environment_type,&
      86             :                                               set_qs_env
      87             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      88             :                                               qs_kind_type
      89             :    USE qs_loc_main,                     ONLY: qs_loc_driver
      90             :    USE qs_loc_methods,                  ONLY: qs_print_cubes
      91             :    USE qs_loc_types,                    ONLY: localized_wfn_control_type,&
      92             :                                               qs_loc_env_create,&
      93             :                                               qs_loc_env_type
      94             :    USE qs_loc_utils,                    ONLY: qs_loc_control_init,&
      95             :                                               qs_loc_env_init,&
      96             :                                               set_loc_centers,&
      97             :                                               set_loc_wfn_lists
      98             :    USE qs_matrix_pools,                 ONLY: mpools_get,&
      99             :                                               qs_matrix_pools_type
     100             :    USE qs_mo_io,                        ONLY: write_mo_set_to_output_unit
     101             :    USE qs_mo_methods,                   ONLY: calculate_subspace_eigenvalues
     102             :    USE qs_mo_types,                     ONLY: get_mo_set,&
     103             :                                               mo_set_type,&
     104             :                                               set_mo_set
     105             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
     106             :    USE qs_operators_ao,                 ONLY: p_xyz_ao,&
     107             :                                               rRc_xyz_ao
     108             :    USE qs_pdos,                         ONLY: calculate_projected_dos
     109             :    USE qs_scf,                          ONLY: scf_env_cleanup
     110             :    USE qs_scf_initialization,           ONLY: qs_scf_env_initialize
     111             :    USE qs_scf_types,                    ONLY: qs_scf_env_type,&
     112             :                                               scf_env_release
     113             :    USE scf_control_types,               ONLY: scf_c_create,&
     114             :                                               scf_c_read_parameters,&
     115             :                                               scf_control_type
     116             :    USE xas_control,                     ONLY: read_xas_control,&
     117             :                                               write_xas_control,&
     118             :                                               xas_control_create,&
     119             :                                               xas_control_type
     120             :    USE xas_env_types,                   ONLY: get_xas_env,&
     121             :                                               set_xas_env,&
     122             :                                               xas_env_create,&
     123             :                                               xas_env_release,&
     124             :                                               xas_environment_type
     125             :    USE xas_restart,                     ONLY: xas_read_restart
     126             :    USE xas_tp_scf,                      ONLY: xas_do_tp_scf,&
     127             :                                               xes_scf_once
     128             : #include "./base/base_uses.f90"
     129             : 
     130             :    IMPLICIT NONE
     131             :    PRIVATE
     132             : 
     133             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_methods'
     134             : 
     135             : ! *** Public subroutines ***
     136             : 
     137             :    PUBLIC :: xas, calc_stogto_overlap
     138             : 
     139             : CONTAINS
     140             : 
     141             : ! **************************************************************************************************
     142             : !> \brief Driver for xas calculations
     143             : !>      The initial mos are prepared
     144             : !>      A loop on the atoms to be excited is started
     145             : !>      For each atom the state to be excited is identified
     146             : !>      An scf optimization using the TP scheme or TD-DFT is used
     147             : !>      to evaluate the spectral energies and oscillator strengths
     148             : !> \param qs_env the qs_env, the xas_env lives in
     149             : !> \param dft_control ...
     150             : !> \par History
     151             : !>      05.2005 created [MI]
     152             : !> \author MI
     153             : !> \note
     154             : !>      the iteration counter is not finalized yet
     155             : !>      only the transition potential approach is active
     156             : !>      the localization can be switched off, otherwise
     157             : !>      it uses by default the berry phase approach
     158             : !>      The number of states to be localized is xas_control%nexc_search
     159             : !>      In general only the core states are needed
     160             : ! **************************************************************************************************
     161          42 :    SUBROUTINE xas(qs_env, dft_control)
     162             : 
     163             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     164             :       TYPE(dft_control_type), POINTER                    :: dft_control
     165             : 
     166             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'xas'
     167             : 
     168             :       INTEGER :: handle, homo, i, iat, iatom, ispin, istate, my_homo(2), my_nelectron(2), my_spin, &
     169             :          nao, nexc_atoms, nexc_search, nmo, nspins, output_unit, state_to_be_excited
     170             :       INTEGER, DIMENSION(2)                              :: added_mos
     171          42 :       INTEGER, DIMENSION(:), POINTER                     :: nexc_states
     172          42 :       INTEGER, DIMENSION(:, :), POINTER                  :: state_of_atom
     173             :       LOGICAL                                            :: ch_method_flags, converged, my_uocc(2), &
     174             :                                                             should_stop, skip_scf, &
     175             :                                                             transition_potential
     176             :       REAL(dp)                                           :: maxocc, occ_estate, tmp, xas_nelectron
     177          42 :       REAL(dp), DIMENSION(:), POINTER                    :: eigenvalues
     178          42 :       REAL(dp), DIMENSION(:, :), POINTER                 :: vecbuffer
     179          42 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     180             :       TYPE(cell_type), POINTER                           :: cell
     181          42 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: groundstate_coeff
     182             :       TYPE(cp_fm_type), POINTER                          :: all_vectors, excvec_coeff, mo_coeff
     183             :       TYPE(cp_logger_type), POINTER                      :: logger
     184          42 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, op_sm, ostrength_sm
     185             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff_b
     186          42 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     187             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     188          42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     189          42 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     190             :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
     191             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     192             :       TYPE(scf_control_type), POINTER                    :: scf_control
     193             :       TYPE(section_vals_type), POINTER                   :: dft_section, loc_section, &
     194             :                                                             print_loc_section, scf_section, &
     195             :                                                             xas_section
     196             :       TYPE(xas_control_type), POINTER                    :: xas_control
     197             :       TYPE(xas_environment_type), POINTER                :: xas_env
     198             : 
     199          42 :       CALL timeset(routineN, handle)
     200             : 
     201          42 :       transition_potential = .FALSE.
     202          42 :       skip_scf = .FALSE.
     203          42 :       converged = .TRUE.
     204          42 :       should_stop = .FALSE.
     205          42 :       ch_method_flags = .FALSE.
     206             : 
     207          42 :       NULLIFY (logger)
     208          42 :       logger => cp_get_default_logger()
     209          42 :       output_unit = cp_logger_get_default_io_unit(logger)
     210             : 
     211          42 :       NULLIFY (xas_env, groundstate_coeff, ostrength_sm, op_sm)
     212          42 :       NULLIFY (excvec_coeff, qs_loc_env, cell, scf_env)
     213          42 :       NULLIFY (matrix_ks)
     214          42 :       NULLIFY (all_vectors, state_of_atom, nexc_states, xas_control)
     215          42 :       NULLIFY (vecbuffer, op_sm, mo_coeff_b)
     216          42 :       NULLIFY (dft_section, xas_section, scf_section, loc_section, print_loc_section)
     217          42 :       dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
     218          42 :       xas_section => section_vals_get_subs_vals(dft_section, "XAS")
     219          42 :       scf_section => section_vals_get_subs_vals(xas_section, "SCF")
     220          42 :       loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE")
     221          42 :       print_loc_section => section_vals_get_subs_vals(loc_section, "PRINT")
     222             : 
     223             :       output_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%PROGRAM_RUN_INFO", &
     224          42 :                                          extension=".Log")
     225          42 :       IF (output_unit > 0) THEN
     226             :          WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T25,A,/,T3,A,/)") &
     227          21 :             REPEAT("=", 77), &
     228          21 :             "START CORE LEVEL SPECTROSCOPY CALCULATION", &
     229          42 :             REPEAT("=", 77)
     230             :       END IF
     231             : 
     232             : !   Create the xas environment
     233          42 :       CALL get_qs_env(qs_env, xas_env=xas_env)
     234          42 :       IF (.NOT. ASSOCIATED(xas_env)) THEN
     235          42 :          IF (output_unit > 0) THEN
     236             :             WRITE (UNIT=output_unit, FMT="(/,T5,A)") &
     237          21 :                "Create and initialize the xas environment"
     238             :          END IF
     239          42 :          ALLOCATE (xas_env)
     240          42 :          CALL xas_env_create(xas_env)
     241          42 :          CALL xas_env_init(xas_env, qs_env, dft_section, logger)
     242          42 :          xas_control => dft_control%xas_control
     243          42 :          CALL set_qs_env(qs_env, xas_env=xas_env)
     244             :       END IF
     245             : 
     246             : !   Initialize the type of calculation
     247          42 :       NULLIFY (atomic_kind_set, qs_kind_set, scf_control, mos, para_env, particle_set)
     248             :       CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
     249             :                       cell=cell, scf_control=scf_control, &
     250             :                       matrix_ks=matrix_ks, mos=mos, para_env=para_env, &
     251          42 :                       particle_set=particle_set)
     252             : 
     253             : !   The eigenstate of the KS Hamiltonian are nedeed
     254          42 :       NULLIFY (mo_coeff, eigenvalues)
     255          42 :       IF (scf_control%use_ot) THEN
     256           2 :          IF (output_unit > 0) THEN
     257             :             WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") &
     258           1 :                "Get eigenstates and eigenvalues from ground state MOs"
     259             :          END IF
     260           6 :          DO ispin = 1, dft_control%nspins
     261             :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
     262           4 :                             eigenvalues=eigenvalues, homo=homo)
     263             :             CALL calculate_subspace_eigenvalues(mo_coeff, &
     264             :                                                 matrix_ks(ispin)%matrix, eigenvalues, &
     265           6 :                                                 do_rotation=.TRUE.)
     266             :          END DO
     267             :       END IF
     268             : !   In xas SCF we need to use the same number of MOS as for GS
     269         126 :       added_mos = scf_control%added_mos
     270          42 :       NULLIFY (scf_control)
     271             : !   Consider to use get function for this
     272          42 :       CALL get_xas_env(xas_env, scf_control=scf_control)
     273         126 :       scf_control%added_mos = added_mos
     274             : 
     275             : !   Set initial occupation numbers, and store the original ones
     276          42 :       my_homo = 0
     277          42 :       my_nelectron = 0
     278         126 :       DO ispin = 1, dft_control%nspins
     279             :          CALL get_mo_set(mos(ispin), nelectron=my_nelectron(ispin), maxocc=maxocc, &
     280         126 :                          homo=my_homo(ispin), uniform_occupation=my_uocc(ispin))
     281             :       END DO
     282             : 
     283          42 :       nspins = dft_control%nspins
     284             : ! at the moment the only implemented method for XAS and XES calculations
     285          42 :       transition_potential = .TRUE. !(xas_control%xas_method==xas_tp_hh).OR.&
     286             :       !                           (xas_control%xas_method==xas_tp_fh).OR.&
     287             :       !                          (xas_control%xas_method==xas_tp_xhh).OR.&
     288             :       !                         (xas_control%xas_method==xas_tp_xfh).OR.&
     289             :       !                        (xas_control%xas_method==xas_dscf)
     290          42 :       IF (nspins == 1 .AND. transition_potential) THEN
     291           0 :          CPABORT("XAS with TP method requires LSD calculations")
     292             :       END IF
     293             : 
     294             :       CALL get_xas_env(xas_env=xas_env, &
     295             :                        all_vectors=all_vectors, &
     296             :                        groundstate_coeff=groundstate_coeff, excvec_coeff=excvec_coeff, &
     297             :                        nexc_atoms=nexc_atoms, &
     298          42 :                        spin_channel=my_spin)
     299             : 
     300             : !   Set of states among which there is the state to be excited
     301          42 :       CALL get_mo_set(mos(my_spin), nao=nao, homo=homo)
     302          42 :       IF (xas_control%nexc_search < 0) xas_control%nexc_search = homo
     303          42 :       nexc_search = xas_control%nexc_search
     304             : 
     305          42 :       CALL set_xas_env(xas_env=xas_env, nexc_search=nexc_search)
     306             : 
     307             :       !Define the qs_loc_env : to find centers, spread and possibly localize them
     308          42 :       CALL get_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env)
     309          42 :       IF (qs_loc_env%do_localize) THEN
     310          42 :          IF (output_unit > 0) THEN
     311             :             WRITE (UNIT=output_unit, FMT="(/,T2,A34,I3,A36/)") &
     312          21 :                "Localize a sub-set of MOs of spin ", my_spin, ","// &
     313          42 :                " to better identify the core states"
     314          21 :             IF ( &
     315             :                qs_loc_env%localized_wfn_control%set_of_states == state_loc_range) THEN
     316          19 :                WRITE (UNIT=output_unit, FMT="( A , I7, A, I7)") " The sub-set contains states from ", &
     317          19 :                   qs_loc_env%localized_wfn_control%lu_bound_states(1, my_spin), " to ", &
     318          38 :                   qs_loc_env%localized_wfn_control%lu_bound_states(2, my_spin)
     319           2 :             ELSEIF (qs_loc_env%localized_wfn_control%set_of_states == state_loc_list) THEN
     320           2 :                WRITE (UNIT=output_unit, FMT="( A )") " The sub-set contains states given in the input list"
     321             :             END IF
     322             : 
     323             :          END IF
     324          42 :          CALL qs_loc_driver(qs_env, qs_loc_env, print_loc_section, myspin=my_spin)
     325             :       END IF
     326             : 
     327          42 :       CPASSERT(ASSOCIATED(groundstate_coeff))
     328         126 :       DO ispin = 1, nspins
     329          84 :          CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, nmo=nmo)
     330          84 :          CALL cp_fm_to_fm(mo_coeff, groundstate_coeff(ispin), nmo, 1, 1)
     331         126 :          IF (ASSOCIATED(mo_coeff_b)) THEN
     332             : 
     333             :          END IF
     334             :       END DO
     335             : 
     336             : !   SCF for only XES using occupied core and empty homo (only one SCF)
     337             : !   Probably better not to do the localization in this case, but only single out the
     338             : !   core orbital for the specific atom for which the spectrum is computed
     339          42 :       IF (xas_control%xas_method == xes_tp_val .AND. &
     340             :           xas_control%xes_core_occupation == 1.0_dp) THEN
     341           4 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,/,T10,A)') &
     342           2 :             "START Core Level Spectroscopy Calculation for the Emission Spectrum"
     343           4 :          IF (xas_control%xes_homo_occupation == 1) THEN
     344           2 :             IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(T10,A,/,A)') &
     345           1 :                "The core state is fully occupied and XES from ground state calculation.", &
     346           2 :                " No SCF is needed, MOS already available"
     347           2 :          ELSE IF (xas_control%xes_homo_occupation == 0) THEN
     348           2 :             IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(T10,A,/,A)') &
     349           1 :                "The core state is fully occupied and the homo is empty", &
     350           2 :                " (final state of the core hole decay). Only one SCF is needed (not one per atom)"
     351             :          END IF
     352           4 :          skip_scf = .TRUE.
     353             : 
     354           4 :          CALL set_xas_env(xas_env=xas_env, xas_estate=-1, homo_occ=xas_control%xes_homo_occupation)
     355           4 :          CALL xes_scf_once(qs_env, xas_env, converged, should_stop)
     356             : 
     357           4 :          IF (converged .AND. .NOT. should_stop .AND. xas_control%xes_homo_occupation == 0) THEN
     358           2 :             IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,T10,A,I6)') &
     359           1 :                "SCF with empty homo converged "
     360           2 :          ELSE IF (.NOT. converged .OR. should_stop) THEN
     361           0 :             IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,T10,A,I6)') &
     362           0 :                "SCF with empty homo NOT converged"
     363             :             ! Release what has to be released
     364             :             IF (ASSOCIATED(vecbuffer)) THEN
     365             :                DEALLOCATE (vecbuffer)
     366             :                DEALLOCATE (op_sm)
     367             :             END IF
     368             : 
     369           0 :             DO ispin = 1, dft_control%nspins
     370             :                CALL set_mo_set(mos(ispin), homo=my_homo(ispin), &
     371           0 :                                uniform_occupation=my_uocc(ispin), nelectron=my_nelectron(ispin))
     372           0 :                CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
     373           0 :                CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
     374             :             END DO
     375             : 
     376           0 :             IF (output_unit > 0) THEN
     377             :                WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T25,A,/,T3,A,/)") &
     378           0 :                   REPEAT("=", 77), &
     379           0 :                   "END CORE LEVEL SPECTROSCOPY CALCULATION", &
     380           0 :                   REPEAT("=", 77)
     381             :             END IF
     382             : 
     383           0 :             CALL xas_env_release(qs_env%xas_env)
     384           0 :             DEALLOCATE (qs_env%xas_env)
     385           0 :             NULLIFY (qs_env%xas_env)
     386             : 
     387             :             CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
     388           0 :                                               "PRINT%PROGRAM_RUN_INFO")
     389           0 :             CALL timestop(handle)
     390           0 :             RETURN
     391             :          END IF
     392             :       END IF
     393             : 
     394             :       ! Assign the character of the selected core states
     395             :       ! through the overlap with atomic-like states
     396             :       CALL cls_assign_core_states(xas_control, xas_env, qs_loc_env%localized_wfn_control, &
     397          42 :                                   qs_env)
     398             :       CALL get_xas_env(xas_env=xas_env, &
     399          42 :                        state_of_atom=state_of_atom, nexc_states=nexc_states)
     400             : 
     401          42 :       IF (skip_scf) THEN
     402           4 :          CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff)
     403             :          CALL cp_fm_to_fm(mo_coeff, all_vectors, ncol=nexc_search, &
     404           4 :                           source_start=1, target_start=1)
     405             :       END IF
     406             : 
     407         126 :       ALLOCATE (vecbuffer(1, nao))
     408         168 :       ALLOCATE (op_sm(3))
     409             : 
     410             :       ! copy the coefficients of the mos in a temporary fm with the right structure
     411             :       IF (transition_potential) THEN
     412             :          ! Calculate the operator
     413          42 :          CALL get_xas_env(xas_env=xas_env, ostrength_sm=ostrength_sm)
     414         168 :          DO i = 1, 3
     415         126 :             NULLIFY (op_sm(i)%matrix)
     416         168 :             op_sm(i)%matrix => ostrength_sm(i)%matrix
     417             :          END DO
     418          42 :          IF (xas_control%dipole_form == xas_dip_vel) THEN
     419          42 :             CALL p_xyz_ao(op_sm, qs_env)
     420             :          END IF
     421             :       END IF
     422             : 
     423             :       ! DO SCF if required
     424         124 :       DO iat = 1, nexc_atoms
     425          82 :          iatom = xas_env%exc_atoms(iat)
     426         212 :          DO istate = 1, nexc_states(iat)
     427             :             ! determine which state has to be excited in the global list
     428          88 :             state_to_be_excited = state_of_atom(iat, istate)
     429             : 
     430             :             ! Take the state_to_be_excited vector from the full set and copy into excvec_coeff
     431          88 :             CALL get_mo_set(mos(my_spin), nmo=nmo)
     432          88 :             CALL get_xas_env(xas_env, occ_estate=occ_estate, xas_nelectron=xas_nelectron)
     433          88 :             tmp = xas_nelectron + 1.0_dp - occ_estate
     434          88 :             IF (nmo < tmp) &
     435           0 :                CPABORT("CLS: the required method needs added_mos to the ground state")
     436             :             ! If the restart file for this atom exists, the mos and the
     437             :             ! occupation numbers are overwritten
     438             :             ! It is necessary that the restart is for the same xas method
     439             :             ! otherwise the number of electrons and the occupation numbers
     440             :             ! may not  be consistent
     441          88 :             IF (xas_control%xas_restart) THEN
     442             :                CALL xas_read_restart(xas_env, xas_section, qs_env, xas_control%xas_method, iatom, &
     443          12 :                                      state_to_be_excited, istate)
     444             :             END IF
     445          88 :             CALL set_xas_env(xas_env=xas_env, xas_estate=state_to_be_excited)
     446          88 :             CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff)
     447          88 :             CPASSERT(ASSOCIATED(excvec_coeff))
     448             :             CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, state_to_be_excited, &
     449          88 :                                      nao, 1, transpose=.TRUE.)
     450             :             CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, &
     451          88 :                                      nao, 1, transpose=.TRUE.)
     452             : 
     453             :             IF (transition_potential) THEN
     454             : 
     455          88 :                IF (.NOT. skip_scf) THEN
     456          80 :                   IF (output_unit > 0) THEN
     457          40 :                      WRITE (UNIT=output_unit, FMT='(/,T5,A)') REPEAT("-", 75)
     458          40 :                      IF (xas_control%xas_method == xas_dscf) THEN
     459             :                         WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') &
     460           0 :                            "START DeltaSCF for the first excited state from the core state of ATOM ", iatom
     461             :                      ELSE
     462             :                         WRITE (UNIT=output_unit, FMT='(/,T10,A,I6)') &
     463          40 :                            "Start Core Level Spectroscopy Calculation with TP approach for ATOM ", iatom
     464             :                         WRITE (UNIT=output_unit, FMT='(/,T10,A,I6,T34,A,T54,I6)') &
     465          40 :                            "Excited state", istate, "out of", nexc_states(iat)
     466          40 :                         WRITE (UNIT=output_unit, FMT='(T10,A,T50,f10.4)') "Occupation of the core orbital", &
     467          80 :                            occ_estate
     468          40 :                         WRITE (UNIT=output_unit, FMT='(T10,A28,I3, T50,F10.4)') "Number of electrons in Spin ", &
     469          80 :                            my_spin, xas_nelectron
     470             :                      END IF
     471             :                   END IF
     472             : 
     473          80 :                   CALL get_xas_env(xas_env=xas_env, scf_env=scf_env)
     474          80 :                   IF (.NOT. ASSOCIATED(scf_env)) THEN
     475          44 :                      CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
     476             :                      !     Moved here from qs_scf_env_initialize to be able to have more scf_env
     477          44 :                      CALL set_xas_env(xas_env, scf_env=scf_env)
     478             :                   ELSE
     479          36 :                      CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
     480             :                   END IF
     481             : 
     482         240 :                   DO ispin = 1, SIZE(mos)
     483         240 :                      IF (ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN !fm->dbcsr
     484             :                         CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
     485         160 :                                               mos(ispin)%mo_coeff_b) !fm->dbcsr
     486             :                      END IF !fm->dbcsr
     487             :                   END DO !fm->dbcsr
     488             : 
     489          80 :                   IF (.NOT. scf_env%skip_diis) THEN
     490          68 :                      IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
     491          40 :                         ALLOCATE (scf_env%scf_diis_buffer)
     492          40 :                         CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
     493             :                      END IF
     494          68 :                      CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
     495             :                   END IF
     496             : 
     497             :                   CALL xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, &
     498          80 :                                      xas_section, scf_section, converged, should_stop)
     499             : 
     500             :                   CALL external_control(should_stop, "CLS", target_time=qs_env%target_time, &
     501          80 :                                         start_time=qs_env%start_time)
     502          80 :                   IF (should_stop) THEN
     503           0 :                      CALL scf_env_cleanup(scf_env)
     504           0 :                      EXIT
     505             :                   END IF
     506             : 
     507             :                END IF
     508             :                ! SCF DONE
     509             : 
     510             :                ! Write last wavefunction to screen
     511          88 :                IF (SIZE(mos) > 1) THEN
     512             :                   CALL write_mo_set_to_output_unit(mos(1), atomic_kind_set, qs_kind_set, particle_set, &
     513          88 :                                                    dft_section, 4, 0, final_mos=.FALSE., spin="XAS ALPHA")
     514             :                   CALL write_mo_set_to_output_unit(mos(2), atomic_kind_set, qs_kind_set, particle_set, &
     515          88 :                                                    dft_section, 4, 0, final_mos=.FALSE., spin="XAS BETA")
     516             :                ELSE
     517             :                   CALL write_mo_set_to_output_unit(mos(1), atomic_kind_set, qs_kind_set, particle_set, &
     518           0 :                                                    dft_section, 4, 0, final_mos=.FALSE., spin="XAS")
     519             :                END IF
     520             : 
     521             :             ELSE
     522             :                ! Core level spectroscopy by TDDFT is not yet implemented
     523             :                ! the states defined by the rotation are the ground state orbitals
     524             :                ! the initial state from which I excite should be localized
     525             :                ! I take the excitations from lumo to nmo
     526             :             END IF
     527             : 
     528          88 :             IF (converged) THEN
     529             :                CALL cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, &
     530          54 :                                            iatom, istate)
     531             :             ELSE
     532          34 :                IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') &
     533          17 :                   "SCF with core hole NOT converged for ATOM ", iatom
     534             :             END IF
     535             : 
     536         258 :             IF (.NOT. skip_scf) THEN
     537             :                ! Reset the initial core orbitals.
     538             :                ! The valence orbitals are taken from the last SCF,
     539             :                ! it should be a better initial guess
     540          80 :                CALL get_qs_env(qs_env, mos=mos)
     541         240 :                DO ispin = 1, nspins
     542         160 :                   CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
     543         240 :                   CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
     544             :                END DO
     545          80 :                IF (iat == nexc_atoms) THEN
     546          44 :                   CALL scf_env_cleanup(scf_env)
     547          44 :                   CALL scf_env_release(xas_env%scf_env)
     548          44 :                   DEALLOCATE (xas_env%scf_env)
     549             :                END IF
     550             :             END IF
     551             : 
     552             :          END DO ! istate
     553             :       END DO ! iat = 1,nexc_atoms
     554             : 
     555             :       ! END of Calculation
     556             : 
     557             :       ! Release what has to be released
     558          42 :       IF (ASSOCIATED(vecbuffer)) THEN
     559          42 :          DEALLOCATE (vecbuffer)
     560          42 :          DEALLOCATE (op_sm)
     561             :       END IF
     562             : 
     563         126 :       DO ispin = 1, dft_control%nspins
     564             :          CALL set_mo_set(mos(ispin), homo=my_homo(ispin), &
     565          84 :                          uniform_occupation=my_uocc(ispin), nelectron=my_nelectron(ispin))
     566          84 :          CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
     567         126 :          CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
     568             :       END DO
     569             : 
     570          42 :       IF (output_unit > 0) THEN
     571             :          WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T25,A,/,T3,A,/)") &
     572          21 :             REPEAT("=", 77), &
     573          21 :             "END CORE LEVEL SPECTROSCOPY CALCULATION", &
     574          42 :             REPEAT("=", 77)
     575             :       END IF
     576             : 
     577          42 :       CALL xas_env_release(qs_env%xas_env)
     578          42 :       DEALLOCATE (qs_env%xas_env)
     579          42 :       NULLIFY (qs_env%xas_env)
     580             : 
     581             :       CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
     582          42 :                                         "PRINT%PROGRAM_RUN_INFO")
     583          42 :       CALL timestop(handle)
     584             : 
     585          84 :    END SUBROUTINE xas
     586             : 
     587             : ! **************************************************************************************************
     588             : !> \brief allocate and initialize the structure needed for the xas calculation
     589             : !> \param xas_env the environment for XAS  calculations
     590             : !> \param qs_env the qs_env, the xas_env lives in
     591             : !> \param dft_section ...
     592             : !> \param logger ...
     593             : !> \par History
     594             : !>      05.2005 created [MI]
     595             : !> \author MI
     596             : ! **************************************************************************************************
     597          42 :    SUBROUTINE xas_env_init(xas_env, qs_env, dft_section, logger)
     598             : 
     599             :       TYPE(xas_environment_type), POINTER                :: xas_env
     600             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     601             :       TYPE(section_vals_type), POINTER                   :: dft_section
     602             :       TYPE(cp_logger_type), POINTER                      :: logger
     603             : 
     604             :       CHARACTER(LEN=default_string_length)               :: name_sto
     605             :       INTEGER :: homo, i, iat, iatom, ik, ikind, ispin, j, l, lfomo, my_spin, n_mo(2), n_rep, nao, &
     606             :          natom, ncubes, nelectron, nexc_atoms, nexc_search, nj, nk, nkind, nmo, nmoloc(2), &
     607             :          nsgf_gto, nsgf_sto, nspins, nvirtual, nvirtual2
     608          42 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: first_sgf, kind_type_tmp, kind_z_tmp, &
     609          42 :                                                             last_sgf
     610             :       INTEGER, DIMENSION(4, 7)                           :: ne
     611          42 :       INTEGER, DIMENSION(:), POINTER                     :: bounds, list, lq, nq, row_blk_sizes
     612             :       LOGICAL                                            :: ihavethis
     613             :       REAL(dp)                                           :: nele, occ_estate, occ_homo, &
     614             :                                                             occ_homo_plus, zatom
     615          42 :       REAL(dp), DIMENSION(:), POINTER                    :: sto_zet
     616          42 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     617             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     618             :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
     619             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     620             :       TYPE(dbcsr_distribution_type), POINTER             :: dbcsr_dist
     621          42 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     622             :       TYPE(dft_control_type), POINTER                    :: dft_control
     623             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     624          42 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     625             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     626             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     627          42 :          POINTER                                         :: sab_orb
     628          42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     629          42 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     630             :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
     631             :       TYPE(qs_matrix_pools_type), POINTER                :: mpools
     632             :       TYPE(scf_control_type), POINTER                    :: scf_control
     633             :       TYPE(section_vals_type), POINTER                   :: loc_section, xas_section
     634             :       TYPE(sto_basis_set_type), POINTER                  :: sto_basis_set
     635             :       TYPE(xas_control_type), POINTER                    :: xas_control
     636             : 
     637          42 :       n_mo(1:2) = 0
     638           0 :       CPASSERT(ASSOCIATED(xas_env))
     639             : 
     640          42 :       NULLIFY (atomic_kind_set, qs_kind_set, dft_control, scf_control, matrix_s, mos, mpools)
     641          42 :       NULLIFY (para_env, particle_set, xas_control)
     642          42 :       NULLIFY (qs_loc_env)
     643          42 :       NULLIFY (sab_orb)
     644             :       CALL get_qs_env(qs_env=qs_env, &
     645             :                       atomic_kind_set=atomic_kind_set, &
     646             :                       qs_kind_set=qs_kind_set, &
     647             :                       dft_control=dft_control, &
     648             :                       mpools=mpools, &
     649             :                       matrix_s=matrix_s, mos=mos, &
     650             :                       para_env=para_env, particle_set=particle_set, &
     651             :                       sab_orb=sab_orb, &
     652          42 :                       dbcsr_dist=dbcsr_dist)
     653             : 
     654          42 :       xas_section => section_vals_get_subs_vals(dft_section, "XAS")
     655          42 :       ALLOCATE (dft_control%xas_control)
     656          42 :       CALL xas_control_create(dft_control%xas_control)
     657          42 :       CALL read_xas_control(dft_control%xas_control, xas_section)
     658          42 :       CALL write_xas_control(dft_control%xas_control, dft_section)
     659          42 :       xas_control => dft_control%xas_control
     660          42 :       ALLOCATE (scf_control)
     661          42 :       CALL scf_c_create(scf_control)
     662          42 :       CALL scf_c_read_parameters(scf_control, xas_section)
     663          42 :       CALL set_xas_env(xas_env, scf_control=scf_control)
     664             : 
     665          42 :       my_spin = xas_control%spin_channel
     666          42 :       nexc_search = xas_control%nexc_search
     667          42 :       IF (nexc_search < 0) THEN
     668             :          ! ground state occupation
     669           2 :          CALL get_mo_set(mos(my_spin), nmo=nmo, lfomo=lfomo)
     670           2 :          nexc_search = lfomo - 1
     671             :       END IF
     672          42 :       nexc_atoms = xas_control%nexc_atoms
     673         126 :       ALLOCATE (xas_env%exc_atoms(nexc_atoms))
     674         206 :       xas_env%exc_atoms = xas_control%exc_atoms
     675             :       CALL set_xas_env(xas_env=xas_env, nexc_search=nexc_search, &
     676          42 :                        nexc_atoms=nexc_atoms, spin_channel=my_spin)
     677             : 
     678          42 :       CALL mpools_get(mpools, ao_mo_fm_pools=xas_env%ao_mo_fm_pools)
     679             : 
     680          42 :       NULLIFY (mo_coeff)
     681          42 :       CALL get_mo_set(mos(my_spin), nao=nao, homo=homo, nmo=nmo, mo_coeff=mo_coeff, nelectron=nelectron)
     682             : 
     683          42 :       nvirtual2 = 0
     684          42 :       IF (xas_control%added_mos .GT. 0) THEN
     685          40 :          nvirtual2 = MIN(xas_control%added_mos, nao - nmo)
     686          40 :          xas_env%unoccupied_eps = xas_control%eps_added
     687          40 :          xas_env%unoccupied_max_iter = xas_control%max_iter_added
     688             :       END IF
     689          42 :       nvirtual = nmo + nvirtual2
     690             : 
     691         126 :       n_mo(1:2) = nmo
     692             : 
     693         126 :       ALLOCATE (xas_env%centers_wfn(3, nexc_search))
     694         126 :       ALLOCATE (xas_env%atom_of_state(nexc_search))
     695         126 :       ALLOCATE (xas_env%type_of_state(nexc_search))
     696         168 :       ALLOCATE (xas_env%state_of_atom(nexc_atoms, nexc_search))
     697         126 :       ALLOCATE (xas_env%nexc_states(nexc_atoms))
     698         126 :       ALLOCATE (xas_env%mykind_of_atom(nexc_atoms))
     699          42 :       nkind = SIZE(atomic_kind_set, 1)
     700         126 :       ALLOCATE (xas_env%mykind_of_kind(nkind))
     701         124 :       xas_env%mykind_of_kind = 0
     702             : 
     703             :       ! create a new matrix structure nao x 1
     704          42 :       NULLIFY (tmp_fm_struct)
     705             :       CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
     706          42 :                                ncol_global=1, para_env=para_env, context=mo_coeff%matrix_struct%context)
     707          42 :       ALLOCATE (xas_env%excvec_coeff)
     708          42 :       CALL cp_fm_create(xas_env%excvec_coeff, tmp_fm_struct)
     709          42 :       CALL cp_fm_struct_release(tmp_fm_struct)
     710             : 
     711          42 :       NULLIFY (tmp_fm_struct)
     712             :       CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=1, &
     713             :                                ncol_global=nexc_search, para_env=para_env, &
     714          42 :                                context=mo_coeff%matrix_struct%context)
     715          42 :       ALLOCATE (xas_env%excvec_overlap)
     716          42 :       CALL cp_fm_create(xas_env%excvec_overlap, tmp_fm_struct)
     717          42 :       CALL cp_fm_struct_release(tmp_fm_struct)
     718             : 
     719          42 :       nspins = SIZE(mos, 1)
     720             : 
     721             :       ! initialize operators for the calculation of the oscillator strengths
     722          42 :       IF (xas_control%xas_method == xas_tp_hh) THEN
     723           8 :          occ_estate = 0.5_dp
     724           8 :          nele = REAL(nelectron, dp) - 0.5_dp
     725           8 :          occ_homo = 1.0_dp
     726           8 :          occ_homo_plus = 0._dp
     727             :       ELSEIF (xas_control%xas_method == xas_tp_xhh) THEN
     728           4 :          occ_estate = 0.5_dp
     729           4 :          nele = REAL(nelectron, dp)
     730           4 :          occ_homo = 1.0_dp
     731           4 :          occ_homo_plus = 0.5_dp
     732             :       ELSEIF (xas_control%xas_method == xas_tp_fh) THEN
     733          10 :          occ_estate = 0.0_dp
     734          10 :          nele = REAL(nelectron, dp) - 1.0_dp
     735          10 :          occ_homo = 1.0_dp
     736          10 :          occ_homo_plus = 0._dp
     737             :       ELSEIF (xas_control%xas_method == xas_tp_xfh) THEN
     738           8 :          occ_estate = 0.0_dp
     739           8 :          nele = REAL(nelectron, dp)
     740           8 :          occ_homo = 1.0_dp
     741           8 :          occ_homo_plus = 1._dp
     742             :       ELSEIF (xas_control%xas_method == xes_tp_val) THEN
     743           6 :          occ_estate = xas_control%xes_core_occupation
     744           6 :          nele = REAL(nelectron, dp) - xas_control%xes_core_occupation
     745           6 :          occ_homo = xas_control%xes_homo_occupation
     746             :       ELSEIF (xas_control%xas_method == xas_dscf) THEN
     747           0 :          occ_estate = 0.0_dp
     748           0 :          nele = REAL(nelectron, dp)
     749           0 :          occ_homo = 1.0_dp
     750           0 :          occ_homo_plus = 1._dp
     751             :       ELSEIF (xas_control%xas_method == xas_tp_flex) THEN
     752           6 :          nele = REAL(xas_control%nel_tot, dp)
     753           6 :          occ_estate = REAL(xas_control%xas_core_occupation, dp)
     754           6 :          IF (nele < 0.0_dp) nele = REAL(nelectron, dp) - (1.0_dp - occ_estate)
     755           6 :          occ_homo = 1.0_dp
     756             :       END IF
     757             :       CALL set_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_nelectron=nele, &
     758          42 :                        nvirtual2=nvirtual2, nvirtual=nvirtual, homo_occ=occ_homo)
     759             : 
     760             :       ! Initialize the list of orbitals for cube files printing
     761          42 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, xas_section, &
     762             :                                            "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
     763           2 :          NULLIFY (bounds, list)
     764             :          CALL section_vals_val_get(xas_section, &
     765             :                                    "PRINT%CLS_FUNCTION_CUBES%CUBES_LU_BOUNDS", &
     766           2 :                                    i_vals=bounds)
     767           2 :          ncubes = bounds(2) - bounds(1) + 1
     768           2 :          IF (ncubes > 0) THEN
     769           0 :             ALLOCATE (xas_control%list_cubes(ncubes))
     770             : 
     771           0 :             DO ik = 1, ncubes
     772           0 :                xas_control%list_cubes(ik) = bounds(1) + (ik - 1)
     773             :             END DO
     774             :          END IF
     775             : 
     776           2 :          IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
     777             :             CALL section_vals_val_get(xas_section, &
     778             :                                       "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
     779           2 :                                       n_rep_val=n_rep)
     780           2 :             ncubes = 0
     781           4 :             DO ik = 1, n_rep
     782           2 :                NULLIFY (list)
     783             :                CALL section_vals_val_get(xas_section, &
     784             :                                          "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
     785           2 :                                          i_rep_val=ik, i_vals=list)
     786           4 :                IF (ASSOCIATED(list)) THEN
     787           2 :                   CALL reallocate(xas_control%list_cubes, 1, ncubes + SIZE(list))
     788           8 :                   DO i = 1, SIZE(list)
     789           8 :                      xas_control%list_cubes(i + ncubes) = list(i)
     790             :                   END DO
     791           2 :                   ncubes = ncubes + SIZE(list)
     792             :                END IF
     793             :             END DO ! ik
     794             :          END IF
     795             : 
     796           2 :          IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
     797           0 :             ncubes = MAX(10, xas_control%added_mos/10)
     798           0 :             ncubes = MIN(ncubes, xas_control%added_mos)
     799           0 :             ALLOCATE (xas_control%list_cubes(ncubes))
     800           0 :             DO ik = 1, ncubes
     801           0 :                xas_control%list_cubes(ik) = homo + ik
     802             :             END DO
     803             :          END IF
     804             :       ELSE
     805          40 :          NULLIFY (xas_control%list_cubes)
     806             :       END IF
     807             : 
     808          42 :       NULLIFY (tmp_fm_struct)
     809         210 :       ALLOCATE (xas_env%groundstate_coeff(nspins))
     810         126 :       DO ispin = 1, nspins
     811          84 :          CALL get_mo_set(mos(ispin), nao=nao, nmo=nmo)
     812             :          CALL fm_pool_create_fm(xas_env%ao_mo_fm_pools(ispin)%pool, &
     813             :                                 xas_env%groundstate_coeff(ispin), &
     814         126 :                                 name="xas_env%mo0"//TRIM(ADJUSTL(cp_to_string(ispin))))
     815             :       END DO ! ispin
     816             : 
     817          42 :       NULLIFY (tmp_fm_struct)
     818             :       CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=1, &
     819             :                                ncol_global=nvirtual, para_env=para_env, &
     820          42 :                                context=mo_coeff%matrix_struct%context)
     821         420 :       ALLOCATE (xas_env%dip_fm_set(2, 3))
     822         168 :       DO i = 1, 3
     823         420 :          DO j = 1, 2
     824         378 :             CALL cp_fm_create(xas_env%dip_fm_set(j, i), tmp_fm_struct)
     825             :          END DO
     826             :       END DO
     827          42 :       CALL cp_fm_struct_release(tmp_fm_struct)
     828             : 
     829             :       !Array to store all the eigenstates: occupied and the required not occupied
     830          42 :       IF (nvirtual2 .GT. 0) THEN
     831         120 :          ALLOCATE (xas_env%unoccupied_evals(nvirtual2))
     832          40 :          NULLIFY (tmp_fm_struct)
     833             :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
     834             :                                   ncol_global=nvirtual2, &
     835          40 :                                   para_env=para_env, context=mo_coeff%matrix_struct%context)
     836          40 :          ALLOCATE (xas_env%unoccupied_orbs)
     837          40 :          CALL cp_fm_create(xas_env%unoccupied_orbs, tmp_fm_struct)
     838          40 :          CALL cp_fm_struct_release(tmp_fm_struct)
     839             :       END IF
     840             : 
     841          42 :       NULLIFY (tmp_fm_struct)
     842             :       CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
     843             :                                ncol_global=nvirtual, &
     844          42 :                                para_env=para_env, context=mo_coeff%matrix_struct%context)
     845          42 :       ALLOCATE (xas_env%all_vectors)
     846          42 :       CALL cp_fm_create(xas_env%all_vectors, tmp_fm_struct)
     847          42 :       CALL cp_fm_struct_release(tmp_fm_struct)
     848             : 
     849             :       ! Array to store all the energies needed  for the spectrum
     850         126 :       ALLOCATE (xas_env%all_evals(nvirtual))
     851             : 
     852          42 :       IF (xas_control%dipole_form == xas_dip_len) THEN
     853           0 :          CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
     854           0 :          DO i = 1, 3
     855           0 :             ALLOCATE (xas_env%ostrength_sm(i)%matrix)
     856             :             CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, matrix_s(1)%matrix, &
     857           0 :                             "xas_env%ostrength_sm-"//TRIM(ADJUSTL(cp_to_string(i))))
     858           0 :             CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
     859             :          END DO
     860          42 :       ELSEIF (xas_control%dipole_form == xas_dip_vel) THEN
     861             :          !
     862             :          ! prepare for allocation
     863          42 :          natom = SIZE(particle_set, 1)
     864         126 :          ALLOCATE (first_sgf(natom))
     865         126 :          ALLOCATE (last_sgf(natom))
     866             :          CALL get_particle_set(particle_set, qs_kind_set, &
     867             :                                first_sgf=first_sgf, &
     868          42 :                                last_sgf=last_sgf)
     869         126 :          ALLOCATE (row_blk_sizes(natom))
     870          42 :          CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
     871          42 :          DEALLOCATE (first_sgf)
     872          42 :          DEALLOCATE (last_sgf)
     873             :          !
     874             :          !
     875          42 :          CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
     876          42 :          ALLOCATE (xas_env%ostrength_sm(1)%matrix)
     877             :          CALL dbcsr_create(matrix=xas_env%ostrength_sm(1)%matrix, &
     878             :                            name="xas_env%ostrength_sm", &
     879             :                            dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
     880             :                            row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
     881          42 :                            nze=0, mutable_work=.TRUE.)
     882          42 :          CALL cp_dbcsr_alloc_block_from_nbl(xas_env%ostrength_sm(1)%matrix, sab_orb)
     883          42 :          CALL dbcsr_set(xas_env%ostrength_sm(1)%matrix, 0.0_dp)
     884         126 :          DO i = 2, 3
     885          84 :             ALLOCATE (xas_env%ostrength_sm(i)%matrix)
     886             :             CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, xas_env%ostrength_sm(1)%matrix, &
     887          84 :                             "xas_env%ostrength_sm-"//TRIM(ADJUSTL(cp_to_string(i))))
     888         126 :             CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
     889             :          END DO
     890             : 
     891          42 :          DEALLOCATE (row_blk_sizes)
     892             :       END IF
     893             : 
     894             :       ! Define the qs_loc_env : to find centers, spread and possibly localize them
     895          42 :       IF (.NOT. (ASSOCIATED(xas_env%qs_loc_env))) THEN
     896          42 :          ALLOCATE (qs_loc_env)
     897          42 :          CALL qs_loc_env_create(qs_loc_env)
     898          42 :          CALL set_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env)
     899          42 :          loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE")
     900             : 
     901             :          CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.TRUE., &
     902          42 :                                   do_xas=.TRUE., nloc_xas=nexc_search, spin_xas=my_spin)
     903             : 
     904          42 :          IF (.NOT. qs_loc_env%do_localize) THEN
     905           0 :             qs_loc_env%localized_wfn_control%localization_method = do_loc_none
     906             : 
     907             :          ELSE
     908         126 :             nmoloc = qs_loc_env%localized_wfn_control%nloc_states
     909          42 :             CALL set_loc_wfn_lists(qs_loc_env%localized_wfn_control, nmoloc, n_mo, nspins, my_spin)
     910          42 :             CALL set_loc_centers(qs_loc_env%localized_wfn_control, nmoloc, nspins)
     911             :             CALL qs_loc_env_init(qs_loc_env, qs_loc_env%localized_wfn_control, &
     912          42 :                                  qs_env, myspin=my_spin, do_localize=qs_loc_env%do_localize)
     913             :          END IF
     914             :       END IF
     915             : 
     916             :       !Type of state
     917          42 :       ALLOCATE (nq(1), lq(1), sto_zet(1))
     918          42 :       IF (xas_control%state_type == xas_1s_type) THEN
     919          40 :          nq(1) = 1
     920          40 :          lq(1) = 0
     921             :       ELSEIF (xas_control%state_type == xas_2s_type) THEN
     922           0 :          nq(1) = 2
     923           0 :          lq(1) = 0
     924             :       ELSEIF (xas_control%state_type == xas_2p_type) THEN
     925           2 :          nq(1) = 2
     926           2 :          lq(1) = 1
     927             :       ELSEIF (xas_control%state_type == xas_3s_type) THEN
     928           0 :          nq(1) = 3
     929           0 :          lq(1) = 0
     930             :       ELSEIF (xas_control%state_type == xas_3p_type) THEN
     931           0 :          nq(1) = 3
     932           0 :          lq(1) = 1
     933             :       ELSEIF (xas_control%state_type == xas_3d_type) THEN
     934           0 :          nq(1) = 3
     935           0 :          lq(1) = 2
     936             :       ELSEIF (xas_control%state_type == xas_4s_type) THEN
     937           0 :          nq(1) = 4
     938           0 :          lq(1) = 0
     939             :       ELSEIF (xas_control%state_type == xas_4p_type) THEN
     940           0 :          nq(1) = 4
     941           0 :          lq(1) = 1
     942             :       ELSEIF (xas_control%state_type == xas_4d_type) THEN
     943           0 :          nq(1) = 4
     944           0 :          lq(1) = 2
     945             :       ELSEIF (xas_control%state_type == xas_4f_type) THEN
     946           0 :          nq(1) = 4
     947           0 :          lq(1) = 3
     948             :       ELSE
     949           0 :          CPABORT("XAS type of state not implemented")
     950             :       END IF
     951             : 
     952             : !   Find core orbitals of right angular momentum
     953         126 :       ALLOCATE (kind_type_tmp(nkind))
     954         126 :       ALLOCATE (kind_z_tmp(nkind))
     955         124 :       kind_type_tmp = 0
     956         124 :       kind_z_tmp = 0
     957             :       nk = 0
     958         124 :       DO iat = 1, nexc_atoms
     959          82 :          iatom = xas_env%exc_atoms(iat)
     960          82 :          NULLIFY (atomic_kind)
     961          82 :          atomic_kind => particle_set(iatom)%atomic_kind
     962          82 :          CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=ikind)
     963          82 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zatom)
     964          82 :          ihavethis = .FALSE.
     965         112 :          DO ik = 1, nk
     966         112 :             IF (ikind == kind_type_tmp(ik)) THEN
     967          10 :                ihavethis = .TRUE.
     968          10 :                xas_env%mykind_of_atom(iat) = ik
     969             :                EXIT
     970             :             END IF
     971             :          END DO
     972         124 :          IF (.NOT. ihavethis) THEN
     973          72 :             nk = nk + 1
     974          72 :             kind_type_tmp(nk) = ikind
     975          72 :             kind_z_tmp(nk) = INT(zatom)
     976          72 :             xas_env%mykind_of_atom(iat) = nk
     977          72 :             xas_env%mykind_of_kind(ikind) = nk
     978             :          END IF
     979             :       END DO ! iat
     980             : 
     981         198 :       ALLOCATE (xas_env%my_gto_basis(nk))
     982         198 :       ALLOCATE (xas_env%stogto_overlap(nk))
     983         114 :       DO ik = 1, nk
     984          72 :          NULLIFY (xas_env%my_gto_basis(ik)%gto_basis_set, sto_basis_set)
     985          72 :          ne = 0
     986         146 :          DO l = 1, lq(1) + 1
     987          74 :             nj = 2*(l - 1) + 1
     988         222 :             DO i = l, nq(1)
     989          76 :                ne(l, i) = ptable(kind_z_tmp(ik))%e_conv(l - 1) - 2*nj*(i - l)
     990          76 :                ne(l, i) = MAX(ne(l, i), 0)
     991         150 :                ne(l, i) = MIN(ne(l, i), 2*nj)
     992             :             END DO
     993             :          END DO
     994             : 
     995          72 :          sto_zet(1) = srules(kind_z_tmp(ik), ne, nq(1), lq(1))
     996          72 :          CALL allocate_sto_basis_set(sto_basis_set)
     997          72 :          name_sto = 'xas_tmp_sto'
     998             :          CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, &
     999          72 :                                 lq=lq, zet=sto_zet, name=name_sto)
    1000             :          CALL create_gto_from_sto_basis(sto_basis_set, &
    1001          72 :                                         xas_env%my_gto_basis(ik)%gto_basis_set, xas_control%ngauss)
    1002          72 :          CALL deallocate_sto_basis_set(sto_basis_set)
    1003          72 :          xas_env%my_gto_basis(ik)%gto_basis_set%norm_type = 2
    1004          72 :          CALL init_orb_basis_set(xas_env%my_gto_basis(ik)%gto_basis_set)
    1005             : 
    1006          72 :          ikind = kind_type_tmp(ik)
    1007          72 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
    1008             : 
    1009          72 :          CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf_gto)
    1010          72 :          CALL get_gto_basis_set(gto_basis_set=xas_env%my_gto_basis(ik)%gto_basis_set, nsgf=nsgf_sto)
    1011         288 :          ALLOCATE (xas_env%stogto_overlap(ik)%array(nsgf_sto, nsgf_gto))
    1012             : 
    1013             :          CALL calc_stogto_overlap(xas_env%my_gto_basis(ik)%gto_basis_set, orb_basis_set, &
    1014         186 :                                   xas_env%stogto_overlap(ik)%array)
    1015             :       END DO
    1016             : 
    1017          42 :       DEALLOCATE (nq, lq, sto_zet)
    1018          42 :       DEALLOCATE (kind_type_tmp, kind_z_tmp)
    1019             : 
    1020         126 :    END SUBROUTINE xas_env_init
    1021             : 
    1022             : ! **************************************************************************************************
    1023             : !> \brief Calculate and write the spectrum relative to the core level excitation
    1024             : !>      of a specific atom. It works for TP approach, because of the definition
    1025             : !>      of the oscillator strengths as  matrix elements of the dipole operator
    1026             : !> \param xas_control ...
    1027             : !> \param xas_env ...
    1028             : !> \param qs_env ...
    1029             : !> \param xas_section ...
    1030             : !> \param iatom index of the excited atom
    1031             : !> \param istate ...
    1032             : !> \par History
    1033             : !>      03.2006 created [MI]
    1034             : !> \author MI
    1035             : !> \note
    1036             : !>      for the tddft calculation should be re-thought
    1037             : ! **************************************************************************************************
    1038          54 :    SUBROUTINE cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, &
    1039             :                                      iatom, istate)
    1040             : 
    1041             :       TYPE(xas_control_type)                             :: xas_control
    1042             :       TYPE(xas_environment_type), POINTER                :: xas_env
    1043             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1044             :       TYPE(section_vals_type), POINTER                   :: xas_section
    1045             :       INTEGER, INTENT(IN)                                :: iatom, istate
    1046             : 
    1047             :       INTEGER                                            :: homo, i, lfomo, my_spin, nabs, nmo, &
    1048             :                                                             nvirtual, output_unit, xas_estate
    1049             :       LOGICAL                                            :: append_cube, length
    1050             :       REAL(dp)                                           :: rc(3)
    1051          54 :       REAL(dp), DIMENSION(:), POINTER                    :: all_evals
    1052             :       REAL(dp), DIMENSION(:, :), POINTER                 :: sp_ab, sp_em
    1053          54 :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: dip_fm_set
    1054             :       TYPE(cp_fm_type), POINTER                          :: all_vectors, excvec_coeff
    1055             :       TYPE(cp_logger_type), POINTER                      :: logger
    1056          54 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: op_sm, ostrength_sm
    1057          54 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1058          54 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1059             : 
    1060          54 :       NULLIFY (logger)
    1061         108 :       logger => cp_get_default_logger()
    1062          54 :       output_unit = cp_logger_get_default_io_unit(logger)
    1063             : 
    1064          54 :       NULLIFY (ostrength_sm, op_sm, dip_fm_set)
    1065          54 :       NULLIFY (all_evals, all_vectors, excvec_coeff)
    1066          54 :       NULLIFY (mos, particle_set, sp_em, sp_ab)
    1067         216 :       ALLOCATE (op_sm(3))
    1068             : 
    1069             :       CALL get_qs_env(qs_env=qs_env, &
    1070          54 :                       mos=mos, particle_set=particle_set)
    1071             : 
    1072             :       CALL get_xas_env(xas_env=xas_env, all_vectors=all_vectors, xas_estate=xas_estate, &
    1073             :                        all_evals=all_evals, dip_fm_set=dip_fm_set, excvec_coeff=excvec_coeff, &
    1074          54 :                        ostrength_sm=ostrength_sm, nvirtual=nvirtual, spin_channel=my_spin)
    1075          54 :       CALL get_mo_set(mos(my_spin), homo=homo, lfomo=lfomo, nmo=nmo)
    1076             : 
    1077          54 :       nabs = nvirtual - lfomo + 1
    1078         162 :       ALLOCATE (sp_em(6, homo))
    1079         162 :       ALLOCATE (sp_ab(6, nabs))
    1080          54 :       CPASSERT(ASSOCIATED(excvec_coeff))
    1081             : 
    1082          54 :       IF (.NOT. xas_control%xas_method == xas_dscf) THEN
    1083             :          ! Calculate the spectrum
    1084          54 :          IF (xas_control%dipole_form == xas_dip_len) THEN
    1085           0 :             rc(1:3) = particle_set(iatom)%r(1:3)
    1086           0 :             DO i = 1, 3
    1087           0 :                NULLIFY (op_sm(i)%matrix)
    1088           0 :                op_sm(i)%matrix => ostrength_sm(i)%matrix
    1089             :             END DO
    1090           0 :             CALL rRc_xyz_ao(op_sm, qs_env, rc, order=1, minimum_image=.TRUE.)
    1091             :             CALL spectrum_dip_vel(dip_fm_set, op_sm, mos, excvec_coeff, &
    1092             :                                   all_vectors, all_evals, &
    1093           0 :                                   sp_em, sp_ab, xas_estate, nvirtual, my_spin)
    1094           0 :             DO i = 1, SIZE(ostrength_sm, 1)
    1095           0 :                CALL dbcsr_set(ostrength_sm(i)%matrix, 0.0_dp)
    1096             :             END DO
    1097             :          ELSE
    1098         216 :             DO i = 1, 3
    1099         162 :                NULLIFY (op_sm(i)%matrix)
    1100         216 :                op_sm(i)%matrix => ostrength_sm(i)%matrix
    1101             :             END DO
    1102             :             CALL spectrum_dip_vel(dip_fm_set, op_sm, mos, excvec_coeff, &
    1103             :                                   all_vectors, all_evals, &
    1104          54 :                                   sp_em, sp_ab, xas_estate, nvirtual, my_spin)
    1105             :          END IF
    1106             :       END IF
    1107             : 
    1108          54 :       CALL get_mo_set(mos(my_spin), lfomo=lfomo)
    1109             :       ! write the spectrum, if the file exists it is appended
    1110          54 :       IF (.NOT. xas_control%xas_method == xas_dscf) THEN
    1111          54 :          length = (.NOT. xas_control%dipole_form == xas_dip_vel)
    1112             :          CALL xas_write(sp_em, sp_ab, xas_estate, &
    1113          54 :                         xas_section, iatom, istate, lfomo, length=length)
    1114             :       END IF
    1115             : 
    1116          54 :       DEALLOCATE (sp_em)
    1117          54 :       DEALLOCATE (sp_ab)
    1118             : 
    1119          54 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, xas_section, &
    1120             :                                            "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
    1121           4 :          append_cube = section_get_lval(xas_section, "PRINT%CLS_FUNCTION_CUBES%APPEND")
    1122             :          CALL xas_print_cubes(xas_control, qs_env, xas_section, mos, all_vectors, &
    1123           4 :                               iatom, append_cube)
    1124             :       END IF
    1125             : 
    1126          54 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, xas_section, &
    1127             :                                            "PRINT%PDOS"), cp_p_file)) THEN
    1128           4 :          CALL xas_pdos(qs_env, xas_section, mos, iatom)
    1129             :       END IF
    1130             : 
    1131          54 :       DEALLOCATE (op_sm)
    1132             : 
    1133         162 :    END SUBROUTINE cls_calculate_spectrum
    1134             : 
    1135             : ! **************************************************************************************************
    1136             : !> \brief write the spectrum for each atom in a different output file
    1137             : !> \param sp_em ...
    1138             : !> \param sp_ab ...
    1139             : !> \param estate ...
    1140             : !> \param xas_section ...
    1141             : !> \param iatom index of the excited atom
    1142             : !> \param state_to_be_excited ...
    1143             : !> \param lfomo ...
    1144             : !> \param length ...
    1145             : !> \par History
    1146             : !>      05.2005 created [MI]
    1147             : !> \author MI
    1148             : !> \note
    1149             : !>      the iteration counter is not finilized yet
    1150             : ! **************************************************************************************************
    1151          54 :    SUBROUTINE xas_write(sp_em, sp_ab, estate, xas_section, iatom, state_to_be_excited, &
    1152             :                         lfomo, length)
    1153             : 
    1154             :       REAL(dp), DIMENSION(:, :), POINTER                 :: sp_em, sp_ab
    1155             :       INTEGER, INTENT(IN)                                :: estate
    1156             :       TYPE(section_vals_type), POINTER                   :: xas_section
    1157             :       INTEGER, INTENT(IN)                                :: iatom, state_to_be_excited, lfomo
    1158             :       LOGICAL, INTENT(IN)                                :: length
    1159             : 
    1160             :       CHARACTER(LEN=default_string_length)               :: mittle_ab, mittle_em, my_act, my_pos
    1161             :       INTEGER                                            :: i, istate, out_sp_ab, out_sp_em
    1162             :       REAL(dp)                                           :: ene2
    1163             :       TYPE(cp_logger_type), POINTER                      :: logger
    1164             : 
    1165          54 :       NULLIFY (logger)
    1166          54 :       logger => cp_get_default_logger()
    1167             : 
    1168          54 :       my_pos = "APPEND"
    1169          54 :       my_act = "WRITE"
    1170             : 
    1171          54 :       mittle_em = "xes_at"//TRIM(ADJUSTL(cp_to_string(iatom)))//"_st"//TRIM(ADJUSTL(cp_to_string(state_to_be_excited)))
    1172             : 
    1173             :       out_sp_em = cp_print_key_unit_nr(logger, xas_section, "PRINT%XES_SPECTRUM", &
    1174             :                                        extension=".spectrum", file_position=my_pos, file_action=my_act, &
    1175          54 :                                        file_form="FORMATTED", middle_name=TRIM(mittle_em))
    1176             : 
    1177          54 :       IF (out_sp_em > 0) THEN
    1178          27 :          WRITE (out_sp_em, '(A,I6,A,I6,A,I6)') " Emission spectrum for atom ", iatom, &
    1179          54 :             ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_em, 2)
    1180          27 :          ene2 = 1.0_dp
    1181         318 :          DO istate = estate, SIZE(sp_em, 2)
    1182         291 :             IF (length) ene2 = sp_em(1, istate)*sp_em(1, istate)
    1183         291 :             WRITE (out_sp_em, '(I6,5F16.8,F10.5)') istate, sp_em(1, istate)*evolt, &
    1184         291 :                sp_em(2, istate)*ene2, sp_em(3, istate)*ene2, &
    1185         609 :                sp_em(4, istate)*ene2, sp_em(5, istate)*ene2, sp_em(6, istate)
    1186             :          END DO
    1187             :       END IF
    1188             :       CALL cp_print_key_finished_output(out_sp_em, logger, xas_section, &
    1189          54 :                                         "PRINT%XES_SPECTRUM")
    1190             : 
    1191          54 :       mittle_ab = "xas_at"//TRIM(ADJUSTL(cp_to_string(iatom)))//"_st"//TRIM(ADJUSTL(cp_to_string(state_to_be_excited)))
    1192             :       out_sp_ab = cp_print_key_unit_nr(logger, xas_section, "PRINT%XAS_SPECTRUM", &
    1193             :                                        extension=".spectrum", file_position=my_pos, file_action=my_act, &
    1194          54 :                                        file_form="FORMATTED", middle_name=TRIM(mittle_ab))
    1195             : 
    1196          54 :       IF (out_sp_ab > 0) THEN
    1197          21 :          WRITE (out_sp_ab, '(A,I6,A,I6,A,I6)') " Absorption spectrum for atom ", iatom, &
    1198          42 :             ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_ab, 2)
    1199          21 :          ene2 = 1.0_dp
    1200         809 :          DO i = 1, SIZE(sp_ab, 2)
    1201         788 :             istate = lfomo - 1 + i
    1202         788 :             IF (length) ene2 = sp_ab(1, i)*sp_ab(1, i)
    1203         788 :             WRITE (out_sp_ab, '(I6,5F16.8,F10.5)') istate, sp_ab(1, i)*evolt, &
    1204         788 :                sp_ab(2, i)*ene2, sp_ab(3, i)*ene2, &
    1205        1597 :                sp_ab(4, i)*ene2, sp_ab(5, i)*ene2, sp_ab(6, i)
    1206             :          END DO
    1207             :       END IF
    1208             : 
    1209             :       CALL cp_print_key_finished_output(out_sp_ab, logger, xas_section, &
    1210          54 :                                         "PRINT%XAS_SPECTRUM")
    1211             : 
    1212          54 :    END SUBROUTINE xas_write
    1213             : 
    1214             : ! **************************************************************************************************
    1215             : !> \brief write the cube files for a set of selected states
    1216             : !> \param xas_control provide number ant indexes of the states to be printed
    1217             : !> \param qs_env ...
    1218             : !> \param xas_section ...
    1219             : !> \param mos mos from which the states to be printed are extracted
    1220             : !> \param all_vectors ...
    1221             : !> \param iatom index of the atom that has been excited
    1222             : !> \param append_cube ...
    1223             : !> \par History
    1224             : !>      08.2005 created [MI]
    1225             : !> \author MI
    1226             : ! **************************************************************************************************
    1227           8 :    SUBROUTINE xas_print_cubes(xas_control, qs_env, xas_section, &
    1228           4 :                               mos, all_vectors, iatom, append_cube)
    1229             : 
    1230             :       TYPE(xas_control_type)                             :: xas_control
    1231             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1232             :       TYPE(section_vals_type), POINTER                   :: xas_section
    1233             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1234             :       TYPE(cp_fm_type), INTENT(IN)                       :: all_vectors
    1235             :       INTEGER, INTENT(IN)                                :: iatom
    1236             :       LOGICAL, INTENT(IN)                                :: append_cube
    1237             : 
    1238             :       CHARACTER(LEN=default_string_length)               :: my_mittle, my_pos
    1239             :       INTEGER                                            :: homo, istate0, my_spin, nspins, nstates
    1240           4 :       REAL(dp), DIMENSION(:, :), POINTER                 :: centers
    1241             :       TYPE(section_vals_type), POINTER                   :: print_key
    1242             : 
    1243           4 :       nspins = SIZE(mos)
    1244             : 
    1245           8 :       print_key => section_vals_get_subs_vals(xas_section, "PRINT%CLS_FUNCTION_CUBES")
    1246           4 :       my_mittle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom)))
    1247           4 :       nstates = SIZE(xas_control%list_cubes, 1)
    1248             : 
    1249           4 :       IF (xas_control%do_centers) THEN
    1250             :          ! one might like to calculate the centers of the xas orbital (without localizing them)
    1251             :       ELSE
    1252          12 :          ALLOCATE (centers(6, nstates))
    1253          88 :          centers = 0.0_dp
    1254             :       END IF
    1255           4 :       my_spin = xas_control%spin_channel
    1256             : 
    1257           4 :       CALL get_mo_set(mos(my_spin), homo=homo)
    1258           4 :       istate0 = 0
    1259             : 
    1260           4 :       my_pos = "REWIND"
    1261           4 :       IF (append_cube) THEN
    1262           0 :          my_pos = "APPEND"
    1263             :       END IF
    1264             : 
    1265             :       CALL qs_print_cubes(qs_env, all_vectors, nstates, xas_control%list_cubes, &
    1266           4 :                           centers, print_key, my_mittle, state0=istate0, file_position=my_pos)
    1267             : 
    1268           4 :       DEALLOCATE (centers)
    1269             : 
    1270           4 :    END SUBROUTINE xas_print_cubes
    1271             : 
    1272             : ! **************************************************************************************************
    1273             : !> \brief write the PDOS after the XAS SCF, i.e., with one excited core
    1274             : !> \param qs_env ...
    1275             : !> \param xas_section ...
    1276             : !> \param mos mos from which the eigenvalues and expansion coeffiecients are obtained
    1277             : !> \param iatom index of the atom that has been excited
    1278             : !> \par History
    1279             : !>      03.2016 created [MI]
    1280             : !> \author MI
    1281             : ! **************************************************************************************************
    1282             : 
    1283           4 :    SUBROUTINE xas_pdos(qs_env, xas_section, mos, iatom)
    1284             : 
    1285             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1286             :       TYPE(section_vals_type), POINTER                   :: xas_section
    1287             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1288             :       INTEGER, INTENT(IN)                                :: iatom
    1289             : 
    1290             :       CHARACTER(LEN=default_string_length)               :: xas_mittle
    1291             :       INTEGER                                            :: ispin
    1292           4 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1293           4 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1294           4 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1295             : 
    1296           4 :       NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
    1297           4 :       xas_mittle = 'xasat'//TRIM(ADJUSTL(cp_to_string(iatom)))//'_'
    1298             : 
    1299             :       CALL get_qs_env(qs_env, &
    1300             :                       atomic_kind_set=atomic_kind_set, &
    1301             :                       particle_set=particle_set, &
    1302           4 :                       qs_kind_set=qs_kind_set)
    1303             : 
    1304          12 :       DO ispin = 1, 2
    1305             :          CALL calculate_projected_dos(mos(ispin), atomic_kind_set, qs_kind_set, particle_set, qs_env, &
    1306          12 :                                       xas_section, ispin, xas_mittle)
    1307             :       END DO
    1308             : 
    1309           4 :    END SUBROUTINE xas_pdos
    1310             : ! **************************************************************************************************
    1311             : !> \brief Calculation of the spectrum when the dipole approximation
    1312             : !>      in the velocity form is used.
    1313             : !> \param fm_set components of the position operator in a full matrix form
    1314             : !>                already multiplied by the coefficiets
    1315             : !>                only the terms <C_i Op C_f> are calculated where
    1316             : !>                C_i are the coefficients of the excited state
    1317             : !> \param op_sm components of the position operator for the dipole
    1318             : !>               in a sparse matrix form (cos and sin)
    1319             : !>               calculated for the basis functions
    1320             : !> \param mos wavefunctions coefficients
    1321             : !> \param excvec coefficients of the excited orbital
    1322             : !> \param all_vectors ...
    1323             : !> \param all_evals ...
    1324             : !> \param sp_em ...
    1325             : !> \param sp_ab ...
    1326             : !> \param estate index of the excited state
    1327             : !> \param nstate ...
    1328             : !> \param my_spin ...
    1329             : !> \par History
    1330             : !>      06.2005 created [MI]
    1331             : !> \author MI
    1332             : ! **************************************************************************************************
    1333          54 :    SUBROUTINE spectrum_dip_vel(fm_set, op_sm, mos, excvec, &
    1334             :                                all_vectors, all_evals, sp_em, sp_ab, estate, nstate, my_spin)
    1335             : 
    1336             :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: fm_set
    1337             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: op_sm
    1338             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1339             :       TYPE(cp_fm_type), INTENT(IN)                       :: excvec, all_vectors
    1340             :       REAL(dp), DIMENSION(:), POINTER                    :: all_evals
    1341             :       REAL(dp), DIMENSION(:, :), POINTER                 :: sp_em, sp_ab
    1342             :       INTEGER, INTENT(IN)                                :: estate, nstate, my_spin
    1343             : 
    1344             :       INTEGER                                            :: homo, i, i_abs, istate, lfomo, nao, nmo
    1345             :       REAL(dp)                                           :: dip(3), ene_f, ene_i
    1346          54 :       REAL(dp), DIMENSION(:), POINTER                    :: eigenvalues, occupation_numbers
    1347             :       TYPE(cp_fm_type)                                   :: fm_work
    1348             : 
    1349           0 :       CPASSERT(ASSOCIATED(fm_set))
    1350          54 :       NULLIFY (eigenvalues, occupation_numbers)
    1351             : 
    1352             :       CALL get_mo_set(mos(my_spin), eigenvalues=eigenvalues, occupation_numbers=occupation_numbers, &
    1353          54 :                       nao=nao, nmo=nmo, homo=homo, lfomo=lfomo)
    1354             : 
    1355          54 :       CALL cp_fm_create(fm_work, all_vectors%matrix_struct)
    1356         216 :       DO i = 1, SIZE(fm_set, 2)
    1357         162 :          CALL cp_fm_set_all(fm_set(my_spin, i), 0.0_dp)
    1358         162 :          CALL cp_fm_set_all(fm_work, 0.0_dp)
    1359         162 :          CALL cp_dbcsr_sm_fm_multiply(op_sm(i)%matrix, all_vectors, fm_work, ncol=nstate)
    1360             :          CALL parallel_gemm("T", "N", 1, nstate, nao, 1.0_dp, excvec, &
    1361         216 :                             fm_work, 0.0_dp, fm_set(my_spin, i), b_first_col=1)
    1362             :       END DO
    1363          54 :       CALL cp_fm_release(fm_work)
    1364             : 
    1365        4282 :       sp_em = 0.0_dp
    1366       11758 :       sp_ab = 0.0_dp
    1367          54 :       ene_i = eigenvalues(estate)
    1368        2312 :       DO istate = 1, nstate
    1369        2258 :          ene_f = all_evals(istate)
    1370        9032 :          DO i = 1, 3
    1371        9032 :             CALL cp_fm_get_element(fm_set(my_spin, i), 1, istate, dip(i))
    1372             :          END DO
    1373        2258 :          IF (istate <= homo) THEN
    1374         604 :             sp_em(1, istate) = ene_f - ene_i
    1375         604 :             sp_em(2, istate) = dip(1)
    1376         604 :             sp_em(3, istate) = dip(2)
    1377         604 :             sp_em(4, istate) = dip(3)
    1378         604 :             sp_em(5, istate) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
    1379         604 :             sp_em(6, istate) = occupation_numbers(istate)
    1380             :          END IF
    1381        2312 :          IF (istate >= lfomo) THEN
    1382        1672 :             i_abs = istate - lfomo + 1
    1383        1672 :             sp_ab(1, i_abs) = ene_f - ene_i
    1384        1672 :             sp_ab(2, i_abs) = dip(1)
    1385        1672 :             sp_ab(3, i_abs) = dip(2)
    1386        1672 :             sp_ab(4, i_abs) = dip(3)
    1387        1672 :             sp_ab(5, i_abs) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
    1388        1672 :             IF (istate <= nmo) sp_ab(6, i_abs) = occupation_numbers(istate)
    1389             :          END IF
    1390             : 
    1391             :       END DO
    1392             : 
    1393          54 :    END SUBROUTINE spectrum_dip_vel
    1394             : 
    1395             : ! **************************************************************************************************
    1396             : !> \brief ...
    1397             : !> \param base_a ...
    1398             : !> \param base_b ...
    1399             : !> \param matrix ...
    1400             : ! **************************************************************************************************
    1401         140 :    SUBROUTINE calc_stogto_overlap(base_a, base_b, matrix)
    1402             : 
    1403             :       TYPE(gto_basis_set_type), POINTER                  :: base_a, base_b
    1404             :       REAL(dp), DIMENSION(:, :), POINTER                 :: matrix
    1405             : 
    1406             :       INTEGER                                            :: iset, jset, ldsab, maxcoa, maxcob, maxl, &
    1407             :                                                             maxla, maxlb, na, nb, nseta, nsetb, &
    1408             :                                                             nsgfa, nsgfb, sgfa, sgfb
    1409         140 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
    1410         140 :                                                             npgfb, nsgfa_set, nsgfb_set
    1411         140 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
    1412             :       REAL(dp)                                           :: rab(3)
    1413         140 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: sab, work
    1414         140 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, scon_a, scon_b, sphi_a, &
    1415         140 :                                                             sphi_b, zeta, zetb
    1416             : 
    1417         140 :       NULLIFY (la_max, la_min, lb_max, lb_min)
    1418         140 :       NULLIFY (npgfa, npgfb, nsgfa_set, nsgfb_set)
    1419         140 :       NULLIFY (first_sgfa, first_sgfb)
    1420         140 :       NULLIFY (rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
    1421             : 
    1422             :       CALL get_gto_basis_set(gto_basis_set=base_a, nsgf=nsgfa, nsgf_set=nsgfa_set, lmax=la_max, &
    1423             :                              lmin=la_min, npgf=npgfa, pgf_radius=rpgfa, &
    1424             :                              sphi=sphi_a, scon=scon_a, zet=zeta, first_sgf=first_sgfa, &
    1425         140 :                              maxco=maxcoa, nset=nseta, maxl=maxla)
    1426             : 
    1427             :       CALL get_gto_basis_set(gto_basis_set=base_b, nsgf=nsgfb, nsgf_set=nsgfb_set, lmax=lb_max, &
    1428             :                              lmin=lb_min, npgf=npgfb, pgf_radius=rpgfb, &
    1429             :                              sphi=sphi_b, scon=scon_b, zet=zetb, first_sgf=first_sgfb, &
    1430         140 :                              maxco=maxcob, nset=nsetb, maxl=maxlb)
    1431             :       ! Initialize and allocate
    1432         140 :       rab = 0.0_dp
    1433        5232 :       matrix = 0.0_dp
    1434             : 
    1435         140 :       ldsab = MAX(maxcoa, maxcob, nsgfa, nsgfb)
    1436         140 :       maxl = MAX(maxla, maxlb)
    1437             : 
    1438         560 :       ALLOCATE (sab(ldsab, ldsab))
    1439         560 :       ALLOCATE (work(ldsab, ldsab))
    1440             : 
    1441         280 :       DO iset = 1, nseta
    1442             : 
    1443         140 :          na = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
    1444         140 :          sgfa = first_sgfa(1, iset)
    1445             : 
    1446        1004 :          DO jset = 1, nsetb
    1447         724 :             nb = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
    1448         724 :             sgfb = first_sgfb(1, jset)
    1449             : 
    1450             :             CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
    1451             :                             lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
    1452         724 :                             rab, sab)
    1453             :             CALL contraction(sab, work, ca=scon_a(:, sgfa:), na=na, ma=nsgfa_set(iset), &
    1454         724 :                              cb=scon_b(:, sgfb:), nb=nb, mb=nsgfb_set(jset))
    1455         864 :             CALL block_add("IN", work, nsgfa_set(iset), nsgfb_set(jset), matrix, sgfa, sgfb)
    1456             : 
    1457             :          END DO ! jset
    1458             :       END DO ! iset
    1459         140 :       DEALLOCATE (sab, work)
    1460             : 
    1461         140 :    END SUBROUTINE calc_stogto_overlap
    1462             : 
    1463             : ! **************************************************************************************************
    1464             : !> \brief Starting from a set of mos, determine on which atom are centered
    1465             : !>      and if they are of the right type (1s,2s ...)
    1466             : !>      to be used in the specific core level spectrum calculation
    1467             : !>      The set of states need to be from the core, otherwise the
    1468             : !>      characterization of the type is not valid, since it assumes that
    1469             : !>      the orbital is localizad on a specific atom
    1470             : !>      It is probably reccomandable to run a localization cycle before
    1471             : !>      proceeding to the assignment of the type
    1472             : !>      The type is determined by computing the overalp with a
    1473             : !>      type specific, minimal, STO bais set
    1474             : !> \param xas_control ...
    1475             : !> \param xas_env ...
    1476             : !> \param localized_wfn_control ...
    1477             : !> \param qs_env ...
    1478             : !> \par History
    1479             : !>      03.2006 created [MI]
    1480             : !> \author MI
    1481             : ! **************************************************************************************************
    1482          42 :    SUBROUTINE cls_assign_core_states(xas_control, xas_env, localized_wfn_control, qs_env)
    1483             : 
    1484             :       TYPE(xas_control_type)                             :: xas_control
    1485             :       TYPE(xas_environment_type), POINTER                :: xas_env
    1486             :       TYPE(localized_wfn_control_type), POINTER          :: localized_wfn_control
    1487             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1488             : 
    1489             :       INTEGER                                            :: chosen_state, homo, i, iat, iatom, &
    1490             :                                                             ikind, isgf, istate, j, my_kind, &
    1491             :                                                             my_spin, nao, natom, nexc_atoms, &
    1492             :                                                             nexc_search, output_unit
    1493             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: first_sgf
    1494             :       INTEGER, DIMENSION(3)                              :: perd0
    1495          42 :       INTEGER, DIMENSION(:), POINTER                     :: atom_of_state, mykind_of_kind, &
    1496          42 :                                                             nexc_states, state_of_mytype, &
    1497          42 :                                                             type_of_state
    1498          42 :       INTEGER, DIMENSION(:, :), POINTER                  :: state_of_atom
    1499             :       REAL(dp)                                           :: component, dist, distmin, maxocc, ra(3), &
    1500             :                                                             rac(3), rc(3)
    1501          42 :       REAL(dp), DIMENSION(:), POINTER                    :: max_overlap, sto_state_overlap
    1502          42 :       REAL(dp), DIMENSION(:, :), POINTER                 :: centers_wfn
    1503             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
    1504             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1505             :       TYPE(cell_type), POINTER                           :: cell
    1506          42 :       TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER        :: stogto_overlap
    1507             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1508             :       TYPE(cp_logger_type), POINTER                      :: logger
    1509          42 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1510          42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1511          42 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1512             : 
    1513          42 :       NULLIFY (cell, mos, particle_set)
    1514          42 :       NULLIFY (atom_of_state, centers_wfn, mykind_of_kind, state_of_atom, nexc_states)
    1515          42 :       NULLIFY (stogto_overlap, type_of_state, max_overlap, qs_kind_set)
    1516          42 :       NULLIFY (state_of_mytype, type_of_state, sto_state_overlap)
    1517             : 
    1518          42 :       NULLIFY (logger)
    1519          84 :       logger => cp_get_default_logger()
    1520          42 :       output_unit = cp_logger_get_default_io_unit(logger)
    1521             : 
    1522             :       CALL get_qs_env(qs_env=qs_env, cell=cell, mos=mos, particle_set=particle_set, &
    1523          42 :                       qs_kind_set=qs_kind_set)
    1524             : 
    1525             :       ! The Berry operator can be used only for periodic systems
    1526             :       ! If an isolated system is used the periodicity is overimposed
    1527         168 :       perd0(1:3) = cell%perd(1:3)
    1528         168 :       cell%perd(1:3) = 1
    1529             : 
    1530             :       CALL get_xas_env(xas_env=xas_env, &
    1531             :                        centers_wfn=centers_wfn, atom_of_state=atom_of_state, &
    1532             :                        mykind_of_kind=mykind_of_kind, &
    1533             :                        type_of_state=type_of_state, state_of_atom=state_of_atom, &
    1534             :                        stogto_overlap=stogto_overlap, nexc_atoms=nexc_atoms, &
    1535          42 :                        spin_channel=my_spin, nexc_search=nexc_search, nexc_states=nexc_states)
    1536             : 
    1537          42 :       CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff, maxocc=maxocc, nao=nao, homo=homo)
    1538             : 
    1539             :       ! scratch array for the state
    1540         126 :       ALLOCATE (vecbuffer(1, nao))
    1541          42 :       natom = SIZE(particle_set)
    1542             : 
    1543         126 :       ALLOCATE (first_sgf(natom))
    1544          42 :       CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
    1545         126 :       ALLOCATE (sto_state_overlap(nexc_search))
    1546         126 :       ALLOCATE (max_overlap(natom))
    1547         166 :       max_overlap = 0.0_dp
    1548         126 :       ALLOCATE (state_of_mytype(natom))
    1549         166 :       state_of_mytype = 0
    1550         216 :       atom_of_state = 0
    1551         124 :       nexc_states = 1
    1552         592 :       state_of_atom = 0
    1553             : 
    1554          42 :       IF (xas_control%orbital_list(1) < 0) THEN !Checks for manually selected orbitals from the localized set
    1555             : 
    1556         198 :          DO istate = 1, nexc_search
    1557         158 :             centers_wfn(1, istate) = localized_wfn_control%centers_set(my_spin)%array(1, istate)
    1558         158 :             centers_wfn(2, istate) = localized_wfn_control%centers_set(my_spin)%array(2, istate)
    1559         158 :             centers_wfn(3, istate) = localized_wfn_control%centers_set(my_spin)%array(3, istate)
    1560             : 
    1561             :             ! Assign the state to the closest atom
    1562         158 :             distmin = 100.0_dp
    1563         518 :             DO iat = 1, nexc_atoms
    1564         360 :                iatom = xas_control%exc_atoms(iat)
    1565        1440 :                ra(1:3) = particle_set(iatom)%r(1:3)
    1566        1440 :                rc(1:3) = centers_wfn(1:3, istate)
    1567         360 :                rac = pbc(ra, rc, cell)
    1568         360 :                dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
    1569             : 
    1570         518 :                IF (dist < distmin) THEN
    1571             : 
    1572         238 :                   atom_of_state(istate) = iatom
    1573         238 :                   distmin = dist
    1574             :                END IF
    1575             :             END DO
    1576         198 :             IF (atom_of_state(istate) /= 0) THEN
    1577             :                !Character of the state
    1578             :                CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, istate, &
    1579         158 :                                         nao, 1, transpose=.TRUE.)
    1580             : 
    1581         158 :                iatom = atom_of_state(istate)
    1582             : 
    1583         158 :                NULLIFY (atomic_kind)
    1584         158 :                atomic_kind => particle_set(iatom)%atomic_kind
    1585             :                CALL get_atomic_kind(atomic_kind=atomic_kind, &
    1586         158 :                                     kind_number=ikind)
    1587             : 
    1588         158 :                my_kind = mykind_of_kind(ikind)
    1589             : 
    1590         158 :                sto_state_overlap(istate) = 0.0_dp
    1591         348 :                DO i = 1, SIZE(stogto_overlap(my_kind)%array, 1)
    1592         190 :                   component = 0.0_dp
    1593        3214 :                   DO j = 1, SIZE(stogto_overlap(my_kind)%array, 2)
    1594        3024 :                      isgf = first_sgf(iatom) + j - 1
    1595        3214 :                      component = component + stogto_overlap(my_kind)%array(i, j)*vecbuffer(1, isgf)
    1596             :                   END DO
    1597             :                   sto_state_overlap(istate) = sto_state_overlap(istate) + &
    1598         348 :                                               component*component
    1599             :                END DO
    1600             : 
    1601         158 :                IF (sto_state_overlap(istate) > max_overlap(iatom)) THEN
    1602          96 :                   state_of_mytype(iatom) = istate
    1603          96 :                   max_overlap(iatom) = sto_state_overlap(istate)
    1604             :                END IF
    1605             :             END IF
    1606             :          END DO ! istate
    1607             : 
    1608             :          ! Includes all states within the chosen threshold relative to the maximum overlap
    1609          40 :          IF (xas_control%overlap_threshold < 1) THEN
    1610           4 :             DO iat = 1, nexc_atoms
    1611           2 :                iatom = xas_control%exc_atoms(iat)
    1612          20 :                DO istate = 1, nexc_search
    1613          18 :                   IF (atom_of_state(istate) == iatom) THEN
    1614             :                      IF (sto_state_overlap(istate) > max_overlap(iatom)*xas_control%overlap_threshold &
    1615          16 :                          .AND. istate /= state_of_mytype(iat)) THEN
    1616           6 :                         nexc_states(iat) = nexc_states(iat) + 1
    1617           6 :                         state_of_atom(iat, nexc_states(iat)) = istate
    1618             :                      END IF
    1619             :                   END IF
    1620             :                END DO
    1621             :             END DO
    1622             :          END IF
    1623             : 
    1624             :          ! In the set of states, assign the index of the state to be excited for iatom
    1625          40 :          IF (output_unit > 0) THEN
    1626             :             WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") &
    1627          20 :                "List the atoms to be excited and the relative of MOs index "
    1628             :          END IF
    1629             : 
    1630         120 :          DO iat = 1, nexc_atoms
    1631          80 :             iatom = xas_env%exc_atoms(iat)
    1632          80 :             state_of_atom(iat, 1) = state_of_mytype(iatom) ! Place the state with maximum overlap first in the list
    1633          80 :             IF (output_unit > 0) THEN
    1634             :                WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A)", advance='NO') &
    1635          40 :                   'Atom: ', iatom, "MO index:"
    1636             :             END IF
    1637         166 :             DO istate = 1, nexc_states(iat)
    1638         166 :                IF (istate < nexc_states(iat)) THEN
    1639           6 :                   IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(I4)", advance='NO') state_of_atom(iat, istate)
    1640             :                ELSE
    1641          80 :                   IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(I4)") state_of_atom(iat, istate)
    1642             :                END IF
    1643             :             END DO
    1644         120 :             IF (state_of_atom(iat, 1) == 0 .OR. state_of_atom(iat, 1) .GT. homo) THEN
    1645           0 :                CPABORT("A wrong state has been selected for excitation, check the Wannier centers")
    1646             :             END IF
    1647             :          END DO
    1648             : 
    1649          40 :          IF (xas_control%overlap_threshold < 1) THEN
    1650           4 :             DO iat = 1, nexc_atoms
    1651           4 :                IF (output_unit > 0) THEN
    1652             :                   WRITE (UNIT=output_unit, FMT="(/,T10,A,I6)") &
    1653           1 :                      'Overlap integrals for Atom: ', iat
    1654           5 :                   DO istate = 1, nexc_states(iat)
    1655             :                      WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A,T38,f10.8)") &
    1656           5 :                         'State: ', state_of_atom(iat, istate), "Overlap:", sto_state_overlap(state_of_atom(iat, istate))
    1657             :                   END DO
    1658             :                END IF
    1659             :             END DO
    1660             :          END IF
    1661             : 
    1662         120 :          CALL reallocate(xas_env%state_of_atom, 1, nexc_atoms, 1, MAXVAL(nexc_states)) ! Scales down the 2d-array to the minimal size
    1663             : 
    1664             :       ELSE ! Manually selected orbital indices
    1665             : 
    1666             :          ! Reallocate nexc_states and state_of_atom to include any atom
    1667           2 :          CALL reallocate(xas_env%nexc_states, 1, natom)
    1668           2 :          CALL reallocate(xas_env%state_of_atom, 1, natom, 1, SIZE(xas_control%orbital_list))
    1669           2 :          CALL get_xas_env(xas_env, nexc_states=nexc_states, state_of_atom=state_of_atom)
    1670             : 
    1671          14 :          nexc_states = 0
    1672          30 :          state_of_atom = 0
    1673           2 :          nexc_atoms = natom !To include all possible atoms in the spectrum calculation
    1674             : 
    1675           6 :          DO istate = 1, SIZE(xas_control%orbital_list)
    1676             : 
    1677           4 :             chosen_state = xas_control%orbital_list(istate)
    1678           4 :             nexc_atoms = 1
    1679           4 :             centers_wfn(1, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(1, chosen_state)
    1680           4 :             centers_wfn(2, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(2, chosen_state)
    1681           4 :             centers_wfn(3, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(3, chosen_state)
    1682             : 
    1683           4 :             distmin = 100.0_dp
    1684          28 :             DO iat = 1, natom
    1685          96 :                ra(1:3) = particle_set(iat)%r(1:3)
    1686          96 :                rc(1:3) = centers_wfn(1:3, chosen_state)
    1687          24 :                rac = pbc(ra, rc, cell)
    1688          24 :                dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
    1689          28 :                IF (dist < distmin) THEN
    1690           6 :                   atom_of_state(chosen_state) = iat !?
    1691           6 :                   distmin = dist
    1692             :                END IF
    1693             :             END DO ! iat
    1694             : 
    1695           4 :             nexc_states(atom_of_state(chosen_state)) = nexc_states(atom_of_state(chosen_state)) + 1
    1696           6 :             state_of_atom(atom_of_state(chosen_state), nexc_states(atom_of_state(chosen_state))) = chosen_state
    1697             : 
    1698             :          END DO !istate
    1699             : 
    1700             :          ! In the set of states, assign the index of the state to be excited for iatom
    1701           2 :          IF (output_unit > 0) THEN
    1702             :             WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") &
    1703           1 :                "List the atoms to be excited and the relative of MOs index "
    1704             :          END IF
    1705             : 
    1706          14 :          DO iat = 1, natom
    1707          12 :             IF (output_unit > 0 .AND. state_of_atom(iat, 1) /= 0) THEN
    1708             :                WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A)", advance='NO') &
    1709           2 :                   'Atom: ', iat, "MO index:"
    1710           4 :                DO i = 1, nexc_states(iat)
    1711           4 :                   IF (i < nexc_states(iat)) THEN
    1712           0 :                      WRITE (UNIT=output_unit, FMT="(I4)", advance='NO') state_of_atom(iat, i)
    1713             :                   ELSE
    1714           2 :                      WRITE (UNIT=output_unit, FMT="(I4)") state_of_atom(iat, i)
    1715             :                   END IF
    1716             :                END DO
    1717             :             END IF
    1718          14 :             IF (state_of_atom(iat, 1) .GT. homo) THEN
    1719           0 :                CPABORT("A wrong state has been selected for excitation, check the Wannier centers")
    1720             :             END IF
    1721             :          END DO
    1722             : 
    1723          14 :          CALL reallocate(xas_env%state_of_atom, 1, natom, 1, MAXVAL(nexc_states)) ! Scales down the 2d-array to the minimal size
    1724             : 
    1725             :       END IF !Checks for manually selected orbitals from the localized set
    1726             : 
    1727             :       ! Set back the correct periodicity
    1728         168 :       cell%perd(1:3) = perd0(1:3)
    1729             : 
    1730          42 :       DEALLOCATE (vecbuffer)
    1731          42 :       DEALLOCATE (first_sgf)
    1732          42 :       DEALLOCATE (sto_state_overlap)
    1733          42 :       DEALLOCATE (max_overlap)
    1734          42 :       DEALLOCATE (state_of_mytype)
    1735             : 
    1736          84 :    END SUBROUTINE cls_assign_core_states
    1737             : 
    1738             : END MODULE xas_methods

Generated by: LCOV version 1.15