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