diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 5f0f61ffd0..9cbbce9438 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -89,7 +89,8 @@ MODULE input_cp2k_dft weight_type_unit, wfi_aspc_nr, wfi_frozen_method_nr, wfi_linear_p_method_nr, & wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, wfi_ps_method_nr, & wfi_use_guess_method_nr, wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, & - wfi_use_prev_wf_method_nr, xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, & + wfi_use_prev_wf_method_nr, xas_1s_type, xas_2p_type, xas_2s_type, xas_3d_type, & + xas_3p_type, xas_3s_type, xas_4d_type, xas_4f_type, xas_4p_type, xas_4s_type, xas_dip_len, & xas_dip_vel, xas_dscf, xas_none, xas_tp_fh, xas_tp_flex, xas_tp_hh, xas_tp_xfh, & xas_tp_xhh, xes_tp_val USE input_cp2k_almo, ONLY: create_almo_scf_section @@ -7585,9 +7586,11 @@ CONTAINS description="Type of the orbitals that are excited for the xas spectra calculation", & usage="STATE_TYPE 1S", & default_i_val=xas_1s_type, & - enum_c_vals=s2a("1S", "2S", "2P"), & - enum_desc=s2a("1s orbitals", "2s orbitals", "2p orbitals"), & - enum_i_vals=(/xas_1s_type, xas_2s_type, xas_2p_type/)) + enum_c_vals=s2a("1S", "2S", "2P", "3S", "3P", "3D", "4S", "4P", "4D", "4F"), & + enum_desc=s2a("1s orbitals", "2s orbitals", "2p orbitals", "3s orbitals", "3p orbitals", & + "3d orbitals", "4s orbitals", "4p orbitals", "4d orbitals", "4f orbitals"), & + enum_i_vals=(/xas_1s_type, xas_2s_type, xas_2p_type, xas_3s_type, xas_3p_type, xas_3d_type, & + xas_4s_type, xas_4p_type, xas_4d_type, xas_4f_type/)) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -7615,6 +7618,23 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, name="OVERLAP_THRESHOLD", & + description="Threshold for including more than one initial core excited state "// & + "per atom. The threshold is taken relative to the maximum overlap.", & + usage="OVERLAP_THRESHOLD 8.e-1", default_r_val=1.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="ORBITAL_LIST", & + variants=(/"ORBITAL_LIST"/), & + description="Indices of the localized orbitals to be excited. "// & + "This keyword can be repeated several times"// & + "(useful if you have to specify many indexes).", & + usage="ORBITAL_LIST {integer} {integer} .. {integer} ", & + n_var=-1, type_of_var=integer_t, repeats=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, name="ADDED_MOS", & description="Number of additional MOS added spin up only", & usage="ADDED_MOS {integer}", default_i_val=-1) diff --git a/src/xas_control.F b/src/xas_control.F index d0432ddc0b..e17c57c441 100644 --- a/src/xas_control.F +++ b/src/xas_control.F @@ -48,11 +48,12 @@ MODULE xas_control INTEGER :: ngauss INTEGER :: stride INTEGER, DIMENSION(:), POINTER :: exc_atoms + INTEGER, DIMENSION(:), POINTER :: orbital_list LOGICAL :: cubes, do_centers LOGICAL :: xas_restart INTEGER, DIMENSION(:), POINTER :: list_cubes ! - REAL(dp) :: eps_added + REAL(dp) :: eps_added, overlap_threshold REAL(dp) :: xes_core_occupation REAL(dp) :: xes_homo_occupation REAL(dp) :: nel_tot, xas_core_occupation @@ -87,7 +88,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'read_xas_control', & routineP = moduleN//':'//routineN - INTEGER :: i, ir, n_rep, nex_at + INTEGER :: i, ir, n_rep, nex_at, nex_st INTEGER, DIMENSION(:), POINTER :: list LOGICAL :: hempty, was_present @@ -172,6 +173,31 @@ CONTAINS CALL section_vals_val_get(xas_section, "NGAUSS", & i_val=xas_control%ngauss) + CALL section_vals_val_get(xas_section, "OVERLAP_THRESHOLD", & + r_val=xas_control%overlap_threshold) + + CALL section_vals_val_get(xas_section, "ORBITAL_LIST", & + n_rep_val=n_rep) + IF (n_rep > 0) THEN + nex_st = 0 + DO ir = 1, n_rep + NULLIFY (list) + CALL section_vals_val_get(xas_section, "ORBITAL_LIST", & + i_rep_val=ir, i_vals=list) + + IF (ASSOCIATED(list)) THEN + CALL reallocate(xas_control%orbital_list, 1, nex_st+SIZE(list)) + DO i = 1, SIZE(list) + xas_control%orbital_list(i+nex_st) = list(i) + END DO + nex_st = nex_st+SIZE(list) + END IF + END DO ! ir + ELSE + ALLOCATE (xas_control%orbital_list(1)) + xas_control%orbital_list(1) = -1 + END IF + END SUBROUTINE read_xas_control ! ************************************************************************************************** @@ -267,6 +293,7 @@ CONTAINS xas_control%xes_core_occupation = 1.0_dp xas_control%xes_homo_occupation = 1.0_dp NULLIFY (xas_control%exc_atoms) + NULLIFY (xas_control%orbital_list) xas_control%cubes = .FALSE. xas_control%do_centers = .FALSE. NULLIFY (xas_control%list_cubes) @@ -291,6 +318,9 @@ CONTAINS IF (ASSOCIATED(xas_control%exc_atoms)) THEN DEALLOCATE (xas_control%exc_atoms) END IF + IF (ASSOCIATED(xas_control%orbital_list)) THEN + DEALLOCATE (xas_control%orbital_list) + END IF IF (ASSOCIATED(xas_control%list_cubes)) THEN DEALLOCATE (xas_control%list_cubes) END IF @@ -299,7 +329,6 @@ CONTAINS END IF END SUBROUTINE xas_control_release - ! ************************************************************************************************** !> \brief ... !> \param xas_control ... diff --git a/src/xas_env_types.F b/src/xas_env_types.F index bf62d563b8..be4aa50555 100644 --- a/src/xas_env_types.F +++ b/src/xas_env_types.F @@ -53,10 +53,12 @@ MODULE xas_env_types !> \param nao number of atomic orbitals in the basis !> \param exc_state state that is now excited (this change atom by atom) !> \param nvirtual number of empy states to take into account for the spectrum -!> \param state_of_atom for each atom the state that have to be excited (global index) -!> dimension is the number of atoms to be excited +!> \param state_of_atom for each atom the states that have to be excited (global index) +!> dimension is the number of atoms to be excited by the largest number of included states !> \param atom_of_state atom to which each state is assigned, -!> dimension is the number of states occupied that might be excited +!> dimension is the number of states occupied that might be excited +!> \param nexc_states number of states to be excited per atom +!> dimension is the number of atoms to be excited !> \param type_of_state character of the state (1s,2s,2p...) !> \param spectrum for each excitation the energy and the oscillator strength !> \param centers_wfn for each wfn the center of charge (optimized by localization) @@ -76,12 +78,15 @@ MODULE xas_env_types INTEGER :: spin_channel INTEGER :: nvirtual, nvirtual2 INTEGER :: unoccupied_max_iter - INTEGER, DIMENSION(:), POINTER :: state_of_atom + INTEGER, DIMENSION(:), POINTER :: atom_of_state INTEGER, DIMENSION(:), POINTER :: type_of_state INTEGER, DIMENSION(:), POINTER :: mykind_of_atom INTEGER, DIMENSION(:), POINTER :: mykind_of_kind INTEGER, DIMENSION(:), POINTER :: exc_atoms + INTEGER, DIMENSION(:), POINTER :: nexc_states + INTEGER, DIMENSION(:, :), POINTER :: state_of_atom + REAL(dp) :: ip_energy, occ_estate, unoccupied_eps, xas_nelectron, homo_occ REAL(dp), DIMENSION(:), POINTER :: all_evals REAL(dp), DIMENSION(:), POINTER :: unoccupied_evals @@ -116,6 +121,7 @@ CONTAINS !> \param centers_wfn ... !> \param atom_of_state ... !> \param exc_atoms ... +!> \param nexc_states ... !> \param type_of_state ... !> \param mykind_of_atom ... !> \param mykind_of_kind ... @@ -145,8 +151,8 @@ CONTAINS !> \param scf_control ... ! ************************************************************************************************** SUBROUTINE get_xas_env(xas_env, iter_count, exc_state, nao, nvirtual, nvirtual2, & - centers_wfn, atom_of_state, exc_atoms, type_of_state, mykind_of_atom, mykind_of_kind, & - state_of_atom, spectrum, groundstate_coeff, ostrength_sm, & + centers_wfn, atom_of_state, exc_atoms, nexc_states, type_of_state, mykind_of_atom, & + mykind_of_kind, state_of_atom, spectrum, groundstate_coeff, ostrength_sm, & dip_fm_set, excvec_coeff, excvec_overlap, & unoccupied_orbs, unoccupied_evals, unoccupied_max_iter, unoccupied_eps, & all_vectors, all_evals, my_gto_basis, qs_loc_env, & @@ -157,9 +163,10 @@ CONTAINS INTEGER, INTENT(OUT), OPTIONAL :: iter_count, exc_state, nao, nvirtual, & nvirtual2 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: centers_wfn - INTEGER, DIMENSION(:), OPTIONAL, POINTER :: atom_of_state, exc_atoms, type_of_state, & - mykind_of_atom, mykind_of_kind, & - state_of_atom + INTEGER, DIMENSION(:), OPTIONAL, POINTER :: atom_of_state, exc_atoms, nexc_states, & + type_of_state, mykind_of_atom, & + mykind_of_kind + INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: state_of_atom REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: spectrum TYPE(cp_fm_p_type), DIMENSION(:), OPTIONAL, & POINTER :: groundstate_coeff @@ -198,6 +205,7 @@ CONTAINS IF (PRESENT(occ_estate)) occ_estate = xas_env%occ_estate IF (PRESENT(xas_estate)) xas_estate = xas_env%xas_estate IF (PRESENT(nexc_search)) nexc_search = xas_env%nexc_search + IF (PRESENT(nexc_states)) nexc_states => xas_env%nexc_states IF (PRESENT(spin_channel)) spin_channel = xas_env%spin_channel IF (PRESENT(nexc_atoms)) nexc_atoms = xas_env%nexc_atoms IF (PRESENT(unoccupied_eps)) unoccupied_eps = xas_env%unoccupied_eps @@ -311,7 +319,9 @@ CONTAINS NULLIFY (xas_env%ao_mo_fm_pools) NULLIFY (xas_env%my_gto_basis) - NULLIFY (xas_env%state_of_atom, xas_env%atom_of_state) + NULLIFY (xas_env%atom_of_state) + NULLIFY (xas_env%nexc_states) + NULLIFY (xas_env%state_of_atom) NULLIFY (xas_env%exc_atoms) NULLIFY (xas_env%excvec_coeff, xas_env%excvec_overlap) NULLIFY (xas_env%type_of_state, xas_env%mykind_of_atom) @@ -345,6 +355,7 @@ CONTAINS IF (xas_env%ref_count == 0) THEN DEALLOCATE (xas_env%state_of_atom, xas_env%atom_of_state) + DEALLOCATE (xas_env%nexc_states) DEALLOCATE (xas_env%type_of_state) DEALLOCATE (xas_env%mykind_of_atom) DEALLOCATE (xas_env%mykind_of_kind) diff --git a/src/xas_methods.F b/src/xas_methods.F index 110b531518..c8c4cb414b 100644 --- a/src/xas_methods.F +++ b/src/xas_methods.F @@ -51,6 +51,7 @@ MODULE xas_methods dbcsr_distribution_type, dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric USE input_constants, ONLY: & do_loc_none, state_loc_list, state_loc_range, xas_1s_type, xas_2p_type, xas_2s_type, & + xas_3d_type, xas_3p_type, xas_3s_type, xas_4d_type, xas_4f_type, xas_4p_type, xas_4s_type, & xas_dip_len, xas_dip_vel, xas_dscf, xas_tp_fh, xas_tp_flex, xas_tp_hh, xas_tp_xfh, & xas_tp_xhh, xes_tp_val USE input_section_types, ONLY: section_get_lval,& @@ -139,7 +140,7 @@ CONTAINS !> 05.2005 created [MI] !> \author MI !> \note -!> the iteration counter is not finilized yet +!> the iteration counter is not finalized yet !> only the transition potential approach is active !> the localization can be switched off, otherwise !> it uses by default the berry phase approach @@ -153,10 +154,11 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'xas', routineP = moduleN//':'//routineN - INTEGER :: handle, homo, i, iat, iatom, ispin, my_homo(2), my_nelectron(2), my_spin, nao, & - nexc_atoms, nexc_search, nmo, nspins, output_unit, state_to_be_excited + INTEGER :: handle, homo, i, iat, iatom, ispin, istate, my_homo(2), my_nelectron(2), my_spin, & + nao, nexc_atoms, nexc_search, nmo, nspins, output_unit, state_to_be_excited INTEGER, DIMENSION(2) :: added_mos - INTEGER, DIMENSION(:), POINTER :: state_of_atom + INTEGER, DIMENSION(:), POINTER :: nexc_states + INTEGER, DIMENSION(:, :), POINTER :: state_of_atom LOGICAL :: ch_method_flags, converged, my_uocc(2), & should_stop, skip_scf, & transition_potential @@ -196,9 +198,9 @@ CONTAINS output_unit = cp_logger_get_default_io_unit(logger) NULLIFY (xas_env, groundstate_coeff, ostrength_sm, op_sm) - NULLIFY (excvec_coeff, state_of_atom, qs_loc_env, cell, scf_env) + NULLIFY (excvec_coeff, qs_loc_env, cell, scf_env) NULLIFY (matrix_ks) - NULLIFY (all_vectors, state_of_atom, xas_control) + NULLIFY (all_vectors, state_of_atom, nexc_states, xas_control) NULLIFY (vecbuffer, op_sm, mo_coeff_b) NULLIFY (dft_section, xas_section, scf_section, loc_section, print_loc_section) dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") @@ -280,9 +282,9 @@ CONTAINS END IF CALL get_xas_env(xas_env=xas_env, & - state_of_atom=state_of_atom, all_vectors=all_vectors, & + all_vectors=all_vectors, & groundstate_coeff=groundstate_coeff, excvec_coeff=excvec_coeff, & - nexc_atoms=nexc_atoms, occ_estate=occ_estate, xas_nelectron=xas_nelectron, & + nexc_atoms=nexc_atoms, & spin_channel=my_spin) ! Set of states among which there is the state to be excited @@ -381,6 +383,8 @@ CONTAINS ! through the overlap with atomic-like states CALL cls_assign_core_states(xas_control, xas_env, qs_loc_env%localized_wfn_control, & qs_env) + CALL get_xas_env(xas_env=xas_env, & + state_of_atom=state_of_atom, nexc_states=nexc_states) IF (skip_scf) THEN CALL get_mo_set(mos(my_spin)%mo_set, mo_coeff=mo_coeff) @@ -393,7 +397,6 @@ CONTAINS ! copy the coefficients of the mos in a temporary fm with the right structure IF (transition_potential) THEN - ! Calculate the operator CALL get_xas_env(xas_env=xas_env, ostrength_sm=ostrength_sm) DO i = 1, 3 @@ -408,125 +411,128 @@ CONTAINS ! DO SCF if required DO iat = 1, nexc_atoms iatom = xas_env%exc_atoms(iat) - ! determine which state has to be excited in the global list - state_to_be_excited = state_of_atom(iat) - - ! Take the state_to_be_excited vector from the full set and copy into excvec_coeff - CALL get_mo_set(mos(my_spin)%mo_set, nmo=nmo) - CALL get_xas_env(xas_env, occ_estate=occ_estate, xas_nelectron=xas_nelectron) - tmp = xas_nelectron+1.0_dp-occ_estate - IF (nmo < tmp) & - CPABORT("CLS: the required method needs added_mos to the ground state") - ! If the restart file for this atom exists, the mos and the - ! occupation numbers are overwritten - ! It is necessary that the restart is for the same xas method - ! otherwise the number of electrons and the occupation numbers - ! may not be consistent - IF (xas_control%xas_restart) THEN - CALL xas_read_restart(xas_env, xas_section, qs_env, xas_control%xas_method, iatom, & - state_to_be_excited) - END IF - - CALL set_xas_env(xas_env=xas_env, xas_estate=state_to_be_excited) - CALL get_mo_set(mos(my_spin)%mo_set, mo_coeff=mo_coeff) - CPASSERT(ASSOCIATED(excvec_coeff)) - CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, state_to_be_excited, & - nao, 1, transpose=.TRUE.) - CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, & - nao, 1, transpose=.TRUE.) - - IF (transition_potential) THEN - - IF (.NOT. skip_scf) THEN - IF (output_unit > 0) THEN - WRITE (UNIT=output_unit, FMT='(/,T5,A)') REPEAT("-", 75) - IF (xas_control%xas_method == xas_dscf) THEN - WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') & - "START DeltaSCF for the first excited state from the core state of ATOM ", iatom - ELSE - WRITE (UNIT=output_unit, FMT='(/,T10,A,I6)') & - "Start Core Level Spectroscopy Calculation with TP approach for ATOM ", iatom - WRITE (UNIT=output_unit, FMT='(T10,A,T50,f10.4)') "Occupation of the core orbital", & - occ_estate - WRITE (UNIT=output_unit, FMT='(T10,A28,I3, T50,F10.4)') "Number of electrons in Spin ", & - my_spin, xas_nelectron - END IF - END IF - - CALL get_xas_env(xas_env=xas_env, scf_env=scf_env) - IF (.NOT. ASSOCIATED(scf_env)) THEN - CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section) -! Moved here from qs_scf_env_initialize to be able to have more scf_env - CALL set_xas_env(xas_env, scf_env=scf_env) - CALL scf_env_release(scf_env) - CALL get_xas_env(xas_env=xas_env, scf_env=scf_env) - ELSE - CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section) - ENDIF - - DO ispin = 1, SIZE(mos) - IF (ASSOCIATED(mos(ispin)%mo_set%mo_coeff_b)) THEN !fm->dbcsr - CALL copy_fm_to_dbcsr(mos(ispin)%mo_set%mo_coeff, & - mos(ispin)%mo_set%mo_coeff_b) !fm->dbcsr - ENDIF !fm->dbcsr - ENDDO !fm->dbcsr - - IF (.NOT. scf_env%skip_diis) THEN - IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN - CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis) - END IF - CALL qs_diis_b_clear(scf_env%scf_diis_buffer) - END IF - - CALL xas_do_tp_scf(dft_control, xas_env, iatom, scf_env, qs_env, & - xas_section, scf_section, converged, should_stop) - - CALL external_control(should_stop, "CLS", target_time=qs_env%target_time, & - start_time=qs_env%start_time) - IF (should_stop) THEN - CALL scf_env_cleanup(scf_env) - EXIT - END IF + DO istate = 1, nexc_states(iat) + ! determine which state has to be excited in the global list + state_to_be_excited = state_of_atom(iat, istate) + ! Take the state_to_be_excited vector from the full set and copy into excvec_coeff + CALL get_mo_set(mos(my_spin)%mo_set, nmo=nmo) + CALL get_xas_env(xas_env, occ_estate=occ_estate, xas_nelectron=xas_nelectron) + tmp = xas_nelectron+1.0_dp-occ_estate + IF (nmo < tmp) & + CPABORT("CLS: the required method needs added_mos to the ground state") + ! If the restart file for this atom exists, the mos and the + ! occupation numbers are overwritten + ! It is necessary that the restart is for the same xas method + ! otherwise the number of electrons and the occupation numbers + ! may not be consistent + IF (xas_control%xas_restart) THEN + CALL xas_read_restart(xas_env, xas_section, qs_env, xas_control%xas_method, iatom, & + state_to_be_excited, istate) END IF - ! SCF DONE + CALL set_xas_env(xas_env=xas_env, xas_estate=state_to_be_excited) + CALL get_mo_set(mos(my_spin)%mo_set, mo_coeff=mo_coeff) + CPASSERT(ASSOCIATED(excvec_coeff)) + CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, state_to_be_excited, & + nao, 1, transpose=.TRUE.) + CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, & + nao, 1, transpose=.TRUE.) + + IF (transition_potential) THEN + + IF (.NOT. skip_scf) THEN + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT='(/,T5,A)') REPEAT("-", 75) + IF (xas_control%xas_method == xas_dscf) THEN + WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') & + "START DeltaSCF for the first excited state from the core state of ATOM ", iatom + ELSE + WRITE (UNIT=output_unit, FMT='(/,T10,A,I6)') & + "Start Core Level Spectroscopy Calculation with TP approach for ATOM ", iatom + WRITE (UNIT=output_unit, FMT='(/,T10,A,I6,T34,A,T54,I6)') & + "Excited state", istate, "out of", nexc_states(iat) + WRITE (UNIT=output_unit, FMT='(T10,A,T50,f10.4)') "Occupation of the core orbital", & + occ_estate + WRITE (UNIT=output_unit, FMT='(T10,A28,I3, T50,F10.4)') "Number of electrons in Spin ", & + my_spin, xas_nelectron + END IF + END IF + + CALL get_xas_env(xas_env=xas_env, scf_env=scf_env) + IF (.NOT. ASSOCIATED(scf_env)) THEN + CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section) + ! Moved here from qs_scf_env_initialize to be able to have more scf_env + CALL set_xas_env(xas_env, scf_env=scf_env) + CALL scf_env_release(scf_env) + CALL get_xas_env(xas_env=xas_env, scf_env=scf_env) + ELSE + CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section) + ENDIF + + DO ispin = 1, SIZE(mos) + IF (ASSOCIATED(mos(ispin)%mo_set%mo_coeff_b)) THEN !fm->dbcsr + CALL copy_fm_to_dbcsr(mos(ispin)%mo_set%mo_coeff, & + mos(ispin)%mo_set%mo_coeff_b) !fm->dbcsr + ENDIF !fm->dbcsr + ENDDO !fm->dbcsr + + IF (.NOT. scf_env%skip_diis) THEN + IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN + CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis) + END IF + CALL qs_diis_b_clear(scf_env%scf_diis_buffer) + END IF + + CALL xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, & + xas_section, scf_section, converged, should_stop) + + CALL external_control(should_stop, "CLS", target_time=qs_env%target_time, & + start_time=qs_env%start_time) + IF (should_stop) THEN + CALL scf_env_cleanup(scf_env) + EXIT + END IF + + END IF + ! SCF DONE ! *** Write last wavefunction to screen *** - DO ispin = 1, nspins - CALL write_mo_set(mos(ispin)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & - dft_section) - ENDDO + DO ispin = 1, nspins + CALL write_mo_set(mos(ispin)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & + dft_section) + ENDDO - ELSE - ! Core level spectroscopy by TDDFT is not yet implemented - ! the states defined by the rotation are the ground state orbitals - ! the initial state from which I excite should be localized - ! I take the excitations from lumo to nmo - END IF - - IF (converged) THEN - CALL cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, & - iatom) - ELSE - IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') & - "SCF with core hole NOT converged for ATOM ", iatom - END IF - - IF (.NOT. skip_scf) THEN - ! Reset the initial core orbitals. - ! The valence orbitals are taken from the last SCF, - ! it should be a better initial guess - CALL get_qs_env(qs_env, mos=mos) - DO ispin = 1, nspins - CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo) - CALL cp_fm_to_fm(groundstate_coeff(ispin)%matrix, mos(ispin)%mo_set%mo_coeff, nmo, 1, 1) - END DO - IF (iat == nexc_atoms) THEN - CALL scf_env_cleanup(scf_env) - CALL scf_env_release(xas_env%scf_env) + ELSE + ! Core level spectroscopy by TDDFT is not yet implemented + ! the states defined by the rotation are the ground state orbitals + ! the initial state from which I excite should be localized + ! I take the excitations from lumo to nmo END IF - END IF + IF (converged) THEN + CALL cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, & + iatom, istate) + ELSE + IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') & + "SCF with core hole NOT converged for ATOM ", iatom + END IF + + IF (.NOT. skip_scf) THEN + ! Reset the initial core orbitals. + ! The valence orbitals are taken from the last SCF, + ! it should be a better initial guess + CALL get_qs_env(qs_env, mos=mos) + DO ispin = 1, nspins + CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo) + CALL cp_fm_to_fm(groundstate_coeff(ispin)%matrix, mos(ispin)%mo_set%mo_coeff, nmo, 1, 1) + END DO + IF (iat == nexc_atoms) THEN + CALL scf_env_cleanup(scf_env) + CALL scf_env_release(xas_env%scf_env) + END IF + END IF + + END DO ! istate END DO ! iat = 1,nexc_atoms ! END of Calculation @@ -581,7 +587,7 @@ CONTAINS CHARACTER(LEN=default_string_length) :: name_sto INTEGER :: homo, i, iat, iatom, ik, ikind, ispin, j, l, lfomo, my_spin, n_mo(2), n_rep, nao, & natom, ncubes, nelectron, nexc_atoms, nexc_search, nj, nk, nkind, nmo, nmoloc(2), & - norb(0:3), nsgf_gto, nsgf_sto, nspins, nvirtual, nvirtual2 + nsgf_gto, nsgf_sto, nspins, nvirtual, nvirtual2 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_type_tmp, kind_z_tmp, & last_sgf INTEGER, DIMENSION(4, 7) :: ne @@ -670,7 +676,8 @@ CONTAINS ALLOCATE (xas_env%centers_wfn(3, nexc_search)) ALLOCATE (xas_env%atom_of_state(nexc_search)) ALLOCATE (xas_env%type_of_state(nexc_search)) - ALLOCATE (xas_env%state_of_atom(nexc_atoms)) + ALLOCATE (xas_env%state_of_atom(nexc_atoms, nexc_search)) + ALLOCATE (xas_env%nexc_states(nexc_atoms)) ALLOCATE (xas_env%mykind_of_atom(nexc_atoms)) nkind = SIZE(atomic_kind_set, 1) ALLOCATE (xas_env%mykind_of_kind(nkind)) @@ -875,10 +882,13 @@ CONTAINS CALL qs_loc_env_release(qs_loc_env) CALL get_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env) loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE") + CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.TRUE., & do_xas=.TRUE., nloc_xas=nexc_search, spin_xas=my_spin) + IF (.NOT. qs_loc_env%do_localize) THEN qs_loc_env%localized_wfn_control%localization_method = do_loc_none + ELSE nmoloc = qs_loc_env%localized_wfn_control%nloc_states CALL set_loc_wfn_lists(qs_loc_env%localized_wfn_control, nmoloc, n_mo, nspins, my_spin) @@ -889,21 +899,37 @@ CONTAINS END IF !Type of state - norb = 0 ALLOCATE (nq(1), lq(1), sto_zet(1)) IF (xas_control%state_type == xas_1s_type) THEN - norb(0) = 1 nq(1) = 1 lq(1) = 0 ELSEIF (xas_control%state_type == xas_2s_type) THEN - norb(0) = 2 nq(1) = 2 lq(1) = 0 ELSEIF (xas_control%state_type == xas_2p_type) THEN - norb(0) = 2 - norb(1) = 1 nq(1) = 2 lq(1) = 1 + ELSEIF (xas_control%state_type == xas_3s_type) THEN + nq(1) = 3 + lq(1) = 0 + ELSEIF (xas_control%state_type == xas_3p_type) THEN + nq(1) = 3 + lq(1) = 1 + ELSEIF (xas_control%state_type == xas_3d_type) THEN + nq(1) = 3 + lq(1) = 2 + ELSEIF (xas_control%state_type == xas_4s_type) THEN + nq(1) = 4 + lq(1) = 0 + ELSEIF (xas_control%state_type == xas_4p_type) THEN + nq(1) = 4 + lq(1) = 1 + ELSEIF (xas_control%state_type == xas_4d_type) THEN + nq(1) = 4 + lq(1) = 2 + ELSEIF (xas_control%state_type == xas_4f_type) THEN + nq(1) = 4 + lq(1) = 3 ELSE CPABORT("XAS type of state not implemented") END IF @@ -942,9 +968,9 @@ CONTAINS DO ik = 1, nk NULLIFY (xas_env%my_gto_basis(ik)%gto_basis_set, sto_basis_set) ne = 0 - DO l = 1, 4 !lq(1)+1 + DO l = 1, lq(1)+1 nj = 2*(l-1)+1 - DO i = l, 7 ! nq(1) + DO i = l, nq(1) ne(l, i) = ptable(kind_z_tmp(ik))%e_conv(l-1)-2*nj*(i-l) ne(l, i) = MAX(ne(l, i), 0) ne(l, i) = MIN(ne(l, i), 2*nj) @@ -987,6 +1013,7 @@ CONTAINS !> \param qs_env ... !> \param xas_section ... !> \param iatom index of the excited atom +!> \param istate ... !> \par History !> 03.2006 created [MI] !> \author MI @@ -994,13 +1021,13 @@ CONTAINS !> for the tddft calculation should be re-thought ! ************************************************************************************************** SUBROUTINE cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, & - iatom) + iatom, istate) TYPE(xas_control_type) :: xas_control TYPE(xas_environment_type), POINTER :: xas_env TYPE(qs_environment_type), POINTER :: qs_env TYPE(section_vals_type), POINTER :: xas_section - INTEGER, INTENT(IN) :: iatom + INTEGER, INTENT(IN) :: iatom, istate CHARACTER(LEN=*), PARAMETER :: routineN = 'cls_calculate_spectrum', & routineP = moduleN//':'//routineN @@ -1067,11 +1094,11 @@ CONTAINS END IF CALL get_mo_set(mos(my_spin)%mo_set, lfomo=lfomo) - ! writw thw spectrum, if the file exists it is appended + ! write the spectrum, if the file exists it is appended IF (.NOT. xas_control%xas_method == xas_dscf) THEN length = (.NOT. xas_control%dipole_form == xas_dip_vel) CALL xas_write(sp_em, sp_ab, xas_estate, & - xas_section, iatom, lfomo, length=length) + xas_section, iatom, istate, lfomo, length=length) END IF DEALLOCATE (sp_em) @@ -1100,6 +1127,7 @@ CONTAINS !> \param estate ... !> \param xas_section ... !> \param iatom index of the excited atom +!> \param state_to_be_excited ... !> \param lfomo ... !> \param length ... !> \par History @@ -1108,13 +1136,13 @@ CONTAINS !> \note !> the iteration counter is not finilized yet ! ************************************************************************************************** - SUBROUTINE xas_write(sp_em, sp_ab, estate, xas_section, iatom, & + SUBROUTINE xas_write(sp_em, sp_ab, estate, xas_section, iatom, state_to_be_excited, & lfomo, length) REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab INTEGER, INTENT(IN) :: estate TYPE(section_vals_type), POINTER :: xas_section - INTEGER, INTENT(IN) :: iatom, lfomo + INTEGER, INTENT(IN) :: iatom, state_to_be_excited, lfomo LOGICAL, INTENT(IN) :: length CHARACTER(len=*), PARAMETER :: routineN = 'xas_write', routineP = moduleN//':'//routineN @@ -1130,7 +1158,8 @@ CONTAINS my_pos = "APPEND" my_act = "WRITE" - mittle_em = "xes_at"//TRIM(ADJUSTL(cp_to_string(iatom))) + mittle_em = "xes_at"//TRIM(ADJUSTL(cp_to_string(iatom)))//"_st"//TRIM(ADJUSTL(cp_to_string(state_to_be_excited))) + out_sp_em = cp_print_key_unit_nr(logger, xas_section, "PRINT%XES_SPECTRUM", & extension=".spectrum", file_position=my_pos, file_action=my_act, & file_form="FORMATTED", middle_name=TRIM(mittle_em)) @@ -1149,7 +1178,7 @@ CONTAINS CALL cp_print_key_finished_output(out_sp_em, logger, xas_section, & "PRINT%XES_SPECTRUM") - mittle_ab = "xas_at"//TRIM(ADJUSTL(cp_to_string(iatom))) + mittle_ab = "xas_at"//TRIM(ADJUSTL(cp_to_string(iatom)))//"_st"//TRIM(ADJUSTL(cp_to_string(state_to_be_excited))) out_sp_ab = cp_print_key_unit_nr(logger, xas_section, "PRINT%XAS_SPECTRUM", & extension=".spectrum", file_position=my_pos, file_action=my_act, & file_form="FORMATTED", middle_name=TRIM(mittle_ab)) @@ -1336,7 +1365,6 @@ CONTAINS sp_em = 0.0_dp sp_ab = 0.0_dp ene_i = eigenvalues(estate) - DO istate = 1, nstate ene_f = all_evals(istate) DO i = 1, 3 @@ -1419,7 +1447,6 @@ CONTAINS sgfa = first_sgfa(1, iset) DO jset = 1, nsetb - nb = npgfb(jset)*(ncoset(lb_max(jset))-ncoset(lb_min(jset)-1)) sgfb = first_sgfb(1, jset) @@ -1432,7 +1459,6 @@ CONTAINS END DO ! jset END DO ! iset - DEALLOCATE (sab, work) END SUBROUTINE calc_stogto_overlap @@ -1466,15 +1492,16 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'cls_assign_core_states', & routineP = moduleN//':'//routineN - INTEGER :: homo, i, iat, iatom, ikind, isgf, & - istate, j, my_kind, my_spin, nao, & - natom, nexc_atoms, nexc_search, & - output_unit + INTEGER :: chosen_state, homo, i, iat, iatom, & + ikind, isgf, istate, j, my_kind, & + my_spin, nao, natom, nexc_atoms, & + nexc_search, output_unit INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf INTEGER, DIMENSION(3) :: perd0 INTEGER, DIMENSION(:), POINTER :: atom_of_state, mykind_of_kind, & - state_of_atom, state_of_mytype, & + nexc_states, state_of_mytype, & type_of_state + INTEGER, DIMENSION(:, :), POINTER :: state_of_atom REAL(dp) :: component, dist, distmin, maxocc, ra(3), & rac(3), rc(3) REAL(dp), DIMENSION(:), POINTER :: max_overlap, sto_state_overlap @@ -1490,7 +1517,7 @@ CONTAINS TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set NULLIFY (cell, mos, particle_set) - NULLIFY (atom_of_state, centers_wfn, mykind_of_kind, state_of_atom) + NULLIFY (atom_of_state, centers_wfn, mykind_of_kind, state_of_atom, nexc_states) NULLIFY (stogto_overlap, type_of_state, max_overlap, qs_kind_set) NULLIFY (state_of_mytype, type_of_state, sto_state_overlap) @@ -1511,7 +1538,7 @@ CONTAINS mykind_of_kind=mykind_of_kind, & type_of_state=type_of_state, state_of_atom=state_of_atom, & stogto_overlap=stogto_overlap, nexc_atoms=nexc_atoms, & - spin_channel=my_spin, nexc_search=nexc_search) + spin_channel=my_spin, nexc_search=nexc_search, nexc_states=nexc_states) CALL get_mo_set(mos(my_spin)%mo_set, mo_coeff=mo_coeff, maxocc=maxocc, nao=nao, homo=homo) @@ -1527,74 +1554,181 @@ CONTAINS ALLOCATE (state_of_mytype(natom)) state_of_mytype = 0 atom_of_state = 0 + nexc_states = 1 + state_of_atom = 0 - DO istate = 1, nexc_search - centers_wfn(1, istate) = localized_wfn_control%centers_set(my_spin)%array(1, istate) - centers_wfn(2, istate) = localized_wfn_control%centers_set(my_spin)%array(2, istate) - centers_wfn(3, istate) = localized_wfn_control%centers_set(my_spin)%array(3, istate) + IF (xas_control%orbital_list(1) < 0) THEN !Checks for manually selected orbitals from the localized set - ! Assign the state to the closest atom - distmin = 100.0_dp + DO istate = 1, nexc_search + centers_wfn(1, istate) = localized_wfn_control%centers_set(my_spin)%array(1, istate) + centers_wfn(2, istate) = localized_wfn_control%centers_set(my_spin)%array(2, istate) + centers_wfn(3, istate) = localized_wfn_control%centers_set(my_spin)%array(3, istate) + + ! Assign the state to the closest atom + distmin = 100.0_dp + DO iat = 1, nexc_atoms + iatom = xas_control%exc_atoms(iat) + ra(1:3) = particle_set(iatom)%r(1:3) + rc(1:3) = centers_wfn(1:3, istate) + rac = pbc(ra, rc, cell) + dist = rac(1)*rac(1)+rac(2)*rac(2)+rac(3)*rac(3) + + IF (dist < distmin) THEN + + atom_of_state(istate) = iatom + distmin = dist + END IF + END DO + IF (atom_of_state(istate) /= 0) THEN + !Character of the state + CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, istate, & + nao, 1, transpose=.TRUE.) + + iatom = atom_of_state(istate) + + NULLIFY (atomic_kind) + atomic_kind => particle_set(iatom)%atomic_kind + CALL get_atomic_kind(atomic_kind=atomic_kind, & + kind_number=ikind) + + my_kind = mykind_of_kind(ikind) + + sto_state_overlap(istate) = 0.0_dp + DO i = 1, SIZE(stogto_overlap(my_kind)%array, 1) + component = 0.0_dp + DO j = 1, SIZE(stogto_overlap(my_kind)%array, 2) + isgf = first_sgf(iatom)+j-1 + component = component+stogto_overlap(my_kind)%array(i, j)*vecbuffer(1, isgf) + END DO + sto_state_overlap(istate) = sto_state_overlap(istate)+ & + component*component + END DO + + IF (sto_state_overlap(istate) > max_overlap(iatom)) THEN + state_of_mytype(iatom) = istate + max_overlap(iatom) = sto_state_overlap(istate) + END IF + END IF + END DO ! istate + + ! Includes all states within the chosen threshold relative to the maximum overlap + IF (xas_control%overlap_threshold < 1) THEN + DO iat = 1, nexc_atoms + iatom = xas_control%exc_atoms(iat) + DO istate = 1, nexc_search + IF (atom_of_state(istate) == iatom) THEN + IF (sto_state_overlap(istate) > max_overlap(iatom)*xas_control%overlap_threshold & + .AND. istate /= state_of_mytype(iat)) THEN + nexc_states(iat) = nexc_states(iat)+1 + state_of_atom(iat, nexc_states(iat)) = istate + END IF + END IF + END DO + END DO + END IF + + ! In the set of states, assign the index of the state to be excited for iatom + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") & + "List the atoms to be excited and the relative of MOs index " + END IF DO iat = 1, nexc_atoms - iatom = xas_control%exc_atoms(iat) - ra(1:3) = particle_set(iatom)%r(1:3) - rc(1:3) = centers_wfn(1:3, istate) - rac = pbc(ra, rc, cell) - dist = rac(1)*rac(1)+rac(2)*rac(2)+rac(3)*rac(3) - - IF (dist < distmin) THEN - atom_of_state(istate) = iatom - distmin = dist + iatom = xas_env%exc_atoms(iat) + state_of_atom(iat, 1) = state_of_mytype(iatom) ! Place the state with maximum overlap first in the list + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A)", advance='NO') & + 'Atom: ', iatom, "MO index:" + END IF + DO istate = 1, nexc_states(iat) + IF (istate < nexc_states(iat)) THEN + IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(I4)", advance='NO') state_of_atom(iat, istate) + ELSE + IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(I4)") state_of_atom(iat, istate) + END IF + END DO + IF (state_of_atom(iat, 1) == 0 .OR. state_of_atom(iat, 1) .GT. homo) THEN + CPABORT("A wrong state has been selected for excitation, check the Wannier centers") END IF END DO - IF (atom_of_state(istate) /= 0) THEN - !Character of the state - CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, istate, & - nao, 1, transpose=.TRUE.) - iatom = atom_of_state(istate) - - NULLIFY (atomic_kind) - atomic_kind => particle_set(iatom)%atomic_kind - CALL get_atomic_kind(atomic_kind=atomic_kind, & - kind_number=ikind) - - my_kind = mykind_of_kind(ikind) - - sto_state_overlap(istate) = 0.0_dp - DO i = 1, SIZE(stogto_overlap(my_kind)%array, 1) - component = 0.0_dp - DO j = 1, SIZE(stogto_overlap(my_kind)%array, 2) - isgf = first_sgf(iatom)+j-1 - component = component+stogto_overlap(my_kind)%array(i, j)*vecbuffer(1, isgf) - END DO - sto_state_overlap(istate) = sto_state_overlap(istate)+ & - component*component + IF (xas_control%overlap_threshold < 1) THEN + DO iat = 1, nexc_atoms + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T10,A,I6)") & + 'Overlap integrals for Atom: ', iat + DO istate = 1, nexc_states(iat) + WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A,T38,f10.8)") & + 'State: ', state_of_atom(iat, istate), "Overlap:", sto_state_overlap(state_of_atom(iat, istate)) + END DO + END IF END DO - IF (sto_state_overlap(istate) > max_overlap(iatom)) THEN - state_of_mytype(iatom) = istate - max_overlap(iatom) = sto_state_overlap(istate) - END IF END IF - END DO ! istate - ! In the set of states, assign the index of the state to be excited for iatom - IF (output_unit > 0) THEN - WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") & - "List the atoms to be excited and the relative of MOs index " - END IF - DO iat = 1, nexc_atoms - iatom = xas_env%exc_atoms(iat) - state_of_atom(iat) = state_of_mytype(iatom) + CALL reallocate(xas_env%state_of_atom, 1, nexc_atoms, 1, MAXVAL(nexc_states)) ! Scales down the 2d-array to the minimal size + + ELSE ! Manually selected orbital indices + + ! Reallocate nexc_states and state_of_atom to include any atom + CALL reallocate(xas_env%nexc_states, 1, natom) + CALL reallocate(xas_env%state_of_atom, 1, natom, 1, SIZE(xas_control%orbital_list)) + CALL get_xas_env(xas_env, nexc_states=nexc_states, state_of_atom=state_of_atom) + + nexc_states = 0 + state_of_atom = 0 + nexc_atoms = natom !To include all possible atoms in the spectrum calculation + + DO istate = 1, SIZE(xas_control%orbital_list) + + chosen_state = xas_control%orbital_list(istate) + nexc_atoms = 1 + centers_wfn(1, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(1, chosen_state) + centers_wfn(2, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(2, chosen_state) + centers_wfn(3, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(3, chosen_state) + + distmin = 100.0_dp + DO iat = 1, natom + ra(1:3) = particle_set(iat)%r(1:3) + rc(1:3) = centers_wfn(1:3, chosen_state) + rac = pbc(ra, rc, cell) + dist = rac(1)*rac(1)+rac(2)*rac(2)+rac(3)*rac(3) + IF (dist < distmin) THEN + atom_of_state(chosen_state) = iat !? + distmin = dist + END IF + END DO ! iat + + nexc_states(atom_of_state(chosen_state)) = nexc_states(atom_of_state(chosen_state))+1 + state_of_atom(atom_of_state(chosen_state), nexc_states(atom_of_state(chosen_state))) = chosen_state + + END DO !istate + + ! In the set of states, assign the index of the state to be excited for iatom IF (output_unit > 0) THEN - WRITE (UNIT=output_unit, FMT="(T10,A,I6,T32,A,I6)") 'Atom: ', iatom, & - "MO index", state_of_atom(iat) + WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") & + "List the atoms to be excited and the relative of MOs index " END IF - IF (state_of_atom(iat) == 0 .OR. state_of_atom(iat) .GT. homo) THEN - CPABORT("A wrong state has been selected for excitation, check the Wannier centers") - END IF - END DO + + DO iat = 1, natom + IF (output_unit > 0 .AND. state_of_atom(iat, 1) /= 0) THEN + WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A)", advance='NO') & + 'Atom: ', iat, "MO index:" + DO i = 1, nexc_states(iat) + IF (i < nexc_states(iat)) THEN + WRITE (UNIT=output_unit, FMT="(I4)", advance='NO') state_of_atom(iat, i) + ELSE + WRITE (UNIT=output_unit, FMT="(I4)") state_of_atom(iat, i) + END IF + END DO + END IF + IF (state_of_atom(iat, 1) .GT. homo) THEN + CPABORT("A wrong state has been selected for excitation, check the Wannier centers") + END IF + END DO + + CALL reallocate(xas_env%state_of_atom, 1, natom, 1, MAXVAL(nexc_states)) ! Scales down the 2d-array to the minimal size + + END IF !Checks for manually selected orbitals from the localized set ! Set back the correct periodicity cell%perd(1:3) = perd0(1:3) diff --git a/src/xas_restart.F b/src/xas_restart.F index d648098c7c..29c707ba66 100644 --- a/src/xas_restart.F +++ b/src/xas_restart.F @@ -89,18 +89,20 @@ CONTAINS !> \param xas_method ... !> \param iatom index of the absorbing atom !> \param estate index of the core-hole orbital +!> \param istate counter of excited states per atom !> error: !> \par History !> 09.2006 created [MI] !> \author MI ! ************************************************************************************************** - SUBROUTINE xas_read_restart(xas_env, xas_section, qs_env, xas_method, iatom, estate) + SUBROUTINE xas_read_restart(xas_env, xas_section, qs_env, xas_method, iatom, estate, istate) TYPE(xas_environment_type), POINTER :: xas_env TYPE(section_vals_type), POINTER :: xas_section TYPE(qs_environment_type), POINTER :: qs_env INTEGER, INTENT(IN) :: xas_method, iatom INTEGER, INTENT(OUT) :: estate + INTEGER, INTENT(IN) :: istate CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_read_restart', & routineP = moduleN//':'//routineN @@ -142,8 +144,8 @@ CONTAINS xas=.TRUE.) CALL xstring(filename, ia, ie) - filename = filename(ia:ie)//'-at'// & - TRIM(ADJUSTL(cp_to_string(iatom)))//'.rst' + filename = filename(ia:ie)//'-at'//TRIM(ADJUSTL(cp_to_string(iatom)))// & + '_st'//TRIM(ADJUSTL(cp_to_string(istate)))//'.rst' INQUIRE (FILE=filename, EXIST=file_exists) ! open file @@ -268,13 +270,14 @@ CONTAINS !> \param qs_env ... !> \param xas_method ... !> \param iatom ... +!> \param istate ... ! ************************************************************************************************** - SUBROUTINE xas_write_restart(xas_env, xas_section, qs_env, xas_method, iatom) + SUBROUTINE xas_write_restart(xas_env, xas_section, qs_env, xas_method, iatom, istate) TYPE(xas_environment_type), POINTER :: xas_env TYPE(section_vals_type), POINTER :: xas_section TYPE(qs_environment_type), POINTER :: qs_env - INTEGER, INTENT(IN) :: xas_method, iatom + INTEGER, INTENT(IN) :: xas_method, iatom, istate CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_write_restart', & routineP = moduleN//':'//routineN @@ -309,7 +312,7 @@ CONTAINS ! Open file rst_unit = -1 - my_middle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom))) + my_middle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom)))//'_st'//TRIM(ADJUSTL(cp_to_string(istate))) rst_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%RESTART", & extension=".rst", file_status="REPLACE", file_action="WRITE", & file_form="UNFORMATTED", middle_name=TRIM(my_middle)) @@ -479,6 +482,7 @@ CONTAINS a_max = ABS(vecbuffer2(1, i)) IF (a_max > b_max) THEN ic_max = i + b_max = a_max ENDIF END DO diff --git a/src/xas_tp_scf.F b/src/xas_tp_scf.F index 70639cb536..c834e17e99 100644 --- a/src/xas_tp_scf.F +++ b/src/xas_tp_scf.F @@ -119,6 +119,7 @@ CONTAINS !> \param dft_control ... !> \param xas_env the environment for XAS calculations !> \param iatom ... +!> \param istate keeps track of the state index for restart writing !> \param scf_env the scf_env where to perform the scf procedure !> \param qs_env the qs_env, the scf_env and xas_env live in !> \param xas_section ... @@ -129,12 +130,12 @@ CONTAINS !> 05.2005 created [MI] !> \author MI ! ************************************************************************************************** - SUBROUTINE xas_do_tp_scf(dft_control, xas_env, iatom, scf_env, qs_env, & + SUBROUTINE xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, & xas_section, scf_section, converged, should_stop) TYPE(dft_control_type), POINTER :: dft_control TYPE(xas_environment_type), POINTER :: xas_env - INTEGER, INTENT(IN) :: iatom + INTEGER, INTENT(IN) :: iatom, istate TYPE(qs_scf_env_type), POINTER :: scf_env TYPE(qs_environment_type), POINTER :: qs_env TYPE(section_vals_type), POINTER :: xas_section, scf_section @@ -235,11 +236,16 @@ CONTAINS cp_to_string(scf_env%method)) CASE (general_diag_method_nr) ! diagonalisation (default) + scf_env%iter_count = iter_count CALL general_eigenproblem(scf_env, mos, matrix_ks, & matrix_s, scf_control, scf_section, & diis_step) + +!dbg + WRITE (*, *) 'tp scf ', xas_env%xas_estate CALL find_excited_core_orbital(xas_env, mos, matrix_s) + ! set occupation for the respective spin channels IF (my_spin == 1) THEN hole_spin = 1 ! hole is generated in channel 1 @@ -323,7 +329,7 @@ CONTAINS ! ** Write restart file ** CALL xas_write_restart(xas_env, xas_section, qs_env, xas_control%xas_method, & - iatom) + iatom, istate) IF (exit_loop) THEN CALL timestop(handle2) @@ -532,6 +538,7 @@ CONTAINS sto_state_overlap = sto_state_overlap+ & component*component END DO ! i size + IF (sto_state_overlap .GT. max_overlap) THEN max_overlap = sto_state_overlap my_state = istate