1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Initialize the XAS orbitals for specific core excitations 8!> Either the GS orbitals are used as initial guess, or the 9!> xas mos are read from a previous calculation. 10!> In the latter case, the core-hole potetial should be the same. 11!> \note 12!> The restart with the same core-hole potential should be checked 13!> and a wrong restart should stop the program 14!> \par History 15!> created 09.2006 16!> \author MI (09.2006) 17! ************************************************************************************************** 18MODULE xas_restart 19 20 USE cp_control_types, ONLY: dft_control_type 21 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply 22 USE cp_files, ONLY: close_file,& 23 open_file 24 USE cp_fm_types, ONLY: cp_fm_create,& 25 cp_fm_get_info,& 26 cp_fm_get_submatrix,& 27 cp_fm_release,& 28 cp_fm_set_all,& 29 cp_fm_set_submatrix,& 30 cp_fm_type,& 31 cp_fm_write_unformatted 32 USE cp_gemm_interface, ONLY: cp_gemm 33 USE cp_log_handling, ONLY: cp_get_default_logger,& 34 cp_logger_type,& 35 cp_to_string 36 USE cp_output_handling, ONLY: cp_p_file,& 37 cp_print_key_finished_output,& 38 cp_print_key_generate_filename,& 39 cp_print_key_should_output,& 40 cp_print_key_unit_nr 41 USE cp_para_types, ONLY: cp_para_env_type 42 USE dbcsr_api, ONLY: dbcsr_p_type 43 USE input_section_types, ONLY: section_vals_type 44 USE kinds, ONLY: default_path_length,& 45 default_string_length,& 46 dp 47 USE message_passing, ONLY: mp_bcast 48 USE qs_environment_types, ONLY: get_qs_env,& 49 qs_environment_type 50 USE qs_ks_types, ONLY: qs_ks_did_change 51 USE qs_mixing_utils, ONLY: mixing_init 52 USE qs_mo_io, ONLY: wfn_restart_file_name 53 USE qs_mo_methods, ONLY: calculate_density_matrix 54 USE qs_mo_occupation, ONLY: set_mo_occupation 55 USE qs_mo_types, ONLY: get_mo_set,& 56 mo_set_p_type,& 57 set_mo_set 58 USE qs_rho_atom_types, ONLY: rho_atom_type 59 USE qs_rho_methods, ONLY: qs_rho_update_rho 60 USE qs_rho_types, ONLY: qs_rho_get,& 61 qs_rho_type 62 USE qs_scf_types, ONLY: qs_scf_env_type 63 USE scf_control_types, ONLY: scf_control_type 64 USE string_utilities, ONLY: xstring 65 USE xas_env_types, ONLY: get_xas_env,& 66 set_xas_env,& 67 xas_environment_type 68#include "./base/base_uses.f90" 69 70 IMPLICIT NONE 71 PRIVATE 72 73 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_restart' 74 75! *** Public subroutines *** 76 77 PUBLIC :: xas_read_restart, xas_write_restart, xas_initialize_rho, find_excited_core_orbital 78 79CONTAINS 80 81! ************************************************************************************************** 82!> \brief Set up for reading the restart 83!> corresponing to the excitation of iatom 84!> If the corresponding restart file does not exist 85!> the GS orbitals are used as initial guess 86!> \param xas_env ... 87!> \param xas_section input section for XAS calculations 88!> qs_env: 89!> \param qs_env ... 90!> \param xas_method ... 91!> \param iatom index of the absorbing atom 92!> \param estate index of the core-hole orbital 93!> \param istate counter of excited states per atom 94!> error: 95!> \par History 96!> 09.2006 created [MI] 97!> \author MI 98! ************************************************************************************************** 99 SUBROUTINE xas_read_restart(xas_env, xas_section, qs_env, xas_method, iatom, estate, istate) 100 101 TYPE(xas_environment_type), POINTER :: xas_env 102 TYPE(section_vals_type), POINTER :: xas_section 103 TYPE(qs_environment_type), POINTER :: qs_env 104 INTEGER, INTENT(IN) :: xas_method, iatom 105 INTEGER, INTENT(OUT) :: estate 106 INTEGER, INTENT(IN) :: istate 107 108 CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_read_restart', & 109 routineP = moduleN//':'//routineN 110 111 CHARACTER(LEN=default_path_length) :: filename 112 INTEGER :: group, handle, i, ia, ie, ispin, my_spin, nao, nao_read, nelectron, nexc_atoms, & 113 nexc_atoms_read, nexc_search, nexc_search_read, nmo, nmo_read, output_unit, rst_unit, & 114 source, xas_estate, xas_estate_read, xas_method_read 115 LOGICAL :: file_exists 116 REAL(dp) :: occ_estate, occ_estate_read, & 117 xas_nelectron, xas_nelectron_read 118 REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers 119 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_read, occ_read 120 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer 121 TYPE(cp_fm_type), POINTER :: mo_coeff 122 TYPE(cp_logger_type), POINTER :: logger 123 TYPE(cp_para_env_type), POINTER :: para_env 124 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 125 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 126 127 CALL timeset(routineN, handle) 128 129 file_exists = .FALSE. 130 rst_unit = -1 131 132 NULLIFY (eigenvalues, matrix_s, mos, occupation_numbers, vecbuffer) 133 NULLIFY (logger) 134 logger => cp_get_default_logger() 135 136 output_unit = cp_print_key_unit_nr(logger, xas_section, & 137 "PRINT%PROGRAM_RUN_INFO", extension=".Log") 138 139 CALL get_qs_env(qs_env=qs_env, para_env=para_env) 140 group = para_env%group 141 source = para_env%source 142 143 IF (para_env%ionode) THEN 144 CALL wfn_restart_file_name(filename, file_exists, xas_section, logger, & 145 xas=.TRUE.) 146 147 CALL xstring(filename, ia, ie) 148 filename = filename(ia:ie)//'-at'//TRIM(ADJUSTL(cp_to_string(iatom)))// & 149 '_st'//TRIM(ADJUSTL(cp_to_string(istate)))//'.rst' 150 151 INQUIRE (FILE=filename, EXIST=file_exists) 152 ! open file 153 IF (file_exists) THEN 154 155 CALL open_file(file_name=TRIM(filename), & 156 file_action="READ", & 157 file_form="UNFORMATTED", & 158 file_position="REWIND", & 159 file_status="OLD", & 160 unit_number=rst_unit) 161 162 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,T20,A,I5,/)") & 163 "Read restart file for atom ", iatom 164 165 ELSE IF (.NOT. file_exists) THEN 166 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,T10,A,I5,A,/)") & 167 "Restart file for atom ", iatom, & 168 " not available. Initialization done with GS orbitals" 169 END IF 170 END IF 171 CALL mp_bcast(file_exists, source, group) 172 173 CALL get_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_estate=xas_estate, & 174 xas_nelectron=xas_nelectron, nexc_search=nexc_search, & 175 nexc_atoms=nexc_atoms, spin_channel=my_spin) 176 177 IF (file_exists) THEN 178 CALL get_qs_env(qs_env=qs_env, mos=mos, matrix_s=matrix_s) 179 180 IF (rst_unit > 0) THEN 181 READ (rst_unit) xas_method_read 182 READ (rst_unit) nexc_search_read, nexc_atoms_read, occ_estate_read, xas_nelectron_read 183 READ (rst_unit) xas_estate_read 184 185 IF (xas_method_read /= xas_method) & 186 CPABORT(" READ XAS RESTART: restart with different XAS method is not possible. ") 187 IF (nexc_atoms_read /= nexc_atoms) & 188 CALL cp_abort(__LOCATION__, & 189 " READ XAS RESTART: restart with different excited atoms "// & 190 " is not possible. Start instead a new XAS run with the new set of atoms. ") 191 ENDIF 192 193 CALL mp_bcast(xas_estate_read, source, group) 194 CALL set_xas_env(xas_env=xas_env, xas_estate=xas_estate_read) 195 estate = xas_estate_read 196 197 CALL get_mo_set(mo_set=mos(my_spin)%mo_set, nao=nao) 198 ALLOCATE (vecbuffer(1, nao)) 199 200 DO ispin = 1, SIZE(mos) 201 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nmo=nmo, eigenvalues=eigenvalues, & 202 occupation_numbers=occupation_numbers, mo_coeff=mo_coeff, nelectron=nelectron) 203 eigenvalues = 0.0_dp 204 occupation_numbers = 0.0_dp 205 CALL cp_fm_set_all(mo_coeff, 0.0_dp) 206 IF (para_env%ionode) THEN 207 READ (rst_unit) nao_read, nmo_read 208 IF (nao /= nao_read) & 209 CPABORT("To change basis is not possible. ") 210 ALLOCATE (eig_read(nmo_read), occ_read(nmo_read)) 211 eig_read = 0.0_dp 212 occ_read = 0.0_dp 213 nmo = MIN(nmo, nmo_read) 214 READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read) 215 eigenvalues(1:nmo) = eig_read(1:nmo) 216 occupation_numbers(1:nmo) = occ_read(1:nmo) 217 IF (nmo_read > nmo) THEN 218 IF (occupation_numbers(nmo) >= EPSILON(0.0_dp)) & 219 CALL cp_warn(__LOCATION__, & 220 "The number of occupied MOs on the restart unit is larger than "// & 221 "the allocated MOs.") 222 223 END IF 224 DEALLOCATE (eig_read, occ_read) 225 ENDIF 226 CALL mp_bcast(eigenvalues, source, group) 227 CALL mp_bcast(occupation_numbers, source, group) 228 229 DO i = 1, nmo 230 IF (para_env%ionode) THEN 231 READ (rst_unit) vecbuffer 232 ELSE 233 vecbuffer(1, :) = 0.0_dp 234 END IF 235 CALL mp_bcast(vecbuffer, source, group) 236 CALL cp_fm_set_submatrix(mo_coeff, & 237 vecbuffer, 1, i, nao, 1, transpose=.TRUE.) 238 END DO 239 ! Skip extra MOs if there any 240 IF (para_env%ionode) THEN 241 DO i = nmo + 1, nmo_read 242 READ (rst_unit) vecbuffer 243 END DO 244 END IF 245 246 END DO ! ispin 247 248 DEALLOCATE (vecbuffer) 249 250! nspin = SIZE(mos,1) 251! DO ispin = 1,nspin 252! ! ortho so that one can restart for different positions (basis sets?) 253! NULLIFY(mo_coeff) 254! CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff,homo=homo) 255! CALL make_basis_sm(mo_coeff,homo,matrix_s(1)%matrix) 256! END DO 257 END IF !file_exist 258 259 IF (para_env%ionode) THEN 260 IF (file_exists) CALL close_file(unit_number=rst_unit) 261 END IF 262 263 CALL timestop(handle) 264 265 END SUBROUTINE xas_read_restart 266 267! ************************************************************************************************** 268!> \brief ... 269!> \param xas_env ... 270!> \param xas_section ... 271!> \param qs_env ... 272!> \param xas_method ... 273!> \param iatom ... 274!> \param istate ... 275! ************************************************************************************************** 276 SUBROUTINE xas_write_restart(xas_env, xas_section, qs_env, xas_method, iatom, istate) 277 278 TYPE(xas_environment_type), POINTER :: xas_env 279 TYPE(section_vals_type), POINTER :: xas_section 280 TYPE(qs_environment_type), POINTER :: qs_env 281 INTEGER, INTENT(IN) :: xas_method, iatom, istate 282 283 CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_write_restart', & 284 routineP = moduleN//':'//routineN 285 286 CHARACTER(LEN=default_path_length) :: filename 287 CHARACTER(LEN=default_string_length) :: my_middle 288 INTEGER :: handle, ispin, nao, nexc_atoms, & 289 nexc_search, nmo, output_unit, & 290 rst_unit, xas_estate 291 REAL(dp) :: occ_estate, xas_nelectron 292 REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers 293 TYPE(cp_fm_type), POINTER :: mo_coeff 294 TYPE(cp_logger_type), POINTER :: logger 295 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 296 TYPE(section_vals_type), POINTER :: print_key 297 298 CALL timeset(routineN, handle) 299 NULLIFY (mos, logger, print_key) 300 logger => cp_get_default_logger() 301 302 CALL get_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_estate=xas_estate, & 303 xas_nelectron=xas_nelectron, nexc_search=nexc_search, nexc_atoms=nexc_atoms) 304 305 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 306 xas_section, "PRINT%RESTART", used_print_key=print_key), & 307 cp_p_file)) THEN 308 309 output_unit = cp_print_key_unit_nr(logger, xas_section, & 310 "PRINT%PROGRAM_RUN_INFO", extension=".Log") 311 312 CALL get_qs_env(qs_env=qs_env, mos=mos) 313 314 ! Open file 315 rst_unit = -1 316 my_middle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom)))//'_st'//TRIM(ADJUSTL(cp_to_string(istate))) 317 rst_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%RESTART", & 318 extension=".rst", file_status="REPLACE", file_action="WRITE", & 319 file_form="UNFORMATTED", middle_name=TRIM(my_middle)) 320 321 filename = cp_print_key_generate_filename(logger, print_key, & 322 middle_name=TRIM(my_middle), extension=".rst", & 323 my_local=.FALSE.) 324 325 IF (output_unit > 0) THEN 326 WRITE (UNIT=output_unit, FMT="(/,T10,A,I5,A,A,/)") & 327 "Xas orbitals for the absorbing atom ", iatom, & 328 " are written in ", TRIM(filename) 329 330 END IF 331 332 ! Write mos 333 IF (rst_unit > 0) THEN 334 WRITE (rst_unit) xas_method 335 WRITE (rst_unit) nexc_search, nexc_atoms, occ_estate, xas_nelectron 336 WRITE (rst_unit) xas_estate 337 ENDIF 338 DO ispin = 1, SIZE(mos) 339 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo, & 340 eigenvalues=eigenvalues, occupation_numbers=occupation_numbers) 341 IF ((rst_unit > 0)) THEN 342 WRITE (rst_unit) nao, nmo 343 WRITE (rst_unit) eigenvalues(1:nmo), & 344 occupation_numbers(1:nmo) 345 END IF 346 CALL cp_fm_write_unformatted(mo_coeff, rst_unit) 347 END DO 348 349! Close file 350 CALL cp_print_key_finished_output(rst_unit, logger, xas_section, & 351 "PRINT%RESTART") 352 END IF 353 CALL timestop(handle) 354 355 END SUBROUTINE xas_write_restart 356 357!****f* xas_restart/xas_initialize_rho [1.0] * 358 359! ************************************************************************************************** 360!> \brief Once the mos and the occupation numbers are initialized 361!> the electronic density of the excited state can be calclated 362!> \param qs_env ... 363!> \param scf_env ... 364!> \param scf_control ... 365!> \par History 366!> 09-2006 MI created 367!> \author MI 368! ************************************************************************************************** 369 SUBROUTINE xas_initialize_rho(qs_env, scf_env, scf_control) 370 371 TYPE(qs_environment_type), POINTER :: qs_env 372 TYPE(qs_scf_env_type), POINTER :: scf_env 373 TYPE(scf_control_type), POINTER :: scf_control 374 375 CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_initialize_rho', & 376 routineP = moduleN//':'//routineN 377 378 INTEGER :: handle, ispin, my_spin, nelectron 379 TYPE(cp_para_env_type), POINTER :: para_env 380 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao 381 TYPE(dft_control_type), POINTER :: dft_control 382 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 383 TYPE(qs_rho_type), POINTER :: rho 384 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom 385 TYPE(xas_environment_type), POINTER :: xas_env 386 387 CALL timeset(routineN, handle) 388 389 NULLIFY (mos, rho, xas_env, para_env, rho_ao) 390 391 CALL get_qs_env(qs_env, & 392 mos=mos, & 393 rho=rho, & 394 xas_env=xas_env, & 395 para_env=para_env) 396 397 my_spin = xas_env%spin_channel 398 CALL qs_rho_get(rho, rho_ao=rho_ao) 399 DO ispin = 1, SIZE(mos) 400 IF (ispin == my_spin) THEN 401 IF (xas_env%homo_occ == 0) THEN 402 CALL get_mo_set(mos(ispin)%mo_set, nelectron=nelectron) 403 nelectron = nelectron - 1 404 CALL set_mo_set(mos(ispin)%mo_set, nelectron=nelectron) 405 END IF 406 CALL set_mo_occupation(mo_set=qs_env%mos(ispin)%mo_set, smear=scf_control%smear, & 407 xas_env=xas_env) 408 ELSE 409 CALL set_mo_occupation(mo_set=qs_env%mos(ispin)%mo_set, smear=scf_control%smear) 410 END IF 411 CALL calculate_density_matrix(mo_set=mos(ispin)%mo_set, & 412 density_matrix=rho_ao(ispin)%matrix) 413 END DO 414 415 CALL qs_rho_update_rho(rho, qs_env=qs_env) 416 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 417 418 IF (scf_env%mixing_method > 1) THEN 419 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 420 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN 421 CPABORT('TB Code not available') 422 ELSE IF (dft_control%qs_control%semi_empirical) THEN 423 CPABORT('SE Code not possible') 424 ELSE 425 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom) 426 CALL mixing_init(scf_env%mixing_method, rho, scf_env%mixing_store, & 427 para_env, rho_atom=rho_atom) 428 END IF 429 END IF 430 431 CALL timestop(handle) 432 433 END SUBROUTINE xas_initialize_rho 434 435! ************************************************************************************************** 436!> \brief Find the index of the core orbital that has been excited by XAS 437!> \param xas_env ... 438!> \param mos ... 439!> \param matrix_s ... 440!> \par History 441!> 03-2010 MI created 442!> \author MI 443! ************************************************************************************************** 444 445 SUBROUTINE find_excited_core_orbital(xas_env, mos, matrix_s) 446 447 TYPE(xas_environment_type), POINTER :: xas_env 448 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 449 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 450 451 CHARACTER(LEN=*), PARAMETER :: routineN = 'find_excited_core_orbital', & 452 routineP = moduleN//':'//routineN 453 454 INTEGER :: i, ic_max, ir_max, m, my_spin, n, nao, & 455 nexc_search, nmo, xas_estate 456 INTEGER, DIMENSION(:), POINTER :: col_indices 457 REAL(dp) :: a_max, b_max, ip_energy, occ_estate 458 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers 459 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer, vecbuffer2 460 TYPE(cp_fm_type), POINTER :: excvec_coeff, excvec_overlap, fm_work, & 461 mo_coeff 462 463 NULLIFY (excvec_coeff, excvec_overlap, fm_work, mo_coeff) 464 ! Some elements from the xas_env 465 CALL get_xas_env(xas_env=xas_env, excvec_coeff=excvec_coeff, & 466 excvec_overlap=excvec_overlap, nexc_search=nexc_search, & 467 xas_estate=xas_estate, occ_estate=occ_estate, spin_channel=my_spin) 468 CPASSERT(ASSOCIATED(excvec_overlap)) 469 470 CALL get_mo_set(mos(my_spin)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo, & 471 eigenvalues=eigenvalues, occupation_numbers=occupation_numbers) 472 ALLOCATE (vecbuffer(1, nao)) 473 vecbuffer = 0.0_dp 474 ALLOCATE (vecbuffer2(1, nexc_search)) 475 vecbuffer2 = 0.0_dp 476 477 ! ** use the maximum overlap criterion to find the index of the excited orbital 478 CALL cp_fm_create(fm_work, mo_coeff%matrix_struct) 479 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, fm_work, ncol=nmo) 480 CALL cp_gemm("T", "N", 1, xas_env%nexc_search, nao, 1.0_dp, excvec_coeff, & 481 fm_work, 0.0_dp, excvec_overlap, b_first_col=1) 482 CALL cp_fm_get_info(matrix=excvec_overlap, col_indices=col_indices, & 483 nrow_global=m, ncol_global=n) 484 CALL cp_fm_get_submatrix(excvec_overlap, vecbuffer2, 1, 1, & 485 1, nexc_search, transpose=.FALSE.) 486 CALL cp_fm_release(fm_work) 487 488 b_max = 0.0_dp 489 ic_max = xas_estate 490 DO i = 1, nexc_search 491 a_max = ABS(vecbuffer2(1, i)) 492 IF (a_max > b_max) THEN 493 ic_max = i 494 495 b_max = a_max 496 ENDIF 497 END DO 498 499 IF (ic_max /= xas_estate) THEN 500 ir_max = xas_estate 501 xas_estate = ic_max 502 occupation_numbers(xas_estate) = occ_estate 503 occupation_numbers(ir_max) = 1.0_dp 504 END IF 505 506 ! Ionization Potential 507 iP_energy = eigenvalues(xas_estate) 508 CALL set_xas_env(xas_env=xas_env, xas_estate=xas_estate, ip_energy=ip_energy) 509 510 CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, xas_estate, & 511 nao, 1, transpose=.TRUE.) 512 CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, & 513 nao, 1, transpose=.TRUE.) 514 515 DEALLOCATE (vecbuffer, vecbuffer2) 516 517 END SUBROUTINE find_excited_core_orbital 518 519END MODULE xas_restart 520