(git:374b731)
Loading...
Searching...
No Matches
xas_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 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! **************************************************************************************************
15MODULE 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
31 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
32 cp_dbcsr_sm_fm_multiply,&
33 dbcsr_allocate_matrix_set
34 USE cp_external_control, ONLY: external_control
35 USE cp_fm_pool_types, ONLY: fm_pool_create_fm
36 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 cp_fm_struct_release,&
38 cp_fm_struct_type
39 USE cp_fm_types, ONLY: cp_fm_create,&
40 cp_fm_get_element,&
41 cp_fm_get_submatrix,&
42 cp_fm_release,&
43 cp_fm_set_all,&
44 cp_fm_set_submatrix,&
45 cp_fm_to_fm,&
46 cp_fm_type
47 USE cp_log_handling, ONLY: cp_get_default_logger,&
48 cp_logger_get_default_io_unit,&
49 cp_logger_type,&
50 cp_to_string
51 USE cp_output_handling, ONLY: cp_p_file,&
52 cp_print_key_finished_output,&
53 cp_print_key_should_output,&
54 cp_print_key_unit_nr
55 USE dbcsr_api, ONLY: dbcsr_convert_offsets_to_sizes,&
56 dbcsr_copy,&
57 dbcsr_create,&
58 dbcsr_distribution_type,&
59 dbcsr_p_type,&
60 dbcsr_set,&
61 dbcsr_type,&
62 dbcsr_type_antisymmetric
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
139CONTAINS
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 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 INTEGER, DIMENSION(:), POINTER :: nexc_states
172 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 REAL(dp), DIMENSION(:), POINTER :: eigenvalues
178 REAL(dp), DIMENSION(:, :), POINTER :: vecbuffer
179 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
180 TYPE(cell_type), POINTER :: cell
181 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, op_sm, ostrength_sm
185 TYPE(dbcsr_type), POINTER :: mo_coeff_b
186 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
187 TYPE(mp_para_env_type), POINTER :: para_env
188 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
189 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 CALL timeset(routinen, handle)
200
201 transition_potential = .false.
202 skip_scf = .false.
203 converged = .true.
204 should_stop = .false.
205 ch_method_flags = .false.
206
207 NULLIFY (logger)
208 logger => cp_get_default_logger()
209 output_unit = cp_logger_get_default_io_unit(logger)
210
211 NULLIFY (xas_env, groundstate_coeff, ostrength_sm, op_sm)
212 NULLIFY (excvec_coeff, qs_loc_env, cell, scf_env)
213 NULLIFY (matrix_ks)
214 NULLIFY (all_vectors, state_of_atom, nexc_states, xas_control)
215 NULLIFY (vecbuffer, op_sm, mo_coeff_b)
216 NULLIFY (dft_section, xas_section, scf_section, loc_section, print_loc_section)
217 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
218 xas_section => section_vals_get_subs_vals(dft_section, "XAS")
219 scf_section => section_vals_get_subs_vals(xas_section, "SCF")
220 loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE")
221 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 extension=".Log")
225 IF (output_unit > 0) THEN
226 WRITE (unit=output_unit, fmt="(/,T3,A,/,T25,A,/,T3,A,/)") &
227 repeat("=", 77), &
228 "START CORE LEVEL SPECTROSCOPY CALCULATION", &
229 repeat("=", 77)
230 END IF
231
232! Create the xas environment
233 CALL get_qs_env(qs_env, xas_env=xas_env)
234 IF (.NOT. ASSOCIATED(xas_env)) THEN
235 IF (output_unit > 0) THEN
236 WRITE (unit=output_unit, fmt="(/,T5,A)") &
237 "Create and initialize the xas environment"
238 END IF
239 ALLOCATE (xas_env)
240 CALL xas_env_create(xas_env)
241 CALL xas_env_init(xas_env, qs_env, dft_section, logger)
242 xas_control => dft_control%xas_control
243 CALL set_qs_env(qs_env, xas_env=xas_env)
244 END IF
245
246! Initialize the type of calculation
247 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 particle_set=particle_set)
252
253! The eigenstate of the KS Hamiltonian are nedeed
254 NULLIFY (mo_coeff, eigenvalues)
255 IF (scf_control%use_ot) THEN
256 IF (output_unit > 0) THEN
257 WRITE (unit=output_unit, fmt="(/,T10,A,/)") &
258 "Get eigenstates and eigenvalues from ground state MOs"
259 END IF
260 DO ispin = 1, dft_control%nspins
261 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
262 eigenvalues=eigenvalues, homo=homo)
263 CALL calculate_subspace_eigenvalues(mo_coeff, &
264 matrix_ks(ispin)%matrix, eigenvalues, &
265 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 added_mos = scf_control%added_mos
270 NULLIFY (scf_control)
271! Consider to use get function for this
272 CALL get_xas_env(xas_env, scf_control=scf_control)
273 scf_control%added_mos = added_mos
274
275! Set initial occupation numbers, and store the original ones
276 my_homo = 0
277 my_nelectron = 0
278 DO ispin = 1, dft_control%nspins
279 CALL get_mo_set(mos(ispin), nelectron=my_nelectron(ispin), maxocc=maxocc, &
280 homo=my_homo(ispin), uniform_occupation=my_uocc(ispin))
281 END DO
282
283 nspins = dft_control%nspins
284! at the moment the only implemented method for XAS and XES calculations
285 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 IF (nspins == 1 .AND. transition_potential) THEN
291 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 spin_channel=my_spin)
299
300! Set of states among which there is the state to be excited
301 CALL get_mo_set(mos(my_spin), nao=nao, homo=homo)
302 IF (xas_control%nexc_search < 0) xas_control%nexc_search = homo
303 nexc_search = xas_control%nexc_search
304
305 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 CALL get_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env)
309 IF (qs_loc_env%do_localize) THEN
310 IF (output_unit > 0) THEN
311 WRITE (unit=output_unit, fmt="(/,T2,A34,I3,A36/)") &
312 "Localize a sub-set of MOs of spin ", my_spin, ","// &
313 " to better identify the core states"
314 IF ( &
315 qs_loc_env%localized_wfn_control%set_of_states == state_loc_range) THEN
316 WRITE (unit=output_unit, fmt="( A , I7, A, I7)") " The sub-set contains states from ", &
317 qs_loc_env%localized_wfn_control%lu_bound_states(1, my_spin), " to ", &
318 qs_loc_env%localized_wfn_control%lu_bound_states(2, my_spin)
319 ELSEIF (qs_loc_env%localized_wfn_control%set_of_states == state_loc_list) THEN
320 WRITE (unit=output_unit, fmt="( A )") " The sub-set contains states given in the input list"
321 END IF
322
323 END IF
324 CALL qs_loc_driver(qs_env, qs_loc_env, print_loc_section, myspin=my_spin)
325 END IF
326
327 cpassert(ASSOCIATED(groundstate_coeff))
328 DO ispin = 1, nspins
329 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, nmo=nmo)
330 CALL cp_fm_to_fm(mo_coeff, groundstate_coeff(ispin), nmo, 1, 1)
331 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 IF (xas_control%xas_method == xes_tp_val .AND. &
340 xas_control%xes_core_occupation == 1.0_dp) THEN
341 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,/,T10,A)') &
342 "START Core Level Spectroscopy Calculation for the Emission Spectrum"
343 IF (xas_control%xes_homo_occupation == 1) THEN
344 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(T10,A,/,A)') &
345 "The core state is fully occupied and XES from ground state calculation.", &
346 " No SCF is needed, MOS already available"
347 ELSE IF (xas_control%xes_homo_occupation == 0) THEN
348 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(T10,A,/,A)') &
349 "The core state is fully occupied and the homo is empty", &
350 " (final state of the core hole decay). Only one SCF is needed (not one per atom)"
351 END IF
352 skip_scf = .true.
353
354 CALL set_xas_env(xas_env=xas_env, xas_estate=-1, homo_occ=xas_control%xes_homo_occupation)
355 CALL xes_scf_once(qs_env, xas_env, converged, should_stop)
356
357 IF (converged .AND. .NOT. should_stop .AND. xas_control%xes_homo_occupation == 0) THEN
358 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,T10,A,I6)') &
359 "SCF with empty homo converged "
360 ELSE IF (.NOT. converged .OR. should_stop) THEN
361 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,T10,A,I6)') &
362 "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 DO ispin = 1, dft_control%nspins
370 CALL set_mo_set(mos(ispin), homo=my_homo(ispin), &
371 uniform_occupation=my_uocc(ispin), nelectron=my_nelectron(ispin))
372 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
373 CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
374 END DO
375
376 IF (output_unit > 0) THEN
377 WRITE (unit=output_unit, fmt="(/,T3,A,/,T25,A,/,T3,A,/)") &
378 repeat("=", 77), &
379 "END CORE LEVEL SPECTROSCOPY CALCULATION", &
380 repeat("=", 77)
381 END IF
382
383 CALL xas_env_release(qs_env%xas_env)
384 DEALLOCATE (qs_env%xas_env)
385 NULLIFY (qs_env%xas_env)
386
387 CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
388 "PRINT%PROGRAM_RUN_INFO")
389 CALL timestop(handle)
390 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 qs_env)
398 CALL get_xas_env(xas_env=xas_env, &
399 state_of_atom=state_of_atom, nexc_states=nexc_states)
400
401 IF (skip_scf) THEN
402 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 source_start=1, target_start=1)
405 END IF
406
407 ALLOCATE (vecbuffer(1, nao))
408 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 CALL get_xas_env(xas_env=xas_env, ostrength_sm=ostrength_sm)
414 DO i = 1, 3
415 NULLIFY (op_sm(i)%matrix)
416 op_sm(i)%matrix => ostrength_sm(i)%matrix
417 END DO
418 IF (xas_control%dipole_form == xas_dip_vel) THEN
419 CALL p_xyz_ao(op_sm, qs_env)
420 END IF
421 END IF
422
423 ! DO SCF if required
424 DO iat = 1, nexc_atoms
425 iatom = xas_env%exc_atoms(iat)
426 DO istate = 1, nexc_states(iat)
427 ! determine which state has to be excited in the global list
428 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 CALL get_mo_set(mos(my_spin), nmo=nmo)
432 CALL get_xas_env(xas_env, occ_estate=occ_estate, xas_nelectron=xas_nelectron)
433 tmp = xas_nelectron + 1.0_dp - occ_estate
434 IF (nmo < tmp) &
435 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 IF (xas_control%xas_restart) THEN
442 CALL xas_read_restart(xas_env, xas_section, qs_env, xas_control%xas_method, iatom, &
443 state_to_be_excited, istate)
444 END IF
445 CALL set_xas_env(xas_env=xas_env, xas_estate=state_to_be_excited)
446 CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff)
447 cpassert(ASSOCIATED(excvec_coeff))
448 CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, state_to_be_excited, &
449 nao, 1, transpose=.true.)
450 CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, &
451 nao, 1, transpose=.true.)
452
453 IF (transition_potential) THEN
454
455 IF (.NOT. skip_scf) THEN
456 IF (output_unit > 0) THEN
457 WRITE (unit=output_unit, fmt='(/,T5,A)') repeat("-", 75)
458 IF (xas_control%xas_method == xas_dscf) THEN
459 WRITE (unit=output_unit, fmt='(/,/,T10,A,I6)') &
460 "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 "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 "Excited state", istate, "out of", nexc_states(iat)
466 WRITE (unit=output_unit, fmt='(T10,A,T50,f10.4)') "Occupation of the core orbital", &
467 occ_estate
468 WRITE (unit=output_unit, fmt='(T10,A28,I3, T50,F10.4)') "Number of electrons in Spin ", &
469 my_spin, xas_nelectron
470 END IF
471 END IF
472
473 CALL get_xas_env(xas_env=xas_env, scf_env=scf_env)
474 IF (.NOT. ASSOCIATED(scf_env)) THEN
475 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 CALL set_xas_env(xas_env, scf_env=scf_env)
478 ELSE
479 CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
480 END IF
481
482 DO ispin = 1, SIZE(mos)
483 IF (ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN !fm->dbcsr
484 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
485 mos(ispin)%mo_coeff_b) !fm->dbcsr
486 END IF !fm->dbcsr
487 END DO !fm->dbcsr
488
489 IF (.NOT. scf_env%skip_diis) THEN
490 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
491 ALLOCATE (scf_env%scf_diis_buffer)
492 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
493 END IF
494 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 xas_section, scf_section, converged, should_stop)
499
500 CALL external_control(should_stop, "CLS", target_time=qs_env%target_time, &
501 start_time=qs_env%start_time)
502 IF (should_stop) THEN
503 CALL scf_env_cleanup(scf_env)
504 EXIT
505 END IF
506
507 END IF
508 ! SCF DONE
509
510 ! Write last wavefunction to screen
511 IF (SIZE(mos) > 1) THEN
512 CALL write_mo_set_to_output_unit(mos(1), atomic_kind_set, qs_kind_set, particle_set, &
513 dft_section, 4, 0, final_mos=.false., spin="XAS ALPHA")
514 CALL write_mo_set_to_output_unit(mos(2), atomic_kind_set, qs_kind_set, particle_set, &
515 dft_section, 4, 0, final_mos=.false., spin="XAS BETA")
516 ELSE
517 CALL write_mo_set_to_output_unit(mos(1), atomic_kind_set, qs_kind_set, particle_set, &
518 dft_section, 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 IF (converged) THEN
529 CALL cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, &
530 iatom, istate)
531 ELSE
532 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,/,T10,A,I6)') &
533 "SCF with core hole NOT converged for ATOM ", iatom
534 END IF
535
536 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 CALL get_qs_env(qs_env, mos=mos)
541 DO ispin = 1, nspins
542 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
543 CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
544 END DO
545 IF (iat == nexc_atoms) THEN
546 CALL scf_env_cleanup(scf_env)
547 CALL scf_env_release(xas_env%scf_env)
548 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 IF (ASSOCIATED(vecbuffer)) THEN
559 DEALLOCATE (vecbuffer)
560 DEALLOCATE (op_sm)
561 END IF
562
563 DO ispin = 1, dft_control%nspins
564 CALL set_mo_set(mos(ispin), homo=my_homo(ispin), &
565 uniform_occupation=my_uocc(ispin), nelectron=my_nelectron(ispin))
566 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
567 CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
568 END DO
569
570 IF (output_unit > 0) THEN
571 WRITE (unit=output_unit, fmt="(/,T3,A,/,T25,A,/,T3,A,/)") &
572 repeat("=", 77), &
573 "END CORE LEVEL SPECTROSCOPY CALCULATION", &
574 repeat("=", 77)
575 END IF
576
577 CALL xas_env_release(qs_env%xas_env)
578 DEALLOCATE (qs_env%xas_env)
579 NULLIFY (qs_env%xas_env)
580
581 CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
582 "PRINT%PROGRAM_RUN_INFO")
583 CALL timestop(handle)
584
585 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_type_tmp, kind_z_tmp, &
609 last_sgf
610 INTEGER, DIMENSION(4, 7) :: ne
611 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 REAL(dp), DIMENSION(:), POINTER :: sto_zet
616 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 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 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 POINTER :: sab_orb
628 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
629 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 n_mo(1:2) = 0
638 cpassert(ASSOCIATED(xas_env))
639
640 NULLIFY (atomic_kind_set, qs_kind_set, dft_control, scf_control, matrix_s, mos, mpools)
641 NULLIFY (para_env, particle_set, xas_control)
642 NULLIFY (qs_loc_env)
643 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 dbcsr_dist=dbcsr_dist)
653
654 xas_section => section_vals_get_subs_vals(dft_section, "XAS")
655 ALLOCATE (dft_control%xas_control)
656 CALL xas_control_create(dft_control%xas_control)
657 CALL read_xas_control(dft_control%xas_control, xas_section)
658 CALL write_xas_control(dft_control%xas_control, dft_section)
659 xas_control => dft_control%xas_control
660 ALLOCATE (scf_control)
661 CALL scf_c_create(scf_control)
662 CALL scf_c_read_parameters(scf_control, xas_section)
663 CALL set_xas_env(xas_env, scf_control=scf_control)
664
665 my_spin = xas_control%spin_channel
666 nexc_search = xas_control%nexc_search
667 IF (nexc_search < 0) THEN
668 ! ground state occupation
669 CALL get_mo_set(mos(my_spin), nmo=nmo, lfomo=lfomo)
670 nexc_search = lfomo - 1
671 END IF
672 nexc_atoms = xas_control%nexc_atoms
673 ALLOCATE (xas_env%exc_atoms(nexc_atoms))
674 xas_env%exc_atoms = xas_control%exc_atoms
675 CALL set_xas_env(xas_env=xas_env, nexc_search=nexc_search, &
676 nexc_atoms=nexc_atoms, spin_channel=my_spin)
677
678 CALL mpools_get(mpools, ao_mo_fm_pools=xas_env%ao_mo_fm_pools)
679
680 NULLIFY (mo_coeff)
681 CALL get_mo_set(mos(my_spin), nao=nao, homo=homo, nmo=nmo, mo_coeff=mo_coeff, nelectron=nelectron)
682
683 nvirtual2 = 0
684 IF (xas_control%added_mos .GT. 0) THEN
685 nvirtual2 = min(xas_control%added_mos, nao - nmo)
686 xas_env%unoccupied_eps = xas_control%eps_added
687 xas_env%unoccupied_max_iter = xas_control%max_iter_added
688 END IF
689 nvirtual = nmo + nvirtual2
690
691 n_mo(1:2) = nmo
692
693 ALLOCATE (xas_env%centers_wfn(3, nexc_search))
694 ALLOCATE (xas_env%atom_of_state(nexc_search))
695 ALLOCATE (xas_env%type_of_state(nexc_search))
696 ALLOCATE (xas_env%state_of_atom(nexc_atoms, nexc_search))
697 ALLOCATE (xas_env%nexc_states(nexc_atoms))
698 ALLOCATE (xas_env%mykind_of_atom(nexc_atoms))
699 nkind = SIZE(atomic_kind_set, 1)
700 ALLOCATE (xas_env%mykind_of_kind(nkind))
701 xas_env%mykind_of_kind = 0
702
703 ! create a new matrix structure nao x 1
704 NULLIFY (tmp_fm_struct)
705 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
706 ncol_global=1, para_env=para_env, context=mo_coeff%matrix_struct%context)
707 ALLOCATE (xas_env%excvec_coeff)
708 CALL cp_fm_create(xas_env%excvec_coeff, tmp_fm_struct)
709 CALL cp_fm_struct_release(tmp_fm_struct)
710
711 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 context=mo_coeff%matrix_struct%context)
715 ALLOCATE (xas_env%excvec_overlap)
716 CALL cp_fm_create(xas_env%excvec_overlap, tmp_fm_struct)
717 CALL cp_fm_struct_release(tmp_fm_struct)
718
719 nspins = SIZE(mos, 1)
720
721 ! initialize operators for the calculation of the oscillator strengths
722 IF (xas_control%xas_method == xas_tp_hh) THEN
723 occ_estate = 0.5_dp
724 nele = real(nelectron, dp) - 0.5_dp
725 occ_homo = 1.0_dp
726 occ_homo_plus = 0._dp
727 ELSEIF (xas_control%xas_method == xas_tp_xhh) THEN
728 occ_estate = 0.5_dp
729 nele = real(nelectron, dp)
730 occ_homo = 1.0_dp
731 occ_homo_plus = 0.5_dp
732 ELSEIF (xas_control%xas_method == xas_tp_fh) THEN
733 occ_estate = 0.0_dp
734 nele = real(nelectron, dp) - 1.0_dp
735 occ_homo = 1.0_dp
736 occ_homo_plus = 0._dp
737 ELSEIF (xas_control%xas_method == xas_tp_xfh) THEN
738 occ_estate = 0.0_dp
739 nele = real(nelectron, dp)
740 occ_homo = 1.0_dp
741 occ_homo_plus = 1._dp
742 ELSEIF (xas_control%xas_method == xes_tp_val) THEN
743 occ_estate = xas_control%xes_core_occupation
744 nele = real(nelectron, dp) - xas_control%xes_core_occupation
745 occ_homo = xas_control%xes_homo_occupation
746 ELSEIF (xas_control%xas_method == xas_dscf) THEN
747 occ_estate = 0.0_dp
748 nele = real(nelectron, dp)
749 occ_homo = 1.0_dp
750 occ_homo_plus = 1._dp
751 ELSEIF (xas_control%xas_method == xas_tp_flex) THEN
752 nele = real(xas_control%nel_tot, dp)
753 occ_estate = real(xas_control%xas_core_occupation, dp)
754 IF (nele < 0.0_dp) nele = real(nelectron, dp) - (1.0_dp - occ_estate)
755 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 nvirtual2=nvirtual2, nvirtual=nvirtual, homo_occ=occ_homo)
759
760 ! Initialize the list of orbitals for cube files printing
761 IF (btest(cp_print_key_should_output(logger%iter_info, xas_section, &
762 "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
763 NULLIFY (bounds, list)
764 CALL section_vals_val_get(xas_section, &
765 "PRINT%CLS_FUNCTION_CUBES%CUBES_LU_BOUNDS", &
766 i_vals=bounds)
767 ncubes = bounds(2) - bounds(1) + 1
768 IF (ncubes > 0) THEN
769 ALLOCATE (xas_control%list_cubes(ncubes))
770
771 DO ik = 1, ncubes
772 xas_control%list_cubes(ik) = bounds(1) + (ik - 1)
773 END DO
774 END IF
775
776 IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
777 CALL section_vals_val_get(xas_section, &
778 "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
779 n_rep_val=n_rep)
780 ncubes = 0
781 DO ik = 1, n_rep
782 NULLIFY (list)
783 CALL section_vals_val_get(xas_section, &
784 "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
785 i_rep_val=ik, i_vals=list)
786 IF (ASSOCIATED(list)) THEN
787 CALL reallocate(xas_control%list_cubes, 1, ncubes + SIZE(list))
788 DO i = 1, SIZE(list)
789 xas_control%list_cubes(i + ncubes) = list(i)
790 END DO
791 ncubes = ncubes + SIZE(list)
792 END IF
793 END DO ! ik
794 END IF
795
796 IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
797 ncubes = max(10, xas_control%added_mos/10)
798 ncubes = min(ncubes, xas_control%added_mos)
799 ALLOCATE (xas_control%list_cubes(ncubes))
800 DO ik = 1, ncubes
801 xas_control%list_cubes(ik) = homo + ik
802 END DO
803 END IF
804 ELSE
805 NULLIFY (xas_control%list_cubes)
806 END IF
807
808 NULLIFY (tmp_fm_struct)
809 ALLOCATE (xas_env%groundstate_coeff(nspins))
810 DO ispin = 1, nspins
811 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 name="xas_env%mo0"//trim(adjustl(cp_to_string(ispin))))
815 END DO ! ispin
816
817 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 context=mo_coeff%matrix_struct%context)
821 ALLOCATE (xas_env%dip_fm_set(2, 3))
822 DO i = 1, 3
823 DO j = 1, 2
824 CALL cp_fm_create(xas_env%dip_fm_set(j, i), tmp_fm_struct)
825 END DO
826 END DO
827 CALL cp_fm_struct_release(tmp_fm_struct)
828
829 !Array to store all the eigenstates: occupied and the required not occupied
830 IF (nvirtual2 .GT. 0) THEN
831 ALLOCATE (xas_env%unoccupied_evals(nvirtual2))
832 NULLIFY (tmp_fm_struct)
833 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
834 ncol_global=nvirtual2, &
835 para_env=para_env, context=mo_coeff%matrix_struct%context)
836 ALLOCATE (xas_env%unoccupied_orbs)
837 CALL cp_fm_create(xas_env%unoccupied_orbs, tmp_fm_struct)
838 CALL cp_fm_struct_release(tmp_fm_struct)
839 END IF
840
841 NULLIFY (tmp_fm_struct)
842 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
843 ncol_global=nvirtual, &
844 para_env=para_env, context=mo_coeff%matrix_struct%context)
845 ALLOCATE (xas_env%all_vectors)
846 CALL cp_fm_create(xas_env%all_vectors, tmp_fm_struct)
847 CALL cp_fm_struct_release(tmp_fm_struct)
848
849 ! Array to store all the energies needed for the spectrum
850 ALLOCATE (xas_env%all_evals(nvirtual))
851
852 IF (xas_control%dipole_form == xas_dip_len) THEN
853 CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
854 DO i = 1, 3
855 ALLOCATE (xas_env%ostrength_sm(i)%matrix)
856 CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, matrix_s(1)%matrix, &
857 "xas_env%ostrength_sm-"//trim(adjustl(cp_to_string(i))))
858 CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
859 END DO
860 ELSEIF (xas_control%dipole_form == xas_dip_vel) THEN
861 !
862 ! prepare for allocation
863 natom = SIZE(particle_set, 1)
864 ALLOCATE (first_sgf(natom))
865 ALLOCATE (last_sgf(natom))
866 CALL get_particle_set(particle_set, qs_kind_set, &
867 first_sgf=first_sgf, &
868 last_sgf=last_sgf)
869 ALLOCATE (row_blk_sizes(natom))
870 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
871 DEALLOCATE (first_sgf)
872 DEALLOCATE (last_sgf)
873 !
874 !
875 CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
876 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 nze=0, mutable_work=.true.)
882 CALL cp_dbcsr_alloc_block_from_nbl(xas_env%ostrength_sm(1)%matrix, sab_orb)
883 CALL dbcsr_set(xas_env%ostrength_sm(1)%matrix, 0.0_dp)
884 DO i = 2, 3
885 ALLOCATE (xas_env%ostrength_sm(i)%matrix)
886 CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, xas_env%ostrength_sm(1)%matrix, &
887 "xas_env%ostrength_sm-"//trim(adjustl(cp_to_string(i))))
888 CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
889 END DO
890
891 DEALLOCATE (row_blk_sizes)
892 END IF
893
894 ! Define the qs_loc_env : to find centers, spread and possibly localize them
895 IF (.NOT. (ASSOCIATED(xas_env%qs_loc_env))) THEN
896 ALLOCATE (qs_loc_env)
897 CALL qs_loc_env_create(qs_loc_env)
898 CALL set_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env)
899 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 do_xas=.true., nloc_xas=nexc_search, spin_xas=my_spin)
903
904 IF (.NOT. qs_loc_env%do_localize) THEN
905 qs_loc_env%localized_wfn_control%localization_method = do_loc_none
906
907 ELSE
908 nmoloc = qs_loc_env%localized_wfn_control%nloc_states
909 CALL set_loc_wfn_lists(qs_loc_env%localized_wfn_control, nmoloc, n_mo, nspins, my_spin)
910 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 qs_env, myspin=my_spin, do_localize=qs_loc_env%do_localize)
913 END IF
914 END IF
915
916 !Type of state
917 ALLOCATE (nq(1), lq(1), sto_zet(1))
918 IF (xas_control%state_type == xas_1s_type) THEN
919 nq(1) = 1
920 lq(1) = 0
921 ELSEIF (xas_control%state_type == xas_2s_type) THEN
922 nq(1) = 2
923 lq(1) = 0
924 ELSEIF (xas_control%state_type == xas_2p_type) THEN
925 nq(1) = 2
926 lq(1) = 1
927 ELSEIF (xas_control%state_type == xas_3s_type) THEN
928 nq(1) = 3
929 lq(1) = 0
930 ELSEIF (xas_control%state_type == xas_3p_type) THEN
931 nq(1) = 3
932 lq(1) = 1
933 ELSEIF (xas_control%state_type == xas_3d_type) THEN
934 nq(1) = 3
935 lq(1) = 2
936 ELSEIF (xas_control%state_type == xas_4s_type) THEN
937 nq(1) = 4
938 lq(1) = 0
939 ELSEIF (xas_control%state_type == xas_4p_type) THEN
940 nq(1) = 4
941 lq(1) = 1
942 ELSEIF (xas_control%state_type == xas_4d_type) THEN
943 nq(1) = 4
944 lq(1) = 2
945 ELSEIF (xas_control%state_type == xas_4f_type) THEN
946 nq(1) = 4
947 lq(1) = 3
948 ELSE
949 cpabort("XAS type of state not implemented")
950 END IF
951
952! Find core orbitals of right angular momentum
953 ALLOCATE (kind_type_tmp(nkind))
954 ALLOCATE (kind_z_tmp(nkind))
955 kind_type_tmp = 0
956 kind_z_tmp = 0
957 nk = 0
958 DO iat = 1, nexc_atoms
959 iatom = xas_env%exc_atoms(iat)
960 NULLIFY (atomic_kind)
961 atomic_kind => particle_set(iatom)%atomic_kind
962 CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=ikind)
963 CALL get_qs_kind(qs_kind_set(ikind), zeff=zatom)
964 ihavethis = .false.
965 DO ik = 1, nk
966 IF (ikind == kind_type_tmp(ik)) THEN
967 ihavethis = .true.
968 xas_env%mykind_of_atom(iat) = ik
969 EXIT
970 END IF
971 END DO
972 IF (.NOT. ihavethis) THEN
973 nk = nk + 1
974 kind_type_tmp(nk) = ikind
975 kind_z_tmp(nk) = int(zatom)
976 xas_env%mykind_of_atom(iat) = nk
977 xas_env%mykind_of_kind(ikind) = nk
978 END IF
979 END DO ! iat
980
981 ALLOCATE (xas_env%my_gto_basis(nk))
982 ALLOCATE (xas_env%stogto_overlap(nk))
983 DO ik = 1, nk
984 NULLIFY (xas_env%my_gto_basis(ik)%gto_basis_set, sto_basis_set)
985 ne = 0
986 DO l = 1, lq(1) + 1
987 nj = 2*(l - 1) + 1
988 DO i = l, nq(1)
989 ne(l, i) = ptable(kind_z_tmp(ik))%e_conv(l - 1) - 2*nj*(i - l)
990 ne(l, i) = max(ne(l, i), 0)
991 ne(l, i) = min(ne(l, i), 2*nj)
992 END DO
993 END DO
994
995 sto_zet(1) = srules(kind_z_tmp(ik), ne, nq(1), lq(1))
996 CALL allocate_sto_basis_set(sto_basis_set)
997 name_sto = 'xas_tmp_sto'
998 CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, &
999 lq=lq, zet=sto_zet, name=name_sto)
1000 CALL create_gto_from_sto_basis(sto_basis_set, &
1001 xas_env%my_gto_basis(ik)%gto_basis_set, xas_control%ngauss)
1002 CALL deallocate_sto_basis_set(sto_basis_set)
1003 xas_env%my_gto_basis(ik)%gto_basis_set%norm_type = 2
1004 CALL init_orb_basis_set(xas_env%my_gto_basis(ik)%gto_basis_set)
1005
1006 ikind = kind_type_tmp(ik)
1007 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1008
1009 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf_gto)
1010 CALL get_gto_basis_set(gto_basis_set=xas_env%my_gto_basis(ik)%gto_basis_set, nsgf=nsgf_sto)
1011 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 xas_env%stogto_overlap(ik)%array)
1015 END DO
1016
1017 DEALLOCATE (nq, lq, sto_zet)
1018 DEALLOCATE (kind_type_tmp, kind_z_tmp)
1019
1020 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 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 REAL(dp), DIMENSION(:), POINTER :: all_evals
1052 REAL(dp), DIMENSION(:, :), POINTER :: sp_ab, sp_em
1053 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_sm, ostrength_sm
1057 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1058 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1059
1060 NULLIFY (logger)
1061 logger => cp_get_default_logger()
1062 output_unit = cp_logger_get_default_io_unit(logger)
1063
1064 NULLIFY (ostrength_sm, op_sm, dip_fm_set)
1065 NULLIFY (all_evals, all_vectors, excvec_coeff)
1066 NULLIFY (mos, particle_set, sp_em, sp_ab)
1067 ALLOCATE (op_sm(3))
1068
1069 CALL get_qs_env(qs_env=qs_env, &
1070 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 ostrength_sm=ostrength_sm, nvirtual=nvirtual, spin_channel=my_spin)
1075 CALL get_mo_set(mos(my_spin), homo=homo, lfomo=lfomo, nmo=nmo)
1076
1077 nabs = nvirtual - lfomo + 1
1078 ALLOCATE (sp_em(6, homo))
1079 ALLOCATE (sp_ab(6, nabs))
1080 cpassert(ASSOCIATED(excvec_coeff))
1081
1082 IF (.NOT. xas_control%xas_method == xas_dscf) THEN
1083 ! Calculate the spectrum
1084 IF (xas_control%dipole_form == xas_dip_len) THEN
1085 rc(1:3) = particle_set(iatom)%r(1:3)
1086 DO i = 1, 3
1087 NULLIFY (op_sm(i)%matrix)
1088 op_sm(i)%matrix => ostrength_sm(i)%matrix
1089 END DO
1090 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 sp_em, sp_ab, xas_estate, nvirtual, my_spin)
1094 DO i = 1, SIZE(ostrength_sm, 1)
1095 CALL dbcsr_set(ostrength_sm(i)%matrix, 0.0_dp)
1096 END DO
1097 ELSE
1098 DO i = 1, 3
1099 NULLIFY (op_sm(i)%matrix)
1100 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 sp_em, sp_ab, xas_estate, nvirtual, my_spin)
1105 END IF
1106 END IF
1107
1108 CALL get_mo_set(mos(my_spin), lfomo=lfomo)
1109 ! write the spectrum, if the file exists it is appended
1110 IF (.NOT. xas_control%xas_method == xas_dscf) THEN
1111 length = (.NOT. xas_control%dipole_form == xas_dip_vel)
1112 CALL xas_write(sp_em, sp_ab, xas_estate, &
1113 xas_section, iatom, istate, lfomo, length=length)
1114 END IF
1115
1116 DEALLOCATE (sp_em)
1117 DEALLOCATE (sp_ab)
1118
1119 IF (btest(cp_print_key_should_output(logger%iter_info, xas_section, &
1120 "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
1121 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 iatom, append_cube)
1124 END IF
1125
1126 IF (btest(cp_print_key_should_output(logger%iter_info, xas_section, &
1127 "PRINT%PDOS"), cp_p_file)) THEN
1128 CALL xas_pdos(qs_env, xas_section, mos, iatom)
1129 END IF
1130
1131 DEALLOCATE (op_sm)
1132
1133 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 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 NULLIFY (logger)
1166 logger => cp_get_default_logger()
1167
1168 my_pos = "APPEND"
1169 my_act = "WRITE"
1170
1171 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 file_form="FORMATTED", middle_name=trim(mittle_em))
1176
1177 IF (out_sp_em > 0) THEN
1178 WRITE (out_sp_em, '(A,I6,A,I6,A,I6)') " Emission spectrum for atom ", iatom, &
1179 ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_em, 2)
1180 ene2 = 1.0_dp
1181 DO istate = estate, SIZE(sp_em, 2)
1182 IF (length) ene2 = sp_em(1, istate)*sp_em(1, istate)
1183 WRITE (out_sp_em, '(I6,5F16.8,F10.5)') istate, sp_em(1, istate)*evolt, &
1184 sp_em(2, istate)*ene2, sp_em(3, istate)*ene2, &
1185 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 "PRINT%XES_SPECTRUM")
1190
1191 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 file_form="FORMATTED", middle_name=trim(mittle_ab))
1195
1196 IF (out_sp_ab > 0) THEN
1197 WRITE (out_sp_ab, '(A,I6,A,I6,A,I6)') " Absorption spectrum for atom ", iatom, &
1198 ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_ab, 2)
1199 ene2 = 1.0_dp
1200 DO i = 1, SIZE(sp_ab, 2)
1201 istate = lfomo - 1 + i
1202 IF (length) ene2 = sp_ab(1, i)*sp_ab(1, i)
1203 WRITE (out_sp_ab, '(I6,5F16.8,F10.5)') istate, sp_ab(1, i)*evolt, &
1204 sp_ab(2, i)*ene2, sp_ab(3, i)*ene2, &
1205 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 "PRINT%XAS_SPECTRUM")
1211
1212 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 SUBROUTINE xas_print_cubes(xas_control, qs_env, xas_section, &
1228 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 REAL(dp), DIMENSION(:, :), POINTER :: centers
1241 TYPE(section_vals_type), POINTER :: print_key
1242
1243 nspins = SIZE(mos)
1244
1245 print_key => section_vals_get_subs_vals(xas_section, "PRINT%CLS_FUNCTION_CUBES")
1246 my_mittle = 'at'//trim(adjustl(cp_to_string(iatom)))
1247 nstates = SIZE(xas_control%list_cubes, 1)
1248
1249 IF (xas_control%do_centers) THEN
1250 ! one might like to calculate the centers of the xas orbital (without localizing them)
1251 ELSE
1252 ALLOCATE (centers(6, nstates))
1253 centers = 0.0_dp
1254 END IF
1255 my_spin = xas_control%spin_channel
1256
1257 CALL get_mo_set(mos(my_spin), homo=homo)
1258 istate0 = 0
1259
1260 my_pos = "REWIND"
1261 IF (append_cube) THEN
1262 my_pos = "APPEND"
1263 END IF
1264
1265 CALL qs_print_cubes(qs_env, all_vectors, nstates, xas_control%list_cubes, &
1266 centers, print_key, my_mittle, state0=istate0, file_position=my_pos)
1267
1268 DEALLOCATE (centers)
1269
1270 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 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 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1293 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1294 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1295
1296 NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
1297 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 qs_kind_set=qs_kind_set)
1303
1304 DO ispin = 1, 2
1305 CALL calculate_projected_dos(mos(ispin), atomic_kind_set, qs_kind_set, particle_set, qs_env, &
1306 xas_section, ispin, xas_mittle)
1307 END DO
1308
1309 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 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 REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
1347 TYPE(cp_fm_type) :: fm_work
1348
1349 cpassert(ASSOCIATED(fm_set))
1350 NULLIFY (eigenvalues, occupation_numbers)
1351
1352 CALL get_mo_set(mos(my_spin), eigenvalues=eigenvalues, occupation_numbers=occupation_numbers, &
1353 nao=nao, nmo=nmo, homo=homo, lfomo=lfomo)
1354
1355 CALL cp_fm_create(fm_work, all_vectors%matrix_struct)
1356 DO i = 1, SIZE(fm_set, 2)
1357 CALL cp_fm_set_all(fm_set(my_spin, i), 0.0_dp)
1358 CALL cp_fm_set_all(fm_work, 0.0_dp)
1359 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 fm_work, 0.0_dp, fm_set(my_spin, i), b_first_col=1)
1362 END DO
1363 CALL cp_fm_release(fm_work)
1364
1365 sp_em = 0.0_dp
1366 sp_ab = 0.0_dp
1367 ene_i = eigenvalues(estate)
1368 DO istate = 1, nstate
1369 ene_f = all_evals(istate)
1370 DO i = 1, 3
1371 CALL cp_fm_get_element(fm_set(my_spin, i), 1, istate, dip(i))
1372 END DO
1373 IF (istate <= homo) THEN
1374 sp_em(1, istate) = ene_f - ene_i
1375 sp_em(2, istate) = dip(1)
1376 sp_em(3, istate) = dip(2)
1377 sp_em(4, istate) = dip(3)
1378 sp_em(5, istate) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
1379 sp_em(6, istate) = occupation_numbers(istate)
1380 END IF
1381 IF (istate >= lfomo) THEN
1382 i_abs = istate - lfomo + 1
1383 sp_ab(1, i_abs) = ene_f - ene_i
1384 sp_ab(2, i_abs) = dip(1)
1385 sp_ab(3, i_abs) = dip(2)
1386 sp_ab(4, i_abs) = dip(3)
1387 sp_ab(5, i_abs) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
1388 IF (istate <= nmo) sp_ab(6, i_abs) = occupation_numbers(istate)
1389 END IF
1390
1391 END DO
1392
1393 END SUBROUTINE spectrum_dip_vel
1394
1395! **************************************************************************************************
1396!> \brief ...
1397!> \param base_a ...
1398!> \param base_b ...
1399!> \param matrix ...
1400! **************************************************************************************************
1401 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1410 npgfb, nsgfa_set, nsgfb_set
1411 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1412 REAL(dp) :: rab(3)
1413 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
1414 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, sphi_a, &
1415 sphi_b, zeta, zetb
1416
1417 NULLIFY (la_max, la_min, lb_max, lb_min)
1418 NULLIFY (npgfa, npgfb, nsgfa_set, nsgfb_set)
1419 NULLIFY (first_sgfa, first_sgfb)
1420 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 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 maxco=maxcob, nset=nsetb, maxl=maxlb)
1431 ! Initialize and allocate
1432 rab = 0.0_dp
1433 matrix = 0.0_dp
1434
1435 ldsab = max(maxcoa, maxcob, nsgfa, nsgfb)
1436 maxl = max(maxla, maxlb)
1437
1438 ALLOCATE (sab(ldsab, ldsab))
1439 ALLOCATE (work(ldsab, ldsab))
1440
1441 DO iset = 1, nseta
1442
1443 na = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1444 sgfa = first_sgfa(1, iset)
1445
1446 DO jset = 1, nsetb
1447 nb = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1448 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 rab, sab)
1453 CALL contraction(sab, work, ca=scon_a(:, sgfa:), na=na, ma=nsgfa_set(iset), &
1454 cb=scon_b(:, sgfb:), nb=nb, mb=nsgfb_set(jset))
1455 CALL block_add("IN", work, nsgfa_set(iset), nsgfb_set(jset), matrix, sgfa, sgfb)
1456
1457 END DO ! jset
1458 END DO ! iset
1459 DEALLOCATE (sab, work)
1460
1461 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 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 INTEGER, DIMENSION(:), POINTER :: atom_of_state, mykind_of_kind, &
1496 nexc_states, state_of_mytype, &
1497 type_of_state
1498 INTEGER, DIMENSION(:, :), POINTER :: state_of_atom
1499 REAL(dp) :: component, dist, distmin, maxocc, ra(3), &
1500 rac(3), rc(3)
1501 REAL(dp), DIMENSION(:), POINTER :: max_overlap, sto_state_overlap
1502 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 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 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1510 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1511 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1512
1513 NULLIFY (cell, mos, particle_set)
1514 NULLIFY (atom_of_state, centers_wfn, mykind_of_kind, state_of_atom, nexc_states)
1515 NULLIFY (stogto_overlap, type_of_state, max_overlap, qs_kind_set)
1516 NULLIFY (state_of_mytype, type_of_state, sto_state_overlap)
1517
1518 NULLIFY (logger)
1519 logger => cp_get_default_logger()
1520 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 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 perd0(1:3) = cell%perd(1:3)
1528 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 spin_channel=my_spin, nexc_search=nexc_search, nexc_states=nexc_states)
1536
1537 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 ALLOCATE (vecbuffer(1, nao))
1541 natom = SIZE(particle_set)
1542
1543 ALLOCATE (first_sgf(natom))
1544 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
1545 ALLOCATE (sto_state_overlap(nexc_search))
1546 ALLOCATE (max_overlap(natom))
1547 max_overlap = 0.0_dp
1548 ALLOCATE (state_of_mytype(natom))
1549 state_of_mytype = 0
1550 atom_of_state = 0
1551 nexc_states = 1
1552 state_of_atom = 0
1553
1554 IF (xas_control%orbital_list(1) < 0) THEN !Checks for manually selected orbitals from the localized set
1555
1556 DO istate = 1, nexc_search
1557 centers_wfn(1, istate) = localized_wfn_control%centers_set(my_spin)%array(1, istate)
1558 centers_wfn(2, istate) = localized_wfn_control%centers_set(my_spin)%array(2, istate)
1559 centers_wfn(3, istate) = localized_wfn_control%centers_set(my_spin)%array(3, istate)
1560
1561 ! Assign the state to the closest atom
1562 distmin = 100.0_dp
1563 DO iat = 1, nexc_atoms
1564 iatom = xas_control%exc_atoms(iat)
1565 ra(1:3) = particle_set(iatom)%r(1:3)
1566 rc(1:3) = centers_wfn(1:3, istate)
1567 rac = pbc(ra, rc, cell)
1568 dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
1569
1570 IF (dist < distmin) THEN
1571
1572 atom_of_state(istate) = iatom
1573 distmin = dist
1574 END IF
1575 END DO
1576 IF (atom_of_state(istate) /= 0) THEN
1577 !Character of the state
1578 CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, istate, &
1579 nao, 1, transpose=.true.)
1580
1581 iatom = atom_of_state(istate)
1582
1583 NULLIFY (atomic_kind)
1584 atomic_kind => particle_set(iatom)%atomic_kind
1585 CALL get_atomic_kind(atomic_kind=atomic_kind, &
1586 kind_number=ikind)
1587
1588 my_kind = mykind_of_kind(ikind)
1589
1590 sto_state_overlap(istate) = 0.0_dp
1591 DO i = 1, SIZE(stogto_overlap(my_kind)%array, 1)
1592 component = 0.0_dp
1593 DO j = 1, SIZE(stogto_overlap(my_kind)%array, 2)
1594 isgf = first_sgf(iatom) + j - 1
1595 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 component*component
1599 END DO
1600
1601 IF (sto_state_overlap(istate) > max_overlap(iatom)) THEN
1602 state_of_mytype(iatom) = istate
1603 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 IF (xas_control%overlap_threshold < 1) THEN
1610 DO iat = 1, nexc_atoms
1611 iatom = xas_control%exc_atoms(iat)
1612 DO istate = 1, nexc_search
1613 IF (atom_of_state(istate) == iatom) THEN
1614 IF (sto_state_overlap(istate) > max_overlap(iatom)*xas_control%overlap_threshold &
1615 .AND. istate /= state_of_mytype(iat)) THEN
1616 nexc_states(iat) = nexc_states(iat) + 1
1617 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 IF (output_unit > 0) THEN
1626 WRITE (unit=output_unit, fmt="(/,T10,A,/)") &
1627 "List the atoms to be excited and the relative of MOs index "
1628 END IF
1629
1630 DO iat = 1, nexc_atoms
1631 iatom = xas_env%exc_atoms(iat)
1632 state_of_atom(iat, 1) = state_of_mytype(iatom) ! Place the state with maximum overlap first in the list
1633 IF (output_unit > 0) THEN
1634 WRITE (unit=output_unit, fmt="(T10,A,I3,T26,A)", advance='NO') &
1635 'Atom: ', iatom, "MO index:"
1636 END IF
1637 DO istate = 1, nexc_states(iat)
1638 IF (istate < nexc_states(iat)) THEN
1639 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(I4)", advance='NO') state_of_atom(iat, istate)
1640 ELSE
1641 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(I4)") state_of_atom(iat, istate)
1642 END IF
1643 END DO
1644 IF (state_of_atom(iat, 1) == 0 .OR. state_of_atom(iat, 1) .GT. homo) THEN
1645 cpabort("A wrong state has been selected for excitation, check the Wannier centers")
1646 END IF
1647 END DO
1648
1649 IF (xas_control%overlap_threshold < 1) THEN
1650 DO iat = 1, nexc_atoms
1651 IF (output_unit > 0) THEN
1652 WRITE (unit=output_unit, fmt="(/,T10,A,I6)") &
1653 'Overlap integrals for Atom: ', iat
1654 DO istate = 1, nexc_states(iat)
1655 WRITE (unit=output_unit, fmt="(T10,A,I3,T26,A,T38,f10.8)") &
1656 '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 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 CALL reallocate(xas_env%nexc_states, 1, natom)
1668 CALL reallocate(xas_env%state_of_atom, 1, natom, 1, SIZE(xas_control%orbital_list))
1669 CALL get_xas_env(xas_env, nexc_states=nexc_states, state_of_atom=state_of_atom)
1670
1671 nexc_states = 0
1672 state_of_atom = 0
1673 nexc_atoms = natom !To include all possible atoms in the spectrum calculation
1674
1675 DO istate = 1, SIZE(xas_control%orbital_list)
1676
1677 chosen_state = xas_control%orbital_list(istate)
1678 nexc_atoms = 1
1679 centers_wfn(1, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(1, chosen_state)
1680 centers_wfn(2, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(2, chosen_state)
1681 centers_wfn(3, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(3, chosen_state)
1682
1683 distmin = 100.0_dp
1684 DO iat = 1, natom
1685 ra(1:3) = particle_set(iat)%r(1:3)
1686 rc(1:3) = centers_wfn(1:3, chosen_state)
1687 rac = pbc(ra, rc, cell)
1688 dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
1689 IF (dist < distmin) THEN
1690 atom_of_state(chosen_state) = iat !?
1691 distmin = dist
1692 END IF
1693 END DO ! iat
1694
1695 nexc_states(atom_of_state(chosen_state)) = nexc_states(atom_of_state(chosen_state)) + 1
1696 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 IF (output_unit > 0) THEN
1702 WRITE (unit=output_unit, fmt="(/,T10,A,/)") &
1703 "List the atoms to be excited and the relative of MOs index "
1704 END IF
1705
1706 DO iat = 1, natom
1707 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 'Atom: ', iat, "MO index:"
1710 DO i = 1, nexc_states(iat)
1711 IF (i < nexc_states(iat)) THEN
1712 WRITE (unit=output_unit, fmt="(I4)", advance='NO') state_of_atom(iat, i)
1713 ELSE
1714 WRITE (unit=output_unit, fmt="(I4)") state_of_atom(iat, i)
1715 END IF
1716 END DO
1717 END IF
1718 IF (state_of_atom(iat, 1) .GT. homo) THEN
1719 cpabort("A wrong state has been selected for excitation, check the Wannier centers")
1720 END IF
1721 END DO
1722
1723 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 cell%perd(1:3) = perd0(1:3)
1729
1730 DEALLOCATE (vecbuffer)
1731 DEALLOCATE (first_sgf)
1732 DEALLOCATE (sto_state_overlap)
1733 DEALLOCATE (max_overlap)
1734 DEALLOCATE (state_of_mytype)
1735
1736 END SUBROUTINE cls_assign_core_states
1737
1738END MODULE xas_methods
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:66
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition ai_contraction.F:18
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:680
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::srules
pure real(dp) function, public srules(z, ne, n, l)
...
Definition basis_set_types.F:2509
basis_set_types::deallocate_sto_basis_set
subroutine, public deallocate_sto_basis_set(sto_basis_set)
...
Definition basis_set_types.F:1964
basis_set_types::allocate_sto_basis_set
subroutine, public allocate_sto_basis_set(sto_basis_set)
...
Definition basis_set_types.F:1950
basis_set_types::create_gto_from_sto_basis
subroutine, public create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho)
...
Definition basis_set_types.F:2254
basis_set_types::set_sto_basis_set
subroutine, public set_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet)
...
Definition basis_set_types.F:2041
basis_set_types::init_orb_basis_set
subroutine, public init_orb_basis_set(gto_basis_set)
Initialise a Gaussian-type orbital (GTO) basis set data set.
Definition basis_set_types.F:1057
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Definition basis_set_types.F:615
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:744
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:118
cp_external_control
Routines to handle the external control of CP2K.
Definition cp_external_control.F:15
cp_external_control::external_control
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
Definition cp_external_control.F:90
cp_fm_pool_types
pool for for elements that are retained and released
Definition cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_create_fm
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
Definition cp_fm_pool_types.F:185
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition cp_fm_types.F:643
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:768
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:901
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::xas_1s_type
integer, parameter, public xas_1s_type
Definition input_constants.F:671
input_constants::xas_4p_type
integer, parameter, public xas_4p_type
Definition input_constants.F:671
input_constants::xas_3s_type
integer, parameter, public xas_3s_type
Definition input_constants.F:671
input_constants::do_loc_none
integer, parameter, public do_loc_none
Definition input_constants.F:452
input_constants::xas_3d_type
integer, parameter, public xas_3d_type
Definition input_constants.F:671
input_constants::xas_dscf
integer, parameter, public xas_dscf
Definition input_constants.F:663
input_constants::xas_4s_type
integer, parameter, public xas_4s_type
Definition input_constants.F:671
input_constants::xas_tp_xhh
integer, parameter, public xas_tp_xhh
Definition input_constants.F:663
input_constants::xas_2p_type
integer, parameter, public xas_2p_type
Definition input_constants.F:671
input_constants::xas_dip_len
integer, parameter, public xas_dip_len
Definition input_constants.F:682
input_constants::xas_dip_vel
integer, parameter, public xas_dip_vel
Definition input_constants.F:682
input_constants::xas_2s_type
integer, parameter, public xas_2s_type
Definition input_constants.F:671
input_constants::xas_tp_xfh
integer, parameter, public xas_tp_xfh
Definition input_constants.F:663
input_constants::xas_3p_type
integer, parameter, public xas_3p_type
Definition input_constants.F:671
input_constants::xas_tp_fh
integer, parameter, public xas_tp_fh
Definition input_constants.F:663
input_constants::xas_tp_flex
integer, parameter, public xas_tp_flex
Definition input_constants.F:663
input_constants::xas_4f_type
integer, parameter, public xas_4f_type
Definition input_constants.F:671
input_constants::state_loc_list
integer, parameter, public state_loc_list
Definition input_constants.F:467
input_constants::state_loc_range
integer, parameter, public state_loc_range
Definition input_constants.F:467
input_constants::xas_4d_type
integer, parameter, public xas_4d_type
Definition input_constants.F:671
input_constants::xas_tp_hh
integer, parameter, public xas_tp_hh
Definition input_constants.F:663
input_constants::xes_tp_val
integer, parameter, public xes_tp_val
Definition input_constants.F:663
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1005
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_diis
Apply the direct inversion in the iterative subspace (DIIS) of Pulay in the framework of an SCF itera...
Definition qs_diis.F:21
qs_diis::qs_diis_b_clear
pure subroutine, public qs_diis_b_clear(diis_buffer)
clears the buffer
Definition qs_diis.F:517
qs_diis::qs_diis_b_create
subroutine, public qs_diis_b_create(diis_buffer, nbuffer)
Allocates an SCF DIIS buffer.
Definition qs_diis.F:101
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1035
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
qs_loc_main
Driver for the localization that should be general for all the methods available and all the definiti...
Definition qs_loc_main.F:22
qs_loc_main::qs_loc_driver
subroutine, public qs_loc_driver(qs_env, qs_loc_env, print_loc_section, myspin, ext_mo_coeff)
set up the calculation of localized orbitals
Definition qs_loc_main.F:96
qs_loc_methods
Driver for the localization that should be general for all the methods available and all the definiti...
Definition qs_loc_methods.F:21
qs_loc_methods::qs_print_cubes
subroutine, public qs_print_cubes(qs_env, mo_coeff, nstates, state_list, centers, print_key, root, ispin, idir, state0, file_position)
write the cube files for a set of selected states
Definition qs_loc_methods.F:784
qs_loc_types
New version of the module for the localization of the molecular orbitals This should be able to use d...
Definition qs_loc_types.F:25
qs_loc_types::qs_loc_env_create
subroutine, public qs_loc_env_create(qs_loc_env)
...
Definition qs_loc_types.F:166
qs_loc_utils
Some utilities for the construction of the localization environment.
Definition qs_loc_utils.F:13
qs_loc_utils::set_loc_wfn_lists
subroutine, public set_loc_wfn_lists(localized_wfn_control, nmoloc, nmo, nspins, my_spin)
create the lists of mos that are taken into account
Definition qs_loc_utils.F:1626
qs_loc_utils::qs_loc_env_init
subroutine, public qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, loc_coeff, mo_loc_history)
allocates the data, and initializes the operators
Definition qs_loc_utils.F:232
qs_loc_utils::set_loc_centers
subroutine, public set_loc_centers(localized_wfn_control, nmoloc, nspins)
create the center and spread array and the file names for the output
Definition qs_loc_utils.F:1601
qs_loc_utils::qs_loc_control_init
subroutine, public qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, do_xas, nloc_xas, spin_xas)
initializes everything needed for localization of the HOMOs
Definition qs_loc_utils.F:1006
qs_matrix_pools
wrapper for the pools of matrixes
Definition qs_matrix_pools.F:14
qs_matrix_pools::mpools_get
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition qs_matrix_pools.F:141
qs_mo_io
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
qs_mo_io::write_mo_set_to_output_unit
subroutine, public write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, dft_section, before, kpoint, final_mos, spin, solver_method, rtp, cpart, sim_step, umo_set)
Write MO information to output file (eigenvalues, occupation numbers, coefficients)
Definition qs_mo_io.F:995
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::set_mo_set
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
Definition qs_mo_types.F:452
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_operators_ao
Definition qs_operators_ao.F:13
qs_operators_ao::p_xyz_ao
subroutine, public p_xyz_ao(op, qs_env, minimum_image)
Calculation of the components of the dipole operator in the velocity form The elements of the sparse ...
Definition qs_operators_ao.F:818
qs_operators_ao::rrc_xyz_ao
subroutine, public rrc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
Calculation of the components of the dipole operator in the length form by taking the relative positi...
Definition qs_operators_ao.F:1084
qs_pdos
Calculation and writing of projected density of states The DOS is computed per angular momentum and p...
Definition qs_pdos.F:15
qs_pdos::calculate_projected_dos
subroutine, public calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, dft_section, ispin, xas_mittle, external_matrix_shalf)
Compute and write projected density of states.
Definition qs_pdos.F:139
qs_scf_initialization
Utility routines for qs_scf.
Definition qs_scf_initialization.F:11
qs_scf_initialization::qs_scf_env_initialize
subroutine, public qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
initializes input parameters if needed or restores values from previous runs to fill scf_env with the...
Definition qs_scf_initialization.F:133
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
qs_scf_types::scf_env_release
subroutine, public scf_env_release(scf_env)
releases an scf_env (see doc/ReferenceCounting.html)
Definition qs_scf_types.F:226
qs_scf
Routines for the Quickstep SCF run.
Definition qs_scf.F:47
qs_scf::scf_env_cleanup
subroutine, public scf_env_cleanup(scf_env)
perform cleanup operations (like releasing temporary storage) at the end of the scf
Definition qs_scf.F:895
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
scf_control_types::scf_c_read_parameters
subroutine, public scf_c_read_parameters(scf_control, inp_section)
reads the parameters of the scf section into the given scf_control
Definition scf_control_types.F:299
scf_control_types::scf_c_create
subroutine, public scf_c_create(scf_control)
allocates and initializes an scf control object with the default values
Definition scf_control_types.F:153
xas_control
Defines control structures, which contain the parameters and the settings for the calculations.
Definition xas_control.F:12
xas_control::read_xas_control
subroutine, public read_xas_control(xas_control, xas_section)
read from input the instructions for a xes/xas calculation
Definition xas_control.F:85
xas_control::xas_control_create
subroutine, public xas_control_create(xas_control)
create retain release the xas_control_type
Definition xas_control.F:269
xas_control::write_xas_control
subroutine, public write_xas_control(xas_control, dft_section)
write on the instructions for a xes/xas calculation
Definition xas_control.F:210
xas_env_types
define create destroy get and put information in xas_env to calculate the x-ray absorption spectra
Definition xas_env_types.F:15
xas_env_types::xas_env_release
subroutine, public xas_env_release(xas_env)
...
Definition xas_env_types.F:329
xas_env_types::set_xas_env
subroutine, public set_xas_env(xas_env, nexc_search, spin_channel, nexc_atoms, nvirtual, nvirtual2, ip_energy, occ_estate, qs_loc_env, xas_estate, xas_nelectron, homo_occ, scf_env, scf_control)
...
Definition xas_env_types.F:243
xas_env_types::xas_env_create
subroutine, public xas_env_create(xas_env)
...
Definition xas_env_types.F:299
xas_env_types::get_xas_env
subroutine, public get_xas_env(xas_env, exc_state, nao, nvirtual, nvirtual2, centers_wfn, atom_of_state, exc_atoms, nexc_states, type_of_state, mykind_of_atom, mykind_of_kind, state_of_atom, spectrum, groundstate_coeff, ostrength_sm, dip_fm_set, excvec_coeff, excvec_overlap, unoccupied_orbs, unoccupied_evals, unoccupied_max_iter, unoccupied_eps, all_vectors, all_evals, my_gto_basis, qs_loc_env, stogto_overlap, occ_estate, xas_nelectron, xas_estate, nexc_atoms, nexc_search, spin_channel, scf_env, scf_control)
...
Definition xas_env_types.F:154
xas_methods
driver for the xas calculation and xas_scf for the tp method
Definition xas_methods.F:15
xas_methods::calc_stogto_overlap
subroutine, public calc_stogto_overlap(base_a, base_b, matrix)
...
Definition xas_methods.F:1402
xas_methods::xas
subroutine, public xas(qs_env, dft_control)
Driver for xas calculations The initial mos are prepared A loop on the atoms to be excited is started...
Definition xas_methods.F:162
xas_restart
Initialize the XAS orbitals for specific core excitations Either the GS orbitals are used as initial ...
Definition xas_restart.F:20
xas_restart::xas_read_restart
subroutine, public xas_read_restart(xas_env, xas_section, qs_env, xas_method, iatom, estate, istate)
Set up for reading the restart corresponding to the excitation of iatom If the corresponding restart ...
Definition xas_restart.F:104
xas_tp_scf
xas_scf for the tp method It is repeaated for every atom that have to be excited
Definition xas_tp_scf.F:15
xas_tp_scf::xas_do_tp_scf
subroutine, public xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, xas_section, scf_section, converged, should_stop)
perform an scf loop to calculate the xas spectrum given by the excitation of a inner state of a selec...
Definition xas_tp_scf.F:135
xas_tp_scf::xes_scf_once
subroutine, public xes_scf_once(qs_env, xas_env, converged, should_stop)
SCF for emission spectra calculations: vacancy in valence.
Definition xas_tp_scf.F:643
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
basis_set_types::sto_basis_set_type
Definition basis_set_types.F:100
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_array_utils::cp_2d_r_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:109
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_loc_types::localized_wfn_control_type
A type that holds controlling information for the calculation of the spread of wfn and the optimizati...
Definition qs_loc_types.F:121
qs_loc_types::qs_loc_env_type
contains all the info needed by quickstep to calculate the spread of a selected set of orbitals and i...
Definition qs_loc_types.F:80
qs_matrix_pools::qs_matrix_pools_type
container for the pools of matrixes used by qs
Definition qs_matrix_pools.F:68
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:95
scf_control_types::scf_control_type
Definition scf_control_types.F:122
xas_control::xas_control_type
A type that holds controlling information for a xas calculation.
Definition xas_control.F:40
xas_env_types::xas_environment_type
Definition xas_env_types.F:69