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