Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : ELSE IF (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) THEN
435 0 : CPABORT("CLS: the required method needs added_mos to the ground state")
436 : END IF
437 : ! If the restart file for this atom exists, the mos and the
438 : ! occupation numbers are overwritten
439 : ! It is necessary that the restart is for the same xas method
440 : ! otherwise the number of electrons and the occupation numbers
441 : ! may not be consistent
442 88 : IF (xas_control%xas_restart) THEN
443 : CALL xas_read_restart(xas_env, xas_section, qs_env, xas_control%xas_method, iatom, &
444 12 : state_to_be_excited, istate)
445 : END IF
446 88 : CALL set_xas_env(xas_env=xas_env, xas_estate=state_to_be_excited)
447 88 : CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff)
448 88 : CPASSERT(ASSOCIATED(excvec_coeff))
449 : CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, state_to_be_excited, &
450 88 : nao, 1, transpose=.TRUE.)
451 : CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, &
452 88 : nao, 1, transpose=.TRUE.)
453 :
454 : IF (transition_potential) THEN
455 :
456 88 : IF (.NOT. skip_scf) THEN
457 80 : IF (output_unit > 0) THEN
458 40 : WRITE (UNIT=output_unit, FMT='(/,T5,A)') REPEAT("-", 75)
459 40 : IF (xas_control%xas_method == xas_dscf) THEN
460 : WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') &
461 0 : "START DeltaSCF for the first excited state from the core state of ATOM ", iatom
462 : ELSE
463 : WRITE (UNIT=output_unit, FMT='(/,T10,A,I6)') &
464 40 : "Start Core Level Spectroscopy Calculation with TP approach for ATOM ", iatom
465 : WRITE (UNIT=output_unit, FMT='(/,T10,A,I6,T34,A,T54,I6)') &
466 40 : "Excited state", istate, "out of", nexc_states(iat)
467 40 : WRITE (UNIT=output_unit, FMT='(T10,A,T50,f10.4)') "Occupation of the core orbital", &
468 80 : occ_estate
469 40 : WRITE (UNIT=output_unit, FMT='(T10,A28,I3, T50,F10.4)') "Number of electrons in Spin ", &
470 80 : my_spin, xas_nelectron
471 : END IF
472 : END IF
473 :
474 80 : CALL get_xas_env(xas_env=xas_env, scf_env=scf_env)
475 80 : IF (.NOT. ASSOCIATED(scf_env)) THEN
476 44 : CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
477 : ! Moved here from qs_scf_env_initialize to be able to have more scf_env
478 44 : CALL set_xas_env(xas_env, scf_env=scf_env)
479 : ELSE
480 36 : CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
481 : END IF
482 :
483 240 : DO ispin = 1, SIZE(mos)
484 240 : IF (ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN !fm->dbcsr
485 : CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
486 160 : mos(ispin)%mo_coeff_b) !fm->dbcsr
487 : END IF !fm->dbcsr
488 : END DO !fm->dbcsr
489 :
490 80 : IF (.NOT. scf_env%skip_diis) THEN
491 68 : IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
492 40 : ALLOCATE (scf_env%scf_diis_buffer)
493 40 : CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
494 : END IF
495 68 : CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
496 : END IF
497 :
498 : CALL xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, &
499 80 : xas_section, scf_section, converged, should_stop)
500 :
501 : CALL external_control(should_stop, "CLS", target_time=qs_env%target_time, &
502 80 : start_time=qs_env%start_time)
503 80 : IF (should_stop) THEN
504 0 : CALL scf_env_cleanup(scf_env)
505 0 : EXIT
506 : END IF
507 :
508 : END IF
509 : ! SCF DONE
510 :
511 : ! Write last wavefunction to screen
512 88 : IF (SIZE(mos) > 1) THEN
513 : CALL write_mo_set_to_output_unit(mos(1), qs_kind_set, particle_set, dft_section, &
514 88 : 4, 0, final_mos=.FALSE., spin="XAS ALPHA")
515 : CALL write_mo_set_to_output_unit(mos(2), qs_kind_set, particle_set, dft_section, &
516 88 : 4, 0, final_mos=.FALSE., spin="XAS BETA")
517 : ELSE
518 : CALL write_mo_set_to_output_unit(mos(1), qs_kind_set, particle_set, dft_section, &
519 0 : 4, 0, final_mos=.FALSE., spin="XAS")
520 : END IF
521 :
522 : ELSE
523 : ! Core level spectroscopy by TDDFT is not yet implemented
524 : ! the states defined by the rotation are the ground state orbitals
525 : ! the initial state from which I excite should be localized
526 : ! I take the excitations from lumo to nmo
527 : END IF
528 :
529 88 : IF (converged) THEN
530 : CALL cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, &
531 54 : iatom, istate)
532 : ELSE
533 34 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,/,T10,A,I6)') &
534 17 : "SCF with core hole NOT converged for ATOM ", iatom
535 : END IF
536 :
537 258 : IF (.NOT. skip_scf) THEN
538 : ! Reset the initial core orbitals.
539 : ! The valence orbitals are taken from the last SCF,
540 : ! it should be a better initial guess
541 80 : CALL get_qs_env(qs_env, mos=mos)
542 240 : DO ispin = 1, nspins
543 160 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
544 240 : CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
545 : END DO
546 80 : IF (iat == nexc_atoms) THEN
547 44 : CALL scf_env_cleanup(scf_env)
548 44 : CALL scf_env_release(xas_env%scf_env)
549 44 : DEALLOCATE (xas_env%scf_env)
550 : END IF
551 : END IF
552 :
553 : END DO ! istate
554 : END DO ! iat = 1,nexc_atoms
555 :
556 : ! END of Calculation
557 :
558 : ! Release what has to be released
559 42 : IF (ASSOCIATED(vecbuffer)) THEN
560 42 : DEALLOCATE (vecbuffer)
561 42 : DEALLOCATE (op_sm)
562 : END IF
563 :
564 126 : DO ispin = 1, dft_control%nspins
565 : CALL set_mo_set(mos(ispin), homo=my_homo(ispin), &
566 84 : uniform_occupation=my_uocc(ispin), nelectron=my_nelectron(ispin))
567 84 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
568 126 : CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
569 : END DO
570 :
571 42 : IF (output_unit > 0) THEN
572 : WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T25,A,/,T3,A,/)") &
573 21 : REPEAT("=", 77), &
574 21 : "END CORE LEVEL SPECTROSCOPY CALCULATION", &
575 42 : REPEAT("=", 77)
576 : END IF
577 :
578 42 : CALL xas_env_release(qs_env%xas_env)
579 42 : DEALLOCATE (qs_env%xas_env)
580 42 : NULLIFY (qs_env%xas_env)
581 :
582 : CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
583 42 : "PRINT%PROGRAM_RUN_INFO")
584 42 : CALL timestop(handle)
585 :
586 84 : END SUBROUTINE xas
587 :
588 : ! **************************************************************************************************
589 : !> \brief allocate and initialize the structure needed for the xas calculation
590 : !> \param xas_env the environment for XAS calculations
591 : !> \param qs_env the qs_env, the xas_env lives in
592 : !> \param dft_section ...
593 : !> \param logger ...
594 : !> \par History
595 : !> 05.2005 created [MI]
596 : !> \author MI
597 : ! **************************************************************************************************
598 42 : SUBROUTINE xas_env_init(xas_env, qs_env, dft_section, logger)
599 :
600 : TYPE(xas_environment_type), POINTER :: xas_env
601 : TYPE(qs_environment_type), POINTER :: qs_env
602 : TYPE(section_vals_type), POINTER :: dft_section
603 : TYPE(cp_logger_type), POINTER :: logger
604 :
605 : CHARACTER(LEN=default_string_length) :: name_sto
606 : INTEGER :: homo, i, iat, iatom, ik, ikind, ispin, j, l, lfomo, my_spin, n_mo(2), n_rep, nao, &
607 : natom, ncubes, nelectron, nexc_atoms, nexc_search, nj, nk, nkind, nmo, nmoloc(2), &
608 : nsgf_gto, nsgf_sto, nspins, nvirtual, nvirtual2
609 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_type_tmp, kind_z_tmp, &
610 42 : last_sgf
611 : INTEGER, DIMENSION(4, 7) :: ne
612 42 : INTEGER, DIMENSION(:), POINTER :: bounds, list, lq, nq, row_blk_sizes
613 : LOGICAL :: ihavethis
614 : REAL(dp) :: nele, occ_estate, occ_homo, &
615 : occ_homo_plus, zatom
616 42 : REAL(dp), DIMENSION(:), POINTER :: sto_zet
617 42 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
618 : TYPE(atomic_kind_type), POINTER :: atomic_kind
619 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
620 : TYPE(cp_fm_type), POINTER :: mo_coeff
621 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
622 42 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
623 : TYPE(dft_control_type), POINTER :: dft_control
624 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
625 42 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
626 : TYPE(mp_para_env_type), POINTER :: para_env
627 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
628 42 : POINTER :: sab_orb
629 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
630 42 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
631 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
632 : TYPE(qs_matrix_pools_type), POINTER :: mpools
633 : TYPE(scf_control_type), POINTER :: scf_control
634 : TYPE(section_vals_type), POINTER :: loc_section, xas_section
635 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
636 : TYPE(xas_control_type), POINTER :: xas_control
637 :
638 42 : n_mo(1:2) = 0
639 0 : CPASSERT(ASSOCIATED(xas_env))
640 :
641 42 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, scf_control, matrix_s, mos, mpools)
642 42 : NULLIFY (para_env, particle_set, xas_control)
643 42 : NULLIFY (qs_loc_env)
644 42 : NULLIFY (sab_orb)
645 : CALL get_qs_env(qs_env=qs_env, &
646 : atomic_kind_set=atomic_kind_set, &
647 : qs_kind_set=qs_kind_set, &
648 : dft_control=dft_control, &
649 : mpools=mpools, &
650 : matrix_s=matrix_s, mos=mos, &
651 : para_env=para_env, particle_set=particle_set, &
652 : sab_orb=sab_orb, &
653 42 : dbcsr_dist=dbcsr_dist)
654 :
655 42 : xas_section => section_vals_get_subs_vals(dft_section, "XAS")
656 42 : ALLOCATE (dft_control%xas_control)
657 42 : CALL xas_control_create(dft_control%xas_control)
658 42 : CALL read_xas_control(dft_control%xas_control, xas_section)
659 42 : CALL write_xas_control(dft_control%xas_control, dft_section)
660 42 : xas_control => dft_control%xas_control
661 1218 : ALLOCATE (scf_control)
662 42 : CALL scf_c_create(scf_control)
663 42 : CALL scf_c_read_parameters(scf_control, xas_section)
664 42 : CALL set_xas_env(xas_env, scf_control=scf_control)
665 :
666 42 : my_spin = xas_control%spin_channel
667 42 : nexc_search = xas_control%nexc_search
668 42 : IF (nexc_search < 0) THEN
669 : ! ground state occupation
670 2 : CALL get_mo_set(mos(my_spin), nmo=nmo, lfomo=lfomo)
671 2 : nexc_search = lfomo - 1
672 : END IF
673 42 : nexc_atoms = xas_control%nexc_atoms
674 126 : ALLOCATE (xas_env%exc_atoms(nexc_atoms))
675 206 : xas_env%exc_atoms = xas_control%exc_atoms
676 : CALL set_xas_env(xas_env=xas_env, nexc_search=nexc_search, &
677 42 : nexc_atoms=nexc_atoms, spin_channel=my_spin)
678 :
679 42 : CALL mpools_get(mpools, ao_mo_fm_pools=xas_env%ao_mo_fm_pools)
680 :
681 42 : NULLIFY (mo_coeff)
682 42 : CALL get_mo_set(mos(my_spin), nao=nao, homo=homo, nmo=nmo, mo_coeff=mo_coeff, nelectron=nelectron)
683 :
684 42 : nvirtual2 = 0
685 42 : IF (xas_control%added_mos > 0) THEN
686 40 : nvirtual2 = MIN(xas_control%added_mos, nao - nmo)
687 40 : xas_env%unoccupied_eps = xas_control%eps_added
688 40 : xas_env%unoccupied_max_iter = xas_control%max_iter_added
689 : END IF
690 42 : nvirtual = nmo + nvirtual2
691 :
692 126 : n_mo(1:2) = nmo
693 :
694 126 : ALLOCATE (xas_env%centers_wfn(3, nexc_search))
695 126 : ALLOCATE (xas_env%atom_of_state(nexc_search))
696 84 : ALLOCATE (xas_env%type_of_state(nexc_search))
697 168 : ALLOCATE (xas_env%state_of_atom(nexc_atoms, nexc_search))
698 84 : ALLOCATE (xas_env%nexc_states(nexc_atoms))
699 84 : ALLOCATE (xas_env%mykind_of_atom(nexc_atoms))
700 42 : nkind = SIZE(atomic_kind_set, 1)
701 126 : ALLOCATE (xas_env%mykind_of_kind(nkind))
702 124 : xas_env%mykind_of_kind = 0
703 :
704 : ! create a new matrix structure nao x 1
705 42 : NULLIFY (tmp_fm_struct)
706 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
707 42 : ncol_global=1, para_env=para_env, context=mo_coeff%matrix_struct%context)
708 42 : ALLOCATE (xas_env%excvec_coeff)
709 42 : CALL cp_fm_create(xas_env%excvec_coeff, tmp_fm_struct)
710 42 : CALL cp_fm_struct_release(tmp_fm_struct)
711 :
712 42 : NULLIFY (tmp_fm_struct)
713 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=1, &
714 : ncol_global=nexc_search, para_env=para_env, &
715 42 : context=mo_coeff%matrix_struct%context)
716 42 : ALLOCATE (xas_env%excvec_overlap)
717 42 : CALL cp_fm_create(xas_env%excvec_overlap, tmp_fm_struct)
718 42 : CALL cp_fm_struct_release(tmp_fm_struct)
719 :
720 42 : nspins = SIZE(mos, 1)
721 :
722 : ! initialize operators for the calculation of the oscillator strengths
723 42 : IF (xas_control%xas_method == xas_tp_hh) THEN
724 8 : occ_estate = 0.5_dp
725 8 : nele = REAL(nelectron, dp) - 0.5_dp
726 8 : occ_homo = 1.0_dp
727 8 : occ_homo_plus = 0._dp
728 : ELSE IF (xas_control%xas_method == xas_tp_xhh) THEN
729 4 : occ_estate = 0.5_dp
730 4 : nele = REAL(nelectron, dp)
731 4 : occ_homo = 1.0_dp
732 4 : occ_homo_plus = 0.5_dp
733 : ELSE IF (xas_control%xas_method == xas_tp_fh) THEN
734 10 : occ_estate = 0.0_dp
735 10 : nele = REAL(nelectron, dp) - 1.0_dp
736 10 : occ_homo = 1.0_dp
737 10 : occ_homo_plus = 0._dp
738 : ELSE IF (xas_control%xas_method == xas_tp_xfh) THEN
739 8 : occ_estate = 0.0_dp
740 8 : nele = REAL(nelectron, dp)
741 8 : occ_homo = 1.0_dp
742 8 : occ_homo_plus = 1._dp
743 : ELSE IF (xas_control%xas_method == xes_tp_val) THEN
744 6 : occ_estate = xas_control%xes_core_occupation
745 6 : nele = REAL(nelectron, dp) - xas_control%xes_core_occupation
746 6 : occ_homo = xas_control%xes_homo_occupation
747 : ELSE IF (xas_control%xas_method == xas_dscf) THEN
748 0 : occ_estate = 0.0_dp
749 0 : nele = REAL(nelectron, dp)
750 0 : occ_homo = 1.0_dp
751 0 : occ_homo_plus = 1._dp
752 : ELSE IF (xas_control%xas_method == xas_tp_flex) THEN
753 6 : nele = REAL(xas_control%nel_tot, dp)
754 6 : occ_estate = REAL(xas_control%xas_core_occupation, dp)
755 6 : IF (nele < 0.0_dp) nele = REAL(nelectron, dp) - (1.0_dp - occ_estate)
756 6 : occ_homo = 1.0_dp
757 : END IF
758 : CALL set_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_nelectron=nele, &
759 42 : nvirtual2=nvirtual2, nvirtual=nvirtual, homo_occ=occ_homo)
760 :
761 : ! Initialize the list of orbitals for cube files printing
762 42 : IF (BTEST(cp_print_key_should_output(logger%iter_info, xas_section, &
763 : "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
764 2 : NULLIFY (bounds, list)
765 : CALL section_vals_val_get(xas_section, &
766 : "PRINT%CLS_FUNCTION_CUBES%CUBES_LU_BOUNDS", &
767 2 : i_vals=bounds)
768 2 : ncubes = bounds(2) - bounds(1) + 1
769 2 : IF (ncubes > 0) THEN
770 0 : ALLOCATE (xas_control%list_cubes(ncubes))
771 :
772 0 : DO ik = 1, ncubes
773 0 : xas_control%list_cubes(ik) = bounds(1) + (ik - 1)
774 : END DO
775 : END IF
776 :
777 2 : IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
778 : CALL section_vals_val_get(xas_section, &
779 : "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
780 2 : n_rep_val=n_rep)
781 2 : ncubes = 0
782 4 : DO ik = 1, n_rep
783 2 : NULLIFY (list)
784 : CALL section_vals_val_get(xas_section, &
785 : "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
786 2 : i_rep_val=ik, i_vals=list)
787 4 : IF (ASSOCIATED(list)) THEN
788 2 : CALL reallocate(xas_control%list_cubes, 1, ncubes + SIZE(list))
789 8 : DO i = 1, SIZE(list)
790 8 : xas_control%list_cubes(i + ncubes) = list(i)
791 : END DO
792 2 : ncubes = ncubes + SIZE(list)
793 : END IF
794 : END DO ! ik
795 : END IF
796 :
797 2 : IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
798 0 : ncubes = MAX(10, xas_control%added_mos/10)
799 0 : ncubes = MIN(ncubes, xas_control%added_mos)
800 0 : ALLOCATE (xas_control%list_cubes(ncubes))
801 0 : DO ik = 1, ncubes
802 0 : xas_control%list_cubes(ik) = homo + ik
803 : END DO
804 : END IF
805 : ELSE
806 40 : NULLIFY (xas_control%list_cubes)
807 : END IF
808 :
809 42 : NULLIFY (tmp_fm_struct)
810 210 : ALLOCATE (xas_env%groundstate_coeff(nspins))
811 126 : DO ispin = 1, nspins
812 84 : CALL get_mo_set(mos(ispin), nao=nao, nmo=nmo)
813 : CALL fm_pool_create_fm(xas_env%ao_mo_fm_pools(ispin)%pool, &
814 : xas_env%groundstate_coeff(ispin), &
815 126 : name="xas_env%mo0"//TRIM(ADJUSTL(cp_to_string(ispin))))
816 : END DO ! ispin
817 :
818 42 : NULLIFY (tmp_fm_struct)
819 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=1, &
820 : ncol_global=nvirtual, para_env=para_env, &
821 42 : context=mo_coeff%matrix_struct%context)
822 420 : ALLOCATE (xas_env%dip_fm_set(2, 3))
823 168 : DO i = 1, 3
824 420 : DO j = 1, 2
825 378 : CALL cp_fm_create(xas_env%dip_fm_set(j, i), tmp_fm_struct)
826 : END DO
827 : END DO
828 42 : CALL cp_fm_struct_release(tmp_fm_struct)
829 :
830 : !Array to store all the eigenstates: occupied and the required not occupied
831 42 : IF (nvirtual2 > 0) THEN
832 120 : ALLOCATE (xas_env%unoccupied_evals(nvirtual2))
833 40 : NULLIFY (tmp_fm_struct)
834 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
835 : ncol_global=nvirtual2, &
836 40 : para_env=para_env, context=mo_coeff%matrix_struct%context)
837 40 : ALLOCATE (xas_env%unoccupied_orbs)
838 40 : CALL cp_fm_create(xas_env%unoccupied_orbs, tmp_fm_struct)
839 40 : CALL cp_fm_struct_release(tmp_fm_struct)
840 : END IF
841 :
842 42 : NULLIFY (tmp_fm_struct)
843 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
844 : ncol_global=nvirtual, &
845 42 : para_env=para_env, context=mo_coeff%matrix_struct%context)
846 42 : ALLOCATE (xas_env%all_vectors)
847 42 : CALL cp_fm_create(xas_env%all_vectors, tmp_fm_struct)
848 42 : CALL cp_fm_struct_release(tmp_fm_struct)
849 :
850 : ! Array to store all the energies needed for the spectrum
851 126 : ALLOCATE (xas_env%all_evals(nvirtual))
852 :
853 42 : IF (xas_control%dipole_form == xas_dip_len) THEN
854 0 : CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
855 0 : DO i = 1, 3
856 0 : ALLOCATE (xas_env%ostrength_sm(i)%matrix)
857 : CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, matrix_s(1)%matrix, &
858 0 : "xas_env%ostrength_sm-"//TRIM(ADJUSTL(cp_to_string(i))))
859 0 : CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
860 : END DO
861 42 : ELSE IF (xas_control%dipole_form == xas_dip_vel) THEN
862 : !
863 : ! prepare for allocation
864 42 : natom = SIZE(particle_set, 1)
865 126 : ALLOCATE (first_sgf(natom))
866 84 : ALLOCATE (last_sgf(natom))
867 : CALL get_particle_set(particle_set, qs_kind_set, &
868 : first_sgf=first_sgf, &
869 42 : last_sgf=last_sgf)
870 84 : ALLOCATE (row_blk_sizes(natom))
871 42 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
872 42 : DEALLOCATE (first_sgf)
873 42 : DEALLOCATE (last_sgf)
874 : !
875 : !
876 42 : CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
877 42 : ALLOCATE (xas_env%ostrength_sm(1)%matrix)
878 : CALL dbcsr_create(matrix=xas_env%ostrength_sm(1)%matrix, &
879 : name="xas_env%ostrength_sm", &
880 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
881 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
882 42 : mutable_work=.TRUE.)
883 42 : CALL cp_dbcsr_alloc_block_from_nbl(xas_env%ostrength_sm(1)%matrix, sab_orb)
884 42 : CALL dbcsr_set(xas_env%ostrength_sm(1)%matrix, 0.0_dp)
885 126 : DO i = 2, 3
886 84 : ALLOCATE (xas_env%ostrength_sm(i)%matrix)
887 : CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, xas_env%ostrength_sm(1)%matrix, &
888 84 : "xas_env%ostrength_sm-"//TRIM(ADJUSTL(cp_to_string(i))))
889 126 : CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
890 : END DO
891 :
892 42 : DEALLOCATE (row_blk_sizes)
893 : END IF
894 :
895 : ! Define the qs_loc_env : to find centers, spread and possibly localize them
896 42 : IF (.NOT. (ASSOCIATED(xas_env%qs_loc_env))) THEN
897 294 : ALLOCATE (qs_loc_env)
898 42 : CALL qs_loc_env_create(qs_loc_env)
899 42 : CALL set_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env)
900 42 : loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE")
901 :
902 : CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.TRUE., &
903 42 : do_xas=.TRUE., nloc_xas=nexc_search, spin_xas=my_spin)
904 :
905 42 : IF (.NOT. qs_loc_env%do_localize) THEN
906 0 : qs_loc_env%localized_wfn_control%localization_method = do_loc_none
907 :
908 : ELSE
909 126 : nmoloc = qs_loc_env%localized_wfn_control%nloc_states
910 42 : CALL set_loc_wfn_lists(qs_loc_env%localized_wfn_control, nmoloc, n_mo, nspins, my_spin)
911 42 : CALL set_loc_centers(qs_loc_env%localized_wfn_control, nmoloc, nspins)
912 : CALL qs_loc_env_init(qs_loc_env, qs_loc_env%localized_wfn_control, &
913 42 : qs_env, myspin=my_spin, do_localize=qs_loc_env%do_localize)
914 : END IF
915 : END IF
916 :
917 : !Type of state
918 42 : ALLOCATE (nq(1), lq(1), sto_zet(1))
919 42 : IF (xas_control%state_type == xas_1s_type) THEN
920 40 : nq(1) = 1
921 40 : lq(1) = 0
922 : ELSE IF (xas_control%state_type == xas_2s_type) THEN
923 0 : nq(1) = 2
924 0 : lq(1) = 0
925 : ELSE IF (xas_control%state_type == xas_2p_type) THEN
926 2 : nq(1) = 2
927 2 : lq(1) = 1
928 : ELSE IF (xas_control%state_type == xas_3s_type) THEN
929 0 : nq(1) = 3
930 0 : lq(1) = 0
931 : ELSE IF (xas_control%state_type == xas_3p_type) THEN
932 0 : nq(1) = 3
933 0 : lq(1) = 1
934 : ELSE IF (xas_control%state_type == xas_3d_type) THEN
935 0 : nq(1) = 3
936 0 : lq(1) = 2
937 : ELSE IF (xas_control%state_type == xas_4s_type) THEN
938 0 : nq(1) = 4
939 0 : lq(1) = 0
940 : ELSE IF (xas_control%state_type == xas_4p_type) THEN
941 0 : nq(1) = 4
942 0 : lq(1) = 1
943 : ELSE IF (xas_control%state_type == xas_4d_type) THEN
944 0 : nq(1) = 4
945 0 : lq(1) = 2
946 : ELSE IF (xas_control%state_type == xas_4f_type) THEN
947 0 : nq(1) = 4
948 0 : lq(1) = 3
949 : ELSE
950 0 : CPABORT("XAS type of state not implemented")
951 : END IF
952 :
953 : ! Find core orbitals of right angular momentum
954 84 : ALLOCATE (kind_type_tmp(nkind))
955 84 : ALLOCATE (kind_z_tmp(nkind))
956 42 : kind_type_tmp = 0
957 42 : kind_z_tmp = 0
958 42 : nk = 0
959 124 : DO iat = 1, nexc_atoms
960 82 : iatom = xas_env%exc_atoms(iat)
961 82 : NULLIFY (atomic_kind)
962 82 : atomic_kind => particle_set(iatom)%atomic_kind
963 82 : CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=ikind)
964 82 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zatom)
965 82 : ihavethis = .FALSE.
966 112 : DO ik = 1, nk
967 112 : IF (ikind == kind_type_tmp(ik)) THEN
968 10 : ihavethis = .TRUE.
969 10 : xas_env%mykind_of_atom(iat) = ik
970 : EXIT
971 : END IF
972 : END DO
973 124 : IF (.NOT. ihavethis) THEN
974 72 : nk = nk + 1
975 72 : kind_type_tmp(nk) = ikind
976 72 : kind_z_tmp(nk) = INT(zatom)
977 72 : xas_env%mykind_of_atom(iat) = nk
978 72 : xas_env%mykind_of_kind(ikind) = nk
979 : END IF
980 : END DO ! iat
981 :
982 198 : ALLOCATE (xas_env%my_gto_basis(nk))
983 198 : ALLOCATE (xas_env%stogto_overlap(nk))
984 114 : DO ik = 1, nk
985 72 : NULLIFY (xas_env%my_gto_basis(ik)%gto_basis_set, sto_basis_set)
986 72 : ne = 0
987 146 : DO l = 1, lq(1) + 1
988 74 : nj = 2*(l - 1) + 1
989 222 : DO i = l, nq(1)
990 76 : ne(l, i) = ptable(kind_z_tmp(ik))%e_conv(l - 1) - 2*nj*(i - l)
991 76 : ne(l, i) = MAX(ne(l, i), 0)
992 150 : ne(l, i) = MIN(ne(l, i), 2*nj)
993 : END DO
994 : END DO
995 :
996 72 : sto_zet(1) = srules(kind_z_tmp(ik), ne, nq(1), lq(1))
997 72 : CALL allocate_sto_basis_set(sto_basis_set)
998 72 : name_sto = 'xas_tmp_sto'
999 : CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, &
1000 72 : lq=lq, zet=sto_zet, name=name_sto)
1001 : CALL create_gto_from_sto_basis(sto_basis_set, &
1002 72 : xas_env%my_gto_basis(ik)%gto_basis_set, xas_control%ngauss)
1003 72 : CALL deallocate_sto_basis_set(sto_basis_set)
1004 72 : xas_env%my_gto_basis(ik)%gto_basis_set%norm_type = 2
1005 72 : CALL init_orb_basis_set(xas_env%my_gto_basis(ik)%gto_basis_set)
1006 :
1007 72 : ikind = kind_type_tmp(ik)
1008 72 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1009 :
1010 72 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf_gto)
1011 72 : CALL get_gto_basis_set(gto_basis_set=xas_env%my_gto_basis(ik)%gto_basis_set, nsgf=nsgf_sto)
1012 288 : ALLOCATE (xas_env%stogto_overlap(ik)%array(nsgf_sto, nsgf_gto))
1013 :
1014 : CALL calc_stogto_overlap(xas_env%my_gto_basis(ik)%gto_basis_set, orb_basis_set, &
1015 186 : xas_env%stogto_overlap(ik)%array)
1016 : END DO
1017 :
1018 42 : DEALLOCATE (nq, lq, sto_zet)
1019 42 : DEALLOCATE (kind_type_tmp, kind_z_tmp)
1020 :
1021 126 : END SUBROUTINE xas_env_init
1022 :
1023 : ! **************************************************************************************************
1024 : !> \brief Calculate and write the spectrum relative to the core level excitation
1025 : !> of a specific atom. It works for TP approach, because of the definition
1026 : !> of the oscillator strengths as matrix elements of the dipole operator
1027 : !> \param xas_control ...
1028 : !> \param xas_env ...
1029 : !> \param qs_env ...
1030 : !> \param xas_section ...
1031 : !> \param iatom index of the excited atom
1032 : !> \param istate ...
1033 : !> \par History
1034 : !> 03.2006 created [MI]
1035 : !> \author MI
1036 : !> \note
1037 : !> for the tddft calculation should be re-thought
1038 : ! **************************************************************************************************
1039 54 : SUBROUTINE cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, &
1040 : iatom, istate)
1041 :
1042 : TYPE(xas_control_type) :: xas_control
1043 : TYPE(xas_environment_type), POINTER :: xas_env
1044 : TYPE(qs_environment_type), POINTER :: qs_env
1045 : TYPE(section_vals_type), POINTER :: xas_section
1046 : INTEGER, INTENT(IN) :: iatom, istate
1047 :
1048 : INTEGER :: homo, i, lfomo, my_spin, nabs, nmo, &
1049 : nvirtual, output_unit, xas_estate
1050 : LOGICAL :: append_cube, length
1051 : REAL(dp) :: rc(3)
1052 54 : REAL(dp), DIMENSION(:), POINTER :: all_evals
1053 : REAL(dp), DIMENSION(:, :), POINTER :: sp_ab, sp_em
1054 54 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dip_fm_set
1055 : TYPE(cp_fm_type), POINTER :: all_vectors, excvec_coeff
1056 : TYPE(cp_logger_type), POINTER :: logger
1057 54 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_sm, ostrength_sm
1058 54 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1059 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1060 :
1061 54 : NULLIFY (logger)
1062 108 : logger => cp_get_default_logger()
1063 54 : output_unit = cp_logger_get_default_io_unit(logger)
1064 :
1065 54 : NULLIFY (ostrength_sm, op_sm, dip_fm_set)
1066 54 : NULLIFY (all_evals, all_vectors, excvec_coeff)
1067 54 : NULLIFY (mos, particle_set, sp_em, sp_ab)
1068 216 : ALLOCATE (op_sm(3))
1069 :
1070 : CALL get_qs_env(qs_env=qs_env, &
1071 54 : mos=mos, particle_set=particle_set)
1072 :
1073 : CALL get_xas_env(xas_env=xas_env, all_vectors=all_vectors, xas_estate=xas_estate, &
1074 : all_evals=all_evals, dip_fm_set=dip_fm_set, excvec_coeff=excvec_coeff, &
1075 54 : ostrength_sm=ostrength_sm, nvirtual=nvirtual, spin_channel=my_spin)
1076 54 : CALL get_mo_set(mos(my_spin), homo=homo, lfomo=lfomo, nmo=nmo)
1077 :
1078 54 : nabs = nvirtual - lfomo + 1
1079 162 : ALLOCATE (sp_em(6, homo))
1080 162 : ALLOCATE (sp_ab(6, nabs))
1081 54 : CPASSERT(ASSOCIATED(excvec_coeff))
1082 :
1083 54 : IF (.NOT. xas_control%xas_method == xas_dscf) THEN
1084 : ! Calculate the spectrum
1085 54 : IF (xas_control%dipole_form == xas_dip_len) THEN
1086 0 : rc(1:3) = particle_set(iatom)%r(1:3)
1087 0 : DO i = 1, 3
1088 0 : NULLIFY (op_sm(i)%matrix)
1089 0 : op_sm(i)%matrix => ostrength_sm(i)%matrix
1090 : END DO
1091 0 : CALL rRc_xyz_ao(op_sm, qs_env, rc, order=1, minimum_image=.TRUE.)
1092 : CALL spectrum_dip_vel(dip_fm_set, op_sm, mos, excvec_coeff, &
1093 : all_vectors, all_evals, &
1094 0 : sp_em, sp_ab, xas_estate, nvirtual, my_spin)
1095 0 : DO i = 1, SIZE(ostrength_sm, 1)
1096 0 : CALL dbcsr_set(ostrength_sm(i)%matrix, 0.0_dp)
1097 : END DO
1098 : ELSE
1099 216 : DO i = 1, 3
1100 162 : NULLIFY (op_sm(i)%matrix)
1101 216 : op_sm(i)%matrix => ostrength_sm(i)%matrix
1102 : END DO
1103 : CALL spectrum_dip_vel(dip_fm_set, op_sm, mos, excvec_coeff, &
1104 : all_vectors, all_evals, &
1105 54 : sp_em, sp_ab, xas_estate, nvirtual, my_spin)
1106 : END IF
1107 : END IF
1108 :
1109 54 : CALL get_mo_set(mos(my_spin), lfomo=lfomo)
1110 : ! write the spectrum, if the file exists it is appended
1111 54 : IF (.NOT. xas_control%xas_method == xas_dscf) THEN
1112 54 : length = (.NOT. xas_control%dipole_form == xas_dip_vel)
1113 : CALL xas_write(sp_em, sp_ab, xas_estate, &
1114 54 : xas_section, iatom, istate, lfomo, length=length)
1115 : END IF
1116 :
1117 54 : DEALLOCATE (sp_em)
1118 54 : DEALLOCATE (sp_ab)
1119 :
1120 54 : IF (BTEST(cp_print_key_should_output(logger%iter_info, xas_section, &
1121 : "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
1122 4 : append_cube = section_get_lval(xas_section, "PRINT%CLS_FUNCTION_CUBES%APPEND")
1123 : CALL xas_print_cubes(xas_control, qs_env, xas_section, mos, all_vectors, &
1124 4 : iatom, append_cube)
1125 : END IF
1126 :
1127 54 : IF (BTEST(cp_print_key_should_output(logger%iter_info, xas_section, &
1128 : "PRINT%PDOS"), cp_p_file)) THEN
1129 4 : CALL xas_pdos(qs_env, xas_section, mos, iatom)
1130 : END IF
1131 :
1132 54 : DEALLOCATE (op_sm)
1133 :
1134 162 : END SUBROUTINE cls_calculate_spectrum
1135 :
1136 : ! **************************************************************************************************
1137 : !> \brief write the spectrum for each atom in a different output file
1138 : !> \param sp_em ...
1139 : !> \param sp_ab ...
1140 : !> \param estate ...
1141 : !> \param xas_section ...
1142 : !> \param iatom index of the excited atom
1143 : !> \param state_to_be_excited ...
1144 : !> \param lfomo ...
1145 : !> \param length ...
1146 : !> \par History
1147 : !> 05.2005 created [MI]
1148 : !> \author MI
1149 : !> \note
1150 : !> the iteration counter is not finilized yet
1151 : ! **************************************************************************************************
1152 54 : SUBROUTINE xas_write(sp_em, sp_ab, estate, xas_section, iatom, state_to_be_excited, &
1153 : lfomo, length)
1154 :
1155 : REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab
1156 : INTEGER, INTENT(IN) :: estate
1157 : TYPE(section_vals_type), POINTER :: xas_section
1158 : INTEGER, INTENT(IN) :: iatom, state_to_be_excited, lfomo
1159 : LOGICAL, INTENT(IN) :: length
1160 :
1161 : CHARACTER(LEN=default_string_length) :: mittle_ab, mittle_em, my_act, my_pos
1162 : INTEGER :: i, istate, out_sp_ab, out_sp_em
1163 : REAL(dp) :: ene2
1164 : TYPE(cp_logger_type), POINTER :: logger
1165 :
1166 54 : NULLIFY (logger)
1167 54 : logger => cp_get_default_logger()
1168 :
1169 54 : my_pos = "APPEND"
1170 54 : my_act = "WRITE"
1171 :
1172 54 : mittle_em = "xes_at"//TRIM(ADJUSTL(cp_to_string(iatom)))//"_st"//TRIM(ADJUSTL(cp_to_string(state_to_be_excited)))
1173 :
1174 : out_sp_em = cp_print_key_unit_nr(logger, xas_section, "PRINT%XES_SPECTRUM", &
1175 : extension=".spectrum", file_position=my_pos, file_action=my_act, &
1176 54 : file_form="FORMATTED", middle_name=TRIM(mittle_em))
1177 :
1178 54 : IF (out_sp_em > 0) THEN
1179 27 : WRITE (out_sp_em, '(A,I6,A,I6,A,I6)') " Emission spectrum for atom ", iatom, &
1180 54 : ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_em, 2)
1181 27 : ene2 = 1.0_dp
1182 318 : DO istate = estate, SIZE(sp_em, 2)
1183 291 : IF (length) ene2 = sp_em(1, istate)*sp_em(1, istate)
1184 291 : WRITE (out_sp_em, '(I6,5F16.8,F10.5)') istate, sp_em(1, istate)*evolt, &
1185 291 : sp_em(2, istate)*ene2, sp_em(3, istate)*ene2, &
1186 609 : sp_em(4, istate)*ene2, sp_em(5, istate)*ene2, sp_em(6, istate)
1187 : END DO
1188 : END IF
1189 : CALL cp_print_key_finished_output(out_sp_em, logger, xas_section, &
1190 54 : "PRINT%XES_SPECTRUM")
1191 :
1192 54 : mittle_ab = "xas_at"//TRIM(ADJUSTL(cp_to_string(iatom)))//"_st"//TRIM(ADJUSTL(cp_to_string(state_to_be_excited)))
1193 : out_sp_ab = cp_print_key_unit_nr(logger, xas_section, "PRINT%XAS_SPECTRUM", &
1194 : extension=".spectrum", file_position=my_pos, file_action=my_act, &
1195 54 : file_form="FORMATTED", middle_name=TRIM(mittle_ab))
1196 :
1197 54 : IF (out_sp_ab > 0) THEN
1198 21 : WRITE (out_sp_ab, '(A,I6,A,I6,A,I6)') " Absorption spectrum for atom ", iatom, &
1199 42 : ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_ab, 2)
1200 21 : ene2 = 1.0_dp
1201 852 : DO i = 1, SIZE(sp_ab, 2)
1202 831 : istate = lfomo - 1 + i
1203 831 : IF (length) ene2 = sp_ab(1, i)*sp_ab(1, i)
1204 831 : WRITE (out_sp_ab, '(I6,5F16.8,F10.5)') istate, sp_ab(1, i)*evolt, &
1205 831 : sp_ab(2, i)*ene2, sp_ab(3, i)*ene2, &
1206 1683 : sp_ab(4, i)*ene2, sp_ab(5, i)*ene2, sp_ab(6, i)
1207 : END DO
1208 : END IF
1209 :
1210 : CALL cp_print_key_finished_output(out_sp_ab, logger, xas_section, &
1211 54 : "PRINT%XAS_SPECTRUM")
1212 :
1213 54 : END SUBROUTINE xas_write
1214 :
1215 : ! **************************************************************************************************
1216 : !> \brief write the cube files for a set of selected states
1217 : !> \param xas_control provide number ant indexes of the states to be printed
1218 : !> \param qs_env ...
1219 : !> \param xas_section ...
1220 : !> \param mos mos from which the states to be printed are extracted
1221 : !> \param all_vectors ...
1222 : !> \param iatom index of the atom that has been excited
1223 : !> \param append_cube ...
1224 : !> \par History
1225 : !> 08.2005 created [MI]
1226 : !> \author MI
1227 : ! **************************************************************************************************
1228 8 : SUBROUTINE xas_print_cubes(xas_control, qs_env, xas_section, &
1229 4 : mos, all_vectors, iatom, append_cube)
1230 :
1231 : TYPE(xas_control_type) :: xas_control
1232 : TYPE(qs_environment_type), POINTER :: qs_env
1233 : TYPE(section_vals_type), POINTER :: xas_section
1234 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1235 : TYPE(cp_fm_type), INTENT(IN) :: all_vectors
1236 : INTEGER, INTENT(IN) :: iatom
1237 : LOGICAL, INTENT(IN) :: append_cube
1238 :
1239 : CHARACTER(LEN=default_string_length) :: my_mittle, my_pos
1240 : INTEGER :: homo, istate0, my_spin, nspins, nstates
1241 4 : REAL(dp), DIMENSION(:, :), POINTER :: centers
1242 : TYPE(section_vals_type), POINTER :: print_key
1243 :
1244 4 : nspins = SIZE(mos)
1245 :
1246 8 : print_key => section_vals_get_subs_vals(xas_section, "PRINT%CLS_FUNCTION_CUBES")
1247 4 : my_mittle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom)))
1248 4 : nstates = SIZE(xas_control%list_cubes, 1)
1249 :
1250 4 : IF (xas_control%do_centers) THEN
1251 : ! one might like to calculate the centers of the xas orbital (without localizing them)
1252 : ELSE
1253 12 : ALLOCATE (centers(6, nstates))
1254 88 : centers = 0.0_dp
1255 : END IF
1256 4 : my_spin = xas_control%spin_channel
1257 :
1258 4 : CALL get_mo_set(mos(my_spin), homo=homo)
1259 4 : istate0 = 0
1260 :
1261 4 : my_pos = "REWIND"
1262 4 : IF (append_cube) THEN
1263 0 : my_pos = "APPEND"
1264 : END IF
1265 :
1266 : CALL qs_print_cubes(qs_env, all_vectors, nstates, xas_control%list_cubes, &
1267 4 : centers, print_key, my_mittle, state0=istate0, file_position=my_pos)
1268 :
1269 4 : DEALLOCATE (centers)
1270 :
1271 4 : END SUBROUTINE xas_print_cubes
1272 :
1273 : ! **************************************************************************************************
1274 : !> \brief write the PDOS after the XAS SCF, i.e., with one excited core
1275 : !> \param qs_env ...
1276 : !> \param xas_section ...
1277 : !> \param mos mos from which the eigenvalues and expansion coeffiecients are obtained
1278 : !> \param iatom index of the atom that has been excited
1279 : !> \par History
1280 : !> 03.2016 created [MI]
1281 : !> \author MI
1282 : ! **************************************************************************************************
1283 :
1284 4 : SUBROUTINE xas_pdos(qs_env, xas_section, mos, iatom)
1285 :
1286 : TYPE(qs_environment_type), POINTER :: qs_env
1287 : TYPE(section_vals_type), POINTER :: xas_section
1288 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1289 : INTEGER, INTENT(IN) :: iatom
1290 :
1291 : CHARACTER(LEN=default_string_length) :: xas_mittle
1292 : INTEGER :: ispin
1293 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1294 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1295 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1296 :
1297 4 : NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
1298 4 : xas_mittle = 'xasat'//TRIM(ADJUSTL(cp_to_string(iatom)))//'_'
1299 :
1300 : CALL get_qs_env(qs_env, &
1301 : atomic_kind_set=atomic_kind_set, &
1302 : particle_set=particle_set, &
1303 4 : qs_kind_set=qs_kind_set)
1304 :
1305 12 : DO ispin = 1, 2
1306 : CALL calculate_projected_dos(mos(ispin), atomic_kind_set, qs_kind_set, particle_set, qs_env, &
1307 12 : xas_section, ispin, xas_mittle)
1308 : END DO
1309 :
1310 4 : END SUBROUTINE xas_pdos
1311 : ! **************************************************************************************************
1312 : !> \brief Calculation of the spectrum when the dipole approximation
1313 : !> in the velocity form is used.
1314 : !> \param fm_set components of the position operator in a full matrix form
1315 : !> already multiplied by the coefficiets
1316 : !> only the terms <C_i Op C_f> are calculated where
1317 : !> C_i are the coefficients of the excited state
1318 : !> \param op_sm components of the position operator for the dipole
1319 : !> in a sparse matrix form (cos and sin)
1320 : !> calculated for the basis functions
1321 : !> \param mos wavefunctions coefficients
1322 : !> \param excvec coefficients of the excited orbital
1323 : !> \param all_vectors ...
1324 : !> \param all_evals ...
1325 : !> \param sp_em ...
1326 : !> \param sp_ab ...
1327 : !> \param estate index of the excited state
1328 : !> \param nstate ...
1329 : !> \param my_spin ...
1330 : !> \par History
1331 : !> 06.2005 created [MI]
1332 : !> \author MI
1333 : ! **************************************************************************************************
1334 54 : SUBROUTINE spectrum_dip_vel(fm_set, op_sm, mos, excvec, &
1335 : all_vectors, all_evals, sp_em, sp_ab, estate, nstate, my_spin)
1336 :
1337 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: fm_set
1338 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_sm
1339 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1340 : TYPE(cp_fm_type), INTENT(IN) :: excvec, all_vectors
1341 : REAL(dp), DIMENSION(:), POINTER :: all_evals
1342 : REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab
1343 : INTEGER, INTENT(IN) :: estate, nstate, my_spin
1344 :
1345 : INTEGER :: homo, i, i_abs, istate, lfomo, nao, nmo
1346 : REAL(dp) :: dip(3), ene_f, ene_i
1347 54 : REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
1348 : TYPE(cp_fm_type) :: fm_work
1349 :
1350 0 : CPASSERT(ASSOCIATED(fm_set))
1351 54 : NULLIFY (eigenvalues, occupation_numbers)
1352 :
1353 : CALL get_mo_set(mos(my_spin), eigenvalues=eigenvalues, occupation_numbers=occupation_numbers, &
1354 54 : nao=nao, nmo=nmo, homo=homo, lfomo=lfomo)
1355 :
1356 54 : CALL cp_fm_create(fm_work, all_vectors%matrix_struct)
1357 216 : DO i = 1, SIZE(fm_set, 2)
1358 162 : CALL cp_fm_set_all(fm_set(my_spin, i), 0.0_dp)
1359 162 : CALL cp_fm_set_all(fm_work, 0.0_dp)
1360 162 : CALL cp_dbcsr_sm_fm_multiply(op_sm(i)%matrix, all_vectors, fm_work, ncol=nstate)
1361 : CALL parallel_gemm("T", "N", 1, nstate, nao, 1.0_dp, excvec, &
1362 216 : fm_work, 0.0_dp, fm_set(my_spin, i), b_first_col=1)
1363 : END DO
1364 54 : CALL cp_fm_release(fm_work)
1365 :
1366 4282 : sp_em = 0.0_dp
1367 12360 : sp_ab = 0.0_dp
1368 54 : ene_i = eigenvalues(estate)
1369 2398 : DO istate = 1, nstate
1370 2344 : ene_f = all_evals(istate)
1371 9376 : DO i = 1, 3
1372 9376 : CALL cp_fm_get_element(fm_set(my_spin, i), 1, istate, dip(i))
1373 : END DO
1374 2344 : IF (istate <= homo) THEN
1375 604 : sp_em(1, istate) = ene_f - ene_i
1376 604 : sp_em(2, istate) = dip(1)
1377 604 : sp_em(3, istate) = dip(2)
1378 604 : sp_em(4, istate) = dip(3)
1379 604 : sp_em(5, istate) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
1380 604 : sp_em(6, istate) = occupation_numbers(istate)
1381 : END IF
1382 2398 : IF (istate >= lfomo) THEN
1383 1758 : i_abs = istate - lfomo + 1
1384 1758 : sp_ab(1, i_abs) = ene_f - ene_i
1385 1758 : sp_ab(2, i_abs) = dip(1)
1386 1758 : sp_ab(3, i_abs) = dip(2)
1387 1758 : sp_ab(4, i_abs) = dip(3)
1388 1758 : sp_ab(5, i_abs) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
1389 1758 : IF (istate <= nmo) sp_ab(6, i_abs) = occupation_numbers(istate)
1390 : END IF
1391 :
1392 : END DO
1393 :
1394 54 : END SUBROUTINE spectrum_dip_vel
1395 :
1396 : ! **************************************************************************************************
1397 : !> \brief ...
1398 : !> \param base_a ...
1399 : !> \param base_b ...
1400 : !> \param matrix ...
1401 : ! **************************************************************************************************
1402 162 : SUBROUTINE calc_stogto_overlap(base_a, base_b, matrix)
1403 :
1404 : TYPE(gto_basis_set_type), POINTER :: base_a, base_b
1405 : REAL(dp), DIMENSION(:, :), POINTER :: matrix
1406 :
1407 : INTEGER :: iset, jset, ldsab, maxcoa, maxcob, maxl, &
1408 : maxla, maxlb, na, nb, nseta, nsetb, &
1409 : nsgfa, nsgfb, sgfa, sgfb
1410 162 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1411 162 : npgfb, nsgfa_set, nsgfb_set
1412 162 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1413 : REAL(dp) :: rab(3)
1414 162 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
1415 162 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, sphi_a, &
1416 162 : sphi_b, zeta, zetb
1417 :
1418 162 : NULLIFY (la_max, la_min, lb_max, lb_min)
1419 162 : NULLIFY (npgfa, npgfb, nsgfa_set, nsgfb_set)
1420 162 : NULLIFY (first_sgfa, first_sgfb)
1421 162 : NULLIFY (rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
1422 :
1423 : CALL get_gto_basis_set(gto_basis_set=base_a, nsgf=nsgfa, nsgf_set=nsgfa_set, lmax=la_max, &
1424 : lmin=la_min, npgf=npgfa, pgf_radius=rpgfa, &
1425 : sphi=sphi_a, scon=scon_a, zet=zeta, first_sgf=first_sgfa, &
1426 162 : maxco=maxcoa, nset=nseta, maxl=maxla)
1427 :
1428 : CALL get_gto_basis_set(gto_basis_set=base_b, nsgf=nsgfb, nsgf_set=nsgfb_set, lmax=lb_max, &
1429 : lmin=lb_min, npgf=npgfb, pgf_radius=rpgfb, &
1430 : sphi=sphi_b, scon=scon_b, zet=zetb, first_sgf=first_sgfb, &
1431 162 : maxco=maxcob, nset=nsetb, maxl=maxlb)
1432 : ! Initialize and allocate
1433 162 : rab = 0.0_dp
1434 5886 : matrix = 0.0_dp
1435 :
1436 162 : ldsab = MAX(maxcoa, maxcob, nsgfa, nsgfb)
1437 162 : maxl = MAX(maxla, maxlb)
1438 :
1439 648 : ALLOCATE (sab(ldsab, ldsab))
1440 486 : ALLOCATE (work(ldsab, ldsab))
1441 :
1442 324 : DO iset = 1, nseta
1443 :
1444 162 : na = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1445 162 : sgfa = first_sgfa(1, iset)
1446 :
1447 1154 : DO jset = 1, nsetb
1448 830 : nb = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1449 830 : sgfb = first_sgfb(1, jset)
1450 :
1451 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1452 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1453 830 : rab, sab)
1454 : CALL contraction(sab, work, ca=scon_a(:, sgfa:), na=na, ma=nsgfa_set(iset), &
1455 830 : cb=scon_b(:, sgfb:), nb=nb, mb=nsgfb_set(jset))
1456 992 : CALL block_add("IN", work, nsgfa_set(iset), nsgfb_set(jset), matrix, sgfa, sgfb)
1457 :
1458 : END DO ! jset
1459 : END DO ! iset
1460 162 : DEALLOCATE (sab, work)
1461 :
1462 162 : END SUBROUTINE calc_stogto_overlap
1463 :
1464 : ! **************************************************************************************************
1465 : !> \brief Starting from a set of mos, determine on which atom are centered
1466 : !> and if they are of the right type (1s,2s ...)
1467 : !> to be used in the specific core level spectrum calculation
1468 : !> The set of states need to be from the core, otherwise the
1469 : !> characterization of the type is not valid, since it assumes that
1470 : !> the orbital is localizad on a specific atom
1471 : !> It is probably reccomandable to run a localization cycle before
1472 : !> proceeding to the assignment of the type
1473 : !> The type is determined by computing the overalp with a
1474 : !> type specific, minimal, STO bais set
1475 : !> \param xas_control ...
1476 : !> \param xas_env ...
1477 : !> \param localized_wfn_control ...
1478 : !> \param qs_env ...
1479 : !> \par History
1480 : !> 03.2006 created [MI]
1481 : !> \author MI
1482 : ! **************************************************************************************************
1483 42 : SUBROUTINE cls_assign_core_states(xas_control, xas_env, localized_wfn_control, qs_env)
1484 :
1485 : TYPE(xas_control_type) :: xas_control
1486 : TYPE(xas_environment_type), POINTER :: xas_env
1487 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1488 : TYPE(qs_environment_type), POINTER :: qs_env
1489 :
1490 : INTEGER :: chosen_state, homo, i, iat, iatom, &
1491 : ikind, isgf, istate, j, my_kind, &
1492 : my_spin, nao, natom, nexc_atoms, &
1493 : nexc_search, output_unit
1494 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf
1495 : INTEGER, DIMENSION(3) :: perd0
1496 42 : INTEGER, DIMENSION(:), POINTER :: atom_of_state, mykind_of_kind, &
1497 42 : nexc_states, state_of_mytype, &
1498 42 : type_of_state
1499 42 : INTEGER, DIMENSION(:, :), POINTER :: state_of_atom
1500 : REAL(dp) :: component, dist, distmin, maxocc, ra(3), &
1501 : rac(3), rc(3)
1502 42 : REAL(dp), DIMENSION(:), POINTER :: max_overlap, sto_state_overlap
1503 42 : REAL(dp), DIMENSION(:, :), POINTER :: centers_wfn
1504 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
1505 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1506 : TYPE(cell_type), POINTER :: cell
1507 42 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: stogto_overlap
1508 : TYPE(cp_fm_type), POINTER :: mo_coeff
1509 : TYPE(cp_logger_type), POINTER :: logger
1510 42 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1511 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1512 42 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1513 :
1514 42 : NULLIFY (cell, mos, particle_set)
1515 42 : NULLIFY (atom_of_state, centers_wfn, mykind_of_kind, state_of_atom, nexc_states)
1516 42 : NULLIFY (stogto_overlap, type_of_state, max_overlap, qs_kind_set)
1517 42 : NULLIFY (state_of_mytype, type_of_state, sto_state_overlap)
1518 :
1519 42 : NULLIFY (logger)
1520 84 : logger => cp_get_default_logger()
1521 42 : output_unit = cp_logger_get_default_io_unit(logger)
1522 :
1523 : CALL get_qs_env(qs_env=qs_env, cell=cell, mos=mos, particle_set=particle_set, &
1524 42 : qs_kind_set=qs_kind_set)
1525 :
1526 : ! The Berry operator can be used only for periodic systems
1527 : ! If an isolated system is used the periodicity is overimposed
1528 168 : perd0(1:3) = cell%perd(1:3)
1529 168 : cell%perd(1:3) = 1
1530 :
1531 : CALL get_xas_env(xas_env=xas_env, &
1532 : centers_wfn=centers_wfn, atom_of_state=atom_of_state, &
1533 : mykind_of_kind=mykind_of_kind, &
1534 : type_of_state=type_of_state, state_of_atom=state_of_atom, &
1535 : stogto_overlap=stogto_overlap, nexc_atoms=nexc_atoms, &
1536 42 : spin_channel=my_spin, nexc_search=nexc_search, nexc_states=nexc_states)
1537 :
1538 42 : CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff, maxocc=maxocc, nao=nao, homo=homo)
1539 :
1540 : ! scratch array for the state
1541 126 : ALLOCATE (vecbuffer(1, nao))
1542 42 : natom = SIZE(particle_set)
1543 :
1544 126 : ALLOCATE (first_sgf(natom))
1545 42 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
1546 126 : ALLOCATE (sto_state_overlap(nexc_search))
1547 126 : ALLOCATE (max_overlap(natom))
1548 166 : max_overlap = 0.0_dp
1549 84 : ALLOCATE (state_of_mytype(natom))
1550 166 : state_of_mytype = 0
1551 216 : atom_of_state = 0
1552 124 : nexc_states = 1
1553 592 : state_of_atom = 0
1554 :
1555 42 : IF (xas_control%orbital_list(1) < 0) THEN !Checks for manually selected orbitals from the localized set
1556 :
1557 198 : DO istate = 1, nexc_search
1558 158 : centers_wfn(1, istate) = localized_wfn_control%centers_set(my_spin)%array(1, istate)
1559 158 : centers_wfn(2, istate) = localized_wfn_control%centers_set(my_spin)%array(2, istate)
1560 158 : centers_wfn(3, istate) = localized_wfn_control%centers_set(my_spin)%array(3, istate)
1561 :
1562 : ! Assign the state to the closest atom
1563 158 : distmin = 100.0_dp
1564 518 : DO iat = 1, nexc_atoms
1565 360 : iatom = xas_control%exc_atoms(iat)
1566 1440 : ra(1:3) = particle_set(iatom)%r(1:3)
1567 1440 : rc(1:3) = centers_wfn(1:3, istate)
1568 360 : rac = pbc(ra, rc, cell)
1569 360 : dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
1570 :
1571 518 : IF (dist < distmin) THEN
1572 :
1573 238 : atom_of_state(istate) = iatom
1574 238 : distmin = dist
1575 : END IF
1576 : END DO
1577 198 : IF (atom_of_state(istate) /= 0) THEN
1578 : !Character of the state
1579 : CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, istate, &
1580 158 : nao, 1, transpose=.TRUE.)
1581 :
1582 158 : iatom = atom_of_state(istate)
1583 :
1584 158 : NULLIFY (atomic_kind)
1585 158 : atomic_kind => particle_set(iatom)%atomic_kind
1586 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1587 158 : kind_number=ikind)
1588 :
1589 158 : my_kind = mykind_of_kind(ikind)
1590 :
1591 158 : sto_state_overlap(istate) = 0.0_dp
1592 348 : DO i = 1, SIZE(stogto_overlap(my_kind)%array, 1)
1593 190 : component = 0.0_dp
1594 3214 : DO j = 1, SIZE(stogto_overlap(my_kind)%array, 2)
1595 3024 : isgf = first_sgf(iatom) + j - 1
1596 3214 : component = component + stogto_overlap(my_kind)%array(i, j)*vecbuffer(1, isgf)
1597 : END DO
1598 : sto_state_overlap(istate) = sto_state_overlap(istate) + &
1599 348 : component*component
1600 : END DO
1601 :
1602 158 : IF (sto_state_overlap(istate) > max_overlap(iatom)) THEN
1603 98 : state_of_mytype(iatom) = istate
1604 98 : max_overlap(iatom) = sto_state_overlap(istate)
1605 : END IF
1606 : END IF
1607 : END DO ! istate
1608 :
1609 : ! Includes all states within the chosen threshold relative to the maximum overlap
1610 40 : IF (xas_control%overlap_threshold < 1) THEN
1611 4 : DO iat = 1, nexc_atoms
1612 2 : iatom = xas_control%exc_atoms(iat)
1613 20 : DO istate = 1, nexc_search
1614 18 : IF (atom_of_state(istate) == iatom) THEN
1615 : IF (sto_state_overlap(istate) > max_overlap(iatom)*xas_control%overlap_threshold &
1616 16 : .AND. istate /= state_of_mytype(iat)) THEN
1617 6 : nexc_states(iat) = nexc_states(iat) + 1
1618 6 : state_of_atom(iat, nexc_states(iat)) = istate
1619 : END IF
1620 : END IF
1621 : END DO
1622 : END DO
1623 : END IF
1624 :
1625 : ! In the set of states, assign the index of the state to be excited for iatom
1626 40 : IF (output_unit > 0) THEN
1627 : WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") &
1628 20 : "List the atoms to be excited and the relative of MOs index "
1629 : END IF
1630 :
1631 120 : DO iat = 1, nexc_atoms
1632 80 : iatom = xas_env%exc_atoms(iat)
1633 80 : state_of_atom(iat, 1) = state_of_mytype(iatom) ! Place the state with maximum overlap first in the list
1634 80 : IF (output_unit > 0) THEN
1635 : WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A)", advance='NO') &
1636 40 : 'Atom: ', iatom, "MO index:"
1637 : END IF
1638 166 : DO istate = 1, nexc_states(iat)
1639 166 : IF (istate < nexc_states(iat)) THEN
1640 6 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(I4)", advance='NO') state_of_atom(iat, istate)
1641 : ELSE
1642 80 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(I4)") state_of_atom(iat, istate)
1643 : END IF
1644 : END DO
1645 120 : IF (state_of_atom(iat, 1) == 0 .OR. state_of_atom(iat, 1) > homo) THEN
1646 0 : CPABORT("A wrong state has been selected for excitation, check the Wannier centers")
1647 : END IF
1648 : END DO
1649 :
1650 40 : IF (xas_control%overlap_threshold < 1) THEN
1651 4 : DO iat = 1, nexc_atoms
1652 4 : IF (output_unit > 0) THEN
1653 : WRITE (UNIT=output_unit, FMT="(/,T10,A,I6)") &
1654 1 : 'Overlap integrals for Atom: ', iat
1655 5 : DO istate = 1, nexc_states(iat)
1656 : WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A,T38,f10.8)") &
1657 5 : 'State: ', state_of_atom(iat, istate), "Overlap:", sto_state_overlap(state_of_atom(iat, istate))
1658 : END DO
1659 : END IF
1660 : END DO
1661 : END IF
1662 :
1663 : ! Scales down the 2d-array to the minimal size
1664 120 : CALL reallocate(xas_env%state_of_atom, 1, nexc_atoms, 1, MAXVAL(nexc_states))
1665 :
1666 : ELSE ! Manually selected orbital indices
1667 :
1668 : ! Reallocate nexc_states and state_of_atom to include any atom
1669 2 : CALL reallocate(xas_env%nexc_states, 1, natom)
1670 2 : CALL reallocate(xas_env%state_of_atom, 1, natom, 1, SIZE(xas_control%orbital_list))
1671 2 : CALL get_xas_env(xas_env, nexc_states=nexc_states, state_of_atom=state_of_atom)
1672 :
1673 14 : nexc_states = 0
1674 30 : state_of_atom = 0
1675 : nexc_atoms = natom !To include all possible atoms in the spectrum calculation
1676 :
1677 6 : DO istate = 1, SIZE(xas_control%orbital_list)
1678 :
1679 4 : chosen_state = xas_control%orbital_list(istate)
1680 4 : nexc_atoms = 1
1681 4 : centers_wfn(1, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(1, chosen_state)
1682 4 : centers_wfn(2, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(2, chosen_state)
1683 4 : centers_wfn(3, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(3, chosen_state)
1684 :
1685 4 : distmin = 100.0_dp
1686 28 : DO iat = 1, natom
1687 96 : ra(1:3) = particle_set(iat)%r(1:3)
1688 96 : rc(1:3) = centers_wfn(1:3, chosen_state)
1689 24 : rac = pbc(ra, rc, cell)
1690 24 : dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
1691 28 : IF (dist < distmin) THEN
1692 6 : atom_of_state(chosen_state) = iat !?
1693 6 : distmin = dist
1694 : END IF
1695 : END DO ! iat
1696 :
1697 4 : nexc_states(atom_of_state(chosen_state)) = nexc_states(atom_of_state(chosen_state)) + 1
1698 6 : state_of_atom(atom_of_state(chosen_state), nexc_states(atom_of_state(chosen_state))) = chosen_state
1699 :
1700 : END DO !istate
1701 :
1702 : ! In the set of states, assign the index of the state to be excited for iatom
1703 2 : IF (output_unit > 0) THEN
1704 : WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") &
1705 1 : "List the atoms to be excited and the relative of MOs index "
1706 : END IF
1707 :
1708 14 : DO iat = 1, natom
1709 12 : IF (output_unit > 0 .AND. state_of_atom(iat, 1) /= 0) THEN
1710 : WRITE (UNIT=output_unit, FMT="(T10,A,I3,T26,A)", advance='NO') &
1711 2 : 'Atom: ', iat, "MO index:"
1712 4 : DO i = 1, nexc_states(iat)
1713 4 : IF (i < nexc_states(iat)) THEN
1714 0 : WRITE (UNIT=output_unit, FMT="(I4)", advance='NO') state_of_atom(iat, i)
1715 : ELSE
1716 2 : WRITE (UNIT=output_unit, FMT="(I4)") state_of_atom(iat, i)
1717 : END IF
1718 : END DO
1719 : END IF
1720 14 : IF (state_of_atom(iat, 1) > homo) THEN
1721 0 : CPABORT("A wrong state has been selected for excitation, check the Wannier centers")
1722 : END IF
1723 : END DO
1724 :
1725 14 : CALL reallocate(xas_env%state_of_atom, 1, natom, 1, MAXVAL(nexc_states)) ! Scales down the 2d-array to the minimal size
1726 :
1727 : END IF !Checks for manually selected orbitals from the localized set
1728 :
1729 : ! Set back the correct periodicity
1730 168 : cell%perd(1:3) = perd0(1:3)
1731 :
1732 42 : DEALLOCATE (vecbuffer)
1733 42 : DEALLOCATE (first_sgf)
1734 42 : DEALLOCATE (sto_state_overlap)
1735 42 : DEALLOCATE (max_overlap)
1736 42 : DEALLOCATE (state_of_mytype)
1737 :
1738 84 : END SUBROUTINE cls_assign_core_states
1739 :
1740 : END MODULE xas_methods
|
