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