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