comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Extension to the XAS code based on the Slater transition potential method.

Browse source 
The extension allows multiple initial core excited states, of any type
(s, p, d, ..). Code by Axel Erbung (Dep. Physics Stockholm)



svn-origin-rev: 18514
This commit is contained in:
Marcella Iannuzzi 2018-07-07 10:51:09 +00:00
parent 5821553b6e
commit 6c66e18b3b
6 changed files with 435 additions and 230 deletions
Download patch file Download diff file Expand all files Collapse all files

28
src/input_cp2k_dft.F
View file

@ -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)

35
src/xas_control.F
View file

@ -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 ...

31
src/xas_env_types.F
View file

@ -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)

542
src/xas_methods.F
View file

@ -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)

16
src/xas_restart.F
View file

@ -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

13
src/xas_tp_scf.F
View file

@ -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