(git:ccc2433)
xas_tdp_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 Methods for X-Ray absorption spectroscopy (XAS) using TDDFPT
10 !> \author AB (11.2017)
11 ! **************************************************************************************************
12 
13 MODULE xas_tdp_methods
14  USE admm_types, ONLY: admm_type
15  USE admm_utils, ONLY: admm_correct_for_eigenvalues,&
16  admm_uncorrect_for_eigenvalues
17  USE atomic_kind_types, ONLY: atomic_kind_type,&
18  get_atomic_kind
19  USE basis_set_types, ONLY: &
20  allocate_sto_basis_set, create_gto_from_sto_basis, deallocate_gto_basis_set, &
21  deallocate_sto_basis_set, get_gto_basis_set, gto_basis_set_type, init_orb_basis_set, &
22  set_sto_basis_set, srules, sto_basis_set_type
23  USE bibliography, ONLY: bussy2021a,&
24  cite_reference
25  USE cell_types, ONLY: cell_type,&
26  pbc
27  USE cp_blacs_env, ONLY: cp_blacs_env_type
28  USE cp_control_types, ONLY: dft_control_type
29  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
30  copy_fm_to_dbcsr,&
31  cp_dbcsr_sm_fm_multiply
32  USE cp_files, ONLY: close_file,&
33  open_file
34  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
35  USE cp_fm_diag, ONLY: cp_fm_geeig,&
36  cp_fm_power,&
37  cp_fm_syevd
38  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
39  cp_fm_struct_release,&
40  cp_fm_struct_type
41  USE cp_fm_types, ONLY: &
42  cp_fm_copy_general, cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, &
43  cp_fm_read_unformatted, cp_fm_release, cp_fm_set_all, cp_fm_to_fm, cp_fm_to_fm_submat, &
44  cp_fm_type, cp_fm_write_unformatted
45  USE cp_log_handling, ONLY: cp_get_default_logger,&
46  cp_logger_get_default_io_unit,&
47  cp_logger_type,&
48  cp_to_string
49  USE cp_output_handling, ONLY: cp_p_file,&
50  cp_print_key_finished_output,&
51  cp_print_key_generate_filename,&
52  cp_print_key_should_output,&
53  cp_print_key_unit_nr,&
54  debug_print_level
55  USE dbcsr_api, ONLY: &
56  dbcsr_add_on_diag, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_filter, &
57  dbcsr_finalize, dbcsr_get_info, dbcsr_get_occupation, dbcsr_multiply, dbcsr_p_type, &
58  dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
59  dbcsr_type_symmetric
60  USE input_constants, ONLY: &
61  do_admm_purify_cauchy_subspace, do_admm_purify_mo_diag, do_admm_purify_none, do_loc_none, &
62  do_potential_coulomb, do_potential_id, do_potential_short, do_potential_truncated, &
63  op_loc_berry, state_loc_list, tddfpt_singlet, tddfpt_spin_cons, tddfpt_spin_flip, &
64  tddfpt_triplet, xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, xas_dip_vel, &
65  xas_not_excited, xas_tdp_by_index, xas_tdp_by_kind
66  USE input_cp2k_loc, ONLY: create_localize_section
67  USE input_section_types, ONLY: section_release,&
68  section_type,&
69  section_vals_create,&
70  section_vals_get_subs_vals,&
71  section_vals_type,&
72  section_vals_val_get,&
73  section_vals_val_set
74  USE kinds, ONLY: default_path_length,&
75  default_string_length,&
76  dp
77  USE libint_wrapper, ONLY: cp_libint_static_init
78  USE machine, ONLY: m_flush
79  USE mathlib, ONLY: get_diag
80  USE memory_utilities, ONLY: reallocate
81  USE message_passing, ONLY: mp_comm_type,&
82  mp_para_env_type
83  USE parallel_gemm_api, ONLY: parallel_gemm
84  USE parallel_rng_types, ONLY: uniform,&
85  rng_stream_type
86  USE particle_methods, ONLY: get_particle_set
87  USE particle_types, ONLY: particle_type
88  USE periodic_table, ONLY: ptable
89  USE physcon, ONLY: a_fine,&
90  angstrom,&
91  evolt
92  USE qs_environment_types, ONLY: get_qs_env,&
93  qs_environment_type
94  USE qs_interactions, ONLY: init_interaction_radii_orb_basis
95  USE qs_kind_types, ONLY: get_qs_kind,&
96  qs_kind_type
97  USE qs_loc_main, ONLY: qs_loc_driver
98  USE qs_loc_methods, ONLY: centers_spreads_berry,&
99  qs_print_cubes
100  USE qs_loc_types, ONLY: get_qs_loc_env,&
101  localized_wfn_control_create,&
102  localized_wfn_control_type,&
103  qs_loc_env_create,&
104  qs_loc_env_release,&
105  qs_loc_env_type
106  USE qs_loc_utils, ONLY: qs_loc_control_init,&
107  qs_loc_env_init,&
108  set_loc_centers
109  USE qs_mo_io, ONLY: write_mo_set_low
110  USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
111  USE qs_mo_types, ONLY: allocate_mo_set,&
112  deallocate_mo_set,&
113  duplicate_mo_set,&
114  get_mo_set,&
115  init_mo_set,&
116  mo_set_type
117  USE qs_operators_ao, ONLY: p_xyz_ao,&
118  rrc_xyz_ao
119  USE qs_pdos, ONLY: calculate_projected_dos
120  USE qs_rho_types, ONLY: qs_rho_get,&
121  qs_rho_type
122  USE qs_scf_types, ONLY: ot_method_nr
123  USE util, ONLY: get_limit,&
124  locate,&
125  sort_unique
126  USE xas_methods, ONLY: calc_stogto_overlap
127  USE xas_tdp_atom, ONLY: init_xas_atom_env,&
128  integrate_fxc_atoms,&
129  integrate_soc_atoms
130  USE xas_tdp_correction, ONLY: gw2x_shift,&
131  get_soc_splitting
132  USE xas_tdp_integrals, ONLY: compute_ri_3c_coulomb,&
133  compute_ri_3c_exchange,&
134  compute_ri_coulomb2_int,&
135  compute_ri_exchange2_int
136  USE xas_tdp_types, ONLY: &
137  donor_state_create, donor_state_type, free_ds_memory, free_exat_memory, &
138  read_xas_tdp_control, set_donor_state, set_xas_tdp_env, xas_atom_env_create, &
139  xas_atom_env_release, xas_atom_env_type, xas_tdp_control_create, xas_tdp_control_release, &
140  xas_tdp_control_type, xas_tdp_env_create, xas_tdp_env_release, xas_tdp_env_type
141  USE xas_tdp_utils, ONLY: include_os_soc,&
142  include_rcs_soc,&
143  setup_xas_tdp_prob,&
144  solve_xas_tdp_prob
145  USE xc_write_output, ONLY: xc_write
146 #include "./base/base_uses.f90"
147 
148  IMPLICIT NONE
149  PRIVATE
150 
151  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_methods'
152 
153  PUBLIC :: xas_tdp, xas_tdp_init
154 
155 CONTAINS
156 
157 ! **************************************************************************************************
158 !> \brief Driver for XAS TDDFT calculations.
159 !> \param qs_env the inherited qs_environment
160 !> \author AB
161 !> \note Empty for now...
162 ! **************************************************************************************************
163  SUBROUTINE xas_tdp(qs_env)
164 
165  TYPE(qs_environment_type), POINTER :: qs_env
166 
167  CHARACTER(len=*), PARAMETER :: routinen = 'xas_tdp'
168 
169  CHARACTER(default_string_length) :: rst_filename
170  INTEGER :: handle, n_rep, output_unit
171  LOGICAL :: do_restart
172  TYPE(section_vals_type), POINTER :: xas_tdp_section
173 
174  CALL timeset(routinen, handle)
175 
176 ! Logger initialization and XAS TDP banner printing
177  NULLIFY (xas_tdp_section)
178  xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
179  output_unit = cp_logger_get_default_io_unit()
180 
181  IF (output_unit > 0) THEN
182  WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,/,T3,A,/,T3,A,/)") &
183  "!===========================================================================!", &
184  "! XAS TDP !", &
185  "! Starting TDDFPT driven X-rays absorption spectroscopy calculations !", &
186  "!===========================================================================!"
187  END IF
188 
189  CALL cite_reference(bussy2021a)
190 
191 ! Check whether this is a restart calculation, i.e. is a restart file is provided
192  CALL section_vals_val_get(xas_tdp_section, "RESTART_FROM_FILE", n_rep_val=n_rep)
193 
194  IF (n_rep < 1) THEN
195  do_restart = .false.
196  ELSE
197  CALL section_vals_val_get(xas_tdp_section, "RESTART_FROM_FILE", c_val=rst_filename)
198  do_restart = .true.
199  END IF
200 
201 ! Restart the calculation if needed
202  IF (do_restart) THEN
203 
204  IF (output_unit > 0) THEN
205  WRITE (unit=output_unit, fmt="(/,T3,A)") &
206  "# This is a RESTART calculation for PDOS and/or CUBE printing"
207  END IF
208 
209  CALL restart_calculation(rst_filename, xas_tdp_section, qs_env)
210 
211 ! or run the core XAS_TDP routine if not
212  ELSE
213  CALL xas_tdp_core(xas_tdp_section, qs_env)
214  END IF
215 
216  IF (output_unit > 0) THEN
217  WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,/,T3,A,/)") &
218  "!===========================================================================!", &
219  "! End of TDDFPT driven X-rays absorption spectroscopy calculations !", &
220  "!===========================================================================!"
221  END IF
222 
223  CALL timestop(handle)
224 
225  END SUBROUTINE xas_tdp
226 
227 ! **************************************************************************************************
228 !> \brief The core workflow of the XAS_TDP method
229 !> \param xas_tdp_section the input values for XAS_TDP
230 !> \param qs_env ...
231 ! **************************************************************************************************
232  SUBROUTINE xas_tdp_core(xas_tdp_section, qs_env)
233 
234  TYPE(section_vals_type), POINTER :: xas_tdp_section
235  TYPE(qs_environment_type), POINTER :: qs_env
236 
237  CHARACTER(LEN=default_string_length) :: kind_name
238  INTEGER :: batch_size, bo(2), current_state_index, iat, iatom, ibatch, ikind, ispin, istate, &
239  nbatch, nex_atom, output_unit, tmp_index
240  INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_atoms, ex_atoms_of_kind
241  INTEGER, DIMENSION(:), POINTER :: atoms_of_kind
242  LOGICAL :: do_os, end_of_batch, unique
243  TYPE(admm_type), POINTER :: admm_env
244  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
245  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
246  TYPE(dft_control_type), POINTER :: dft_control
247  TYPE(donor_state_type), POINTER :: current_state
248  TYPE(gto_basis_set_type), POINTER :: tmp_basis
249  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
250  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
251  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
252  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
253 
254  NULLIFY (xas_tdp_env, xas_tdp_control, atomic_kind_set, atoms_of_kind, current_state)
255  NULLIFY (xas_atom_env, dft_control, matrix_ks, admm_env, qs_kind_set, tmp_basis)
256 
257 ! Initialization
258  output_unit = cp_logger_get_default_io_unit()
259 
260  IF (output_unit > 0) THEN
261  WRITE (unit=output_unit, fmt="(/,T3,A)") &
262  "# Create and initialize the XAS_TDP environment"
263  END IF
264  CALL get_qs_env(qs_env, dft_control=dft_control)
265  CALL xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
266  CALL print_info(output_unit, xas_tdp_control, qs_env)
267 
268  IF (output_unit > 0) THEN
269  IF (xas_tdp_control%check_only) THEN
270  cpwarn("This is a CHECK_ONLY run for donor MOs verification")
271  END IF
272  END IF
273 
274 ! Localization of the core orbitals if requested (used for better identification of donor states)
275  IF (xas_tdp_control%do_loc) THEN
276  IF (output_unit > 0) THEN
277  WRITE (unit=output_unit, fmt="(/,T3,A,/)") &
278  "# Localizing core orbitals for better identification"
279  END IF
280 ! closed shell or ROKS => myspin=1
281  IF (xas_tdp_control%do_uks) THEN
282  DO ispin = 1, dft_control%nspins
283  CALL qs_loc_driver(qs_env, xas_tdp_env%qs_loc_env, &
284  xas_tdp_control%print_loc_subsection, myspin=ispin)
285  END DO
286  ELSE
287  CALL qs_loc_driver(qs_env, xas_tdp_env%qs_loc_env, &
288  xas_tdp_control%print_loc_subsection, myspin=1)
289  END IF
290  END IF
291 
292 ! Find the MO centers
293  CALL find_mo_centers(xas_tdp_env, xas_tdp_control, qs_env)
294 
295 ! Assign lowest energy orbitals to excited atoms
296  CALL assign_mos_to_ex_atoms(xas_tdp_env, xas_tdp_control, qs_env)
297 
298 ! Once assigned, diagonalize the MOs wrt the KS matrix in the subspace associated to each atom
299  IF (xas_tdp_control%do_loc) THEN
300  IF (output_unit > 0) THEN
301  WRITE (unit=output_unit, fmt="(/,T3,A,/,T5,A)") &
302  "# Diagonalize localized MOs wrt the KS matrix in the subspace of each excited", &
303  "atom for better donor state identification."
304  END IF
305  CALL diagonalize_assigned_mo_subset(xas_tdp_env, xas_tdp_control, qs_env)
306  ! update MO centers
307  CALL find_mo_centers(xas_tdp_env, xas_tdp_control, qs_env)
308  END IF
309 
310  IF (output_unit > 0) THEN
311  WRITE (unit=output_unit, fmt="(/,T3,A,I4,A,/)") &
312  "# Assign the relevant subset of the ", xas_tdp_control%n_search, &
313  " lowest energy MOs to excited atoms"
314  END IF
315  CALL write_mos_to_ex_atoms_association(xas_tdp_env, xas_tdp_control, qs_env)
316 
317 ! If CHECK_ONLY run, check the donor MOs
318  IF (xas_tdp_control%check_only) CALL print_checks(xas_tdp_env, xas_tdp_control, qs_env)
319 
320 ! If not simply exact exchange, setup a xas_atom_env and compute the xc integrals on the atomic grids
321 ! Also needed if SOC is included or XPS GW2X(). Done before looping on atoms as it's all done at once
322  IF ((xas_tdp_control%do_xc .OR. xas_tdp_control%do_soc .OR. xas_tdp_control%do_gw2x) &
323  .AND. .NOT. xas_tdp_control%check_only) THEN
324 
325  IF (output_unit > 0 .AND. xas_tdp_control%do_xc) THEN
326  WRITE (unit=output_unit, fmt="(/,T3,A,I4,A)") &
327  "# Integrating the xc kernel on the atomic grids ..."
328  CALL m_flush(output_unit)
329  END IF
330 
331  CALL xas_atom_env_create(xas_atom_env)
332  CALL init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env)
333  do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
334 
335  IF (xas_tdp_control%do_xc .AND. (.NOT. xas_tdp_control%xps_only)) THEN
336  CALL integrate_fxc_atoms(xas_tdp_env%ri_fxc, xas_atom_env, xas_tdp_control, qs_env)
337  END IF
338 
339  IF (xas_tdp_control%do_soc .OR. xas_tdp_control%do_gw2x) THEN
340  CALL integrate_soc_atoms(xas_tdp_env%orb_soc, xas_atom_env=xas_atom_env, qs_env=qs_env)
341  END IF
342 
343  CALL xas_atom_env_release(xas_atom_env)
344  END IF
345 
346 ! Compute the 3-center Coulomb integrals for the whole system
347  IF ((.NOT. (xas_tdp_control%check_only .OR. xas_tdp_control%xps_only)) .AND. &
348  (xas_tdp_control%do_xc .OR. xas_tdp_control%do_coulomb)) THEN
349  IF (output_unit > 0) THEN
350  WRITE (unit=output_unit, fmt="(/,T3,A,I4,A)") &
351  "# Computing the RI 3-center Coulomb integrals ..."
352  CALL m_flush(output_unit)
353  END IF
354  CALL compute_ri_3c_coulomb(xas_tdp_env, qs_env)
355 
356  END IF
357 
358 ! Loop over donor states for calculation
359  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
360  current_state_index = 1
361 
362 ! Loop over atomic kinds
363  DO ikind = 1, SIZE(atomic_kind_set)
364 
365  IF (xas_tdp_control%check_only) EXIT
366  IF (.NOT. any(xas_tdp_env%ex_kind_indices == ikind)) cycle
367 
368  CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), name=kind_name, &
369  atom_list=atoms_of_kind)
370 
371  ! compute the RI coulomb2 inverse for this kind, and RI exchange2 if needed
372  CALL compute_ri_coulomb2_int(ikind, xas_tdp_env, xas_tdp_control, qs_env)
373  IF (xas_tdp_control%do_hfx) THEN
374  CALL compute_ri_exchange2_int(ikind, xas_tdp_env, xas_tdp_control, qs_env)
375  END IF
376 
377  !Randomly distribute excited atoms of current kinds into batches for optimal load balance
378  !of RI 3c exchange integrals. Take batch sizes of 2 to avoid taxing memory too much, while
379  !greatly improving load balance
380  batch_size = 2
381  CALL get_ri_3c_batches(ex_atoms_of_kind, nbatch, batch_size, atoms_of_kind, xas_tdp_env)
382  nex_atom = SIZE(ex_atoms_of_kind)
383 
384  !Loop over batches
385  DO ibatch = 1, nbatch
386 
387  bo = get_limit(nex_atom, nbatch, ibatch - 1)
388  batch_size = bo(2) - bo(1) + 1
389  ALLOCATE (batch_atoms(batch_size))
390  iatom = 0
391  DO iat = bo(1), bo(2)
392  iatom = iatom + 1
393  batch_atoms(iatom) = ex_atoms_of_kind(iat)
394  END DO
395  CALL sort_unique(batch_atoms, unique)
396 
397  !compute RI 3c exchange integrals on batch, if so required
398  IF (xas_tdp_control%do_hfx) THEN
399  IF (output_unit > 0) THEN
400  WRITE (unit=output_unit, fmt="(/,T3,A,I4,A,I4,A,I1,A,A)") &
401  "# Computing the RI 3-center Exchange integrals for batch ", ibatch, "(/", nbatch, ") of ", &
402  batch_size, " atoms of kind: ", trim(kind_name)
403  CALL m_flush(output_unit)
404  END IF
405  CALL compute_ri_3c_exchange(batch_atoms, xas_tdp_env, xas_tdp_control, qs_env)
406  END IF
407 
408 ! Loop over atoms of batch
409  DO iat = 1, batch_size
410  iatom = batch_atoms(iat)
411 
412  tmp_index = locate(xas_tdp_env%ex_atom_indices, iatom)
413 
414  !if dipole in length rep, compute the dipole in the AO basis for this atom
415  !if quadrupole is required, compute it there too (in length rep)
416  IF (xas_tdp_control%dipole_form == xas_dip_len .OR. xas_tdp_control%do_quad) THEN
417  CALL compute_lenrep_multipole(iatom, xas_tdp_env, xas_tdp_control, qs_env)
418  END IF
419 
420 ! Loop over states of excited atom of kind
421  DO istate = 1, SIZE(xas_tdp_env%state_types, 1)
422 
423  IF (xas_tdp_env%state_types(istate, tmp_index) == xas_not_excited) cycle
424 
425  current_state => xas_tdp_env%donor_states(current_state_index)
426  CALL set_donor_state(current_state, at_index=iatom, &
427  at_symbol=kind_name, kind_index=ikind, &
428  state_type=xas_tdp_env%state_types(istate, tmp_index))
429 
430 ! Initial write for the donor state of interest
431  IF (output_unit > 0) THEN
432  WRITE (unit=output_unit, fmt="(/,T3,A,A2,A,I4,A,A,/)") &
433  "# Start of calculations for donor state of type ", &
434  xas_tdp_env%state_type_char(current_state%state_type), " for atom", &
435  current_state%at_index, " of kind ", trim(current_state%at_symbol)
436  CALL m_flush(output_unit)
437  END IF
438 
439 ! Assign best fitting MO(s) to current core donnor state
440  CALL assign_mos_to_donor_state(current_state, xas_tdp_env, xas_tdp_control, qs_env)
441 
442 ! Perform MO restricted Mulliken pop analysis for verification
443  CALL perform_mulliken_on_donor_state(current_state, qs_env)
444 
445 ! GW2X correction
446  IF (xas_tdp_control%do_gw2x) THEN
447  CALL gw2x_shift(current_state, xas_tdp_env, xas_tdp_control, qs_env)
448  END IF
449 
450 ! Do main XAS calculations here
451  IF (.NOT. xas_tdp_control%xps_only) THEN
452  CALL setup_xas_tdp_prob(current_state, qs_env, xas_tdp_env, xas_tdp_control)
453 
454  IF (xas_tdp_control%do_spin_cons) THEN
455  CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
456  ex_type=tddfpt_spin_cons)
457  CALL compute_dipole_fosc(current_state, xas_tdp_control, xas_tdp_env)
458  IF (xas_tdp_control%do_quad) CALL compute_quadrupole_fosc(current_state, &
459  xas_tdp_control, xas_tdp_env)
460  CALL xas_tdp_post(tddfpt_spin_cons, current_state, xas_tdp_env, &
461  xas_tdp_section, qs_env)
462  CALL write_donor_state_restart(tddfpt_spin_cons, current_state, xas_tdp_section, qs_env)
463  END IF
464 
465  IF (xas_tdp_control%do_spin_flip) THEN
466  CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
467  ex_type=tddfpt_spin_flip)
468  !no dipole in spin-flip (spin-forbidden)
469  CALL xas_tdp_post(tddfpt_spin_flip, current_state, xas_tdp_env, &
470  xas_tdp_section, qs_env)
471  CALL write_donor_state_restart(tddfpt_spin_flip, current_state, xas_tdp_section, qs_env)
472  END IF
473 
474  IF (xas_tdp_control%do_singlet) THEN
475  CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
476  ex_type=tddfpt_singlet)
477  CALL compute_dipole_fosc(current_state, xas_tdp_control, xas_tdp_env)
478  IF (xas_tdp_control%do_quad) CALL compute_quadrupole_fosc(current_state, &
479  xas_tdp_control, xas_tdp_env)
480  CALL xas_tdp_post(tddfpt_singlet, current_state, xas_tdp_env, &
481  xas_tdp_section, qs_env)
482  CALL write_donor_state_restart(tddfpt_singlet, current_state, xas_tdp_section, qs_env)
483  END IF
484 
485  IF (xas_tdp_control%do_triplet) THEN
486  CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
487  ex_type=tddfpt_triplet)
488  !no dipole for triplets by construction
489  CALL xas_tdp_post(tddfpt_triplet, current_state, xas_tdp_env, &
490  xas_tdp_section, qs_env)
491  CALL write_donor_state_restart(tddfpt_triplet, current_state, xas_tdp_section, qs_env)
492  END IF
493 
494 ! Include the SOC if required, only for 2p donor stataes
495  IF (xas_tdp_control%do_soc .AND. current_state%state_type == xas_2p_type) THEN
496  IF (xas_tdp_control%do_singlet .AND. xas_tdp_control%do_triplet) THEN
497  CALL include_rcs_soc(current_state, xas_tdp_env, xas_tdp_control, qs_env)
498  END IF
499  IF (xas_tdp_control%do_spin_cons .AND. xas_tdp_control%do_spin_flip) THEN
500  CALL include_os_soc(current_state, xas_tdp_env, xas_tdp_control, qs_env)
501  END IF
502  END IF
503 
504 ! Print the requested properties
505  CALL print_xas_tdp_to_file(current_state, xas_tdp_env, xas_tdp_control, xas_tdp_section)
506  END IF !xps_only
507  IF (xas_tdp_control%do_gw2x) CALL print_xps(current_state, xas_tdp_env, xas_tdp_control, qs_env)
508 
509 ! Free some unneeded attributes of current_state
510  CALL free_ds_memory(current_state)
511 
512  current_state_index = current_state_index + 1
513  NULLIFY (current_state)
514 
515  END DO ! state type
516 
517  end_of_batch = .false.
518  IF (iat == batch_size) end_of_batch = .true.
519  CALL free_exat_memory(xas_tdp_env, iatom, end_of_batch)
520  END DO ! atom of batch
521  DEALLOCATE (batch_atoms)
522  END DO !ibatch
523  DEALLOCATE (ex_atoms_of_kind)
524  END DO ! kind
525 
526 ! Return to ususal KS matrix
527  IF (dft_control%do_admm) THEN
528  CALL get_qs_env(qs_env, matrix_ks=matrix_ks, admm_env=admm_env)
529  DO ispin = 1, dft_control%nspins
530  CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks(ispin)%matrix)
531  END DO
532  END IF
533 
534 ! Return to initial basis set radii
535  IF (xas_tdp_control%eps_pgf > 0.0_dp) THEN
536  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
537  DO ikind = 1, SIZE(atomic_kind_set)
538  CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="ORB")
539  CALL init_interaction_radii_orb_basis(tmp_basis, eps_pgf_orb=dft_control%qs_control%eps_pgf_orb)
540  END DO
541  END IF
542 
543 ! Clean-up
544  CALL xas_tdp_env_release(xas_tdp_env)
545  CALL xas_tdp_control_release(xas_tdp_control)
546 
547  END SUBROUTINE xas_tdp_core
548 
549 ! **************************************************************************************************
550 !> \brief Overall control and environment types initialization
551 !> \param xas_tdp_env the environment type to initialize
552 !> \param xas_tdp_control the control type to initialize
553 !> \param qs_env the inherited qs environement type
554 ! **************************************************************************************************
555  SUBROUTINE xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
556 
557  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
558  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
559  TYPE(qs_environment_type), POINTER :: qs_env
560 
561  CHARACTER(LEN=default_string_length) :: kind_name
562  INTEGER :: at_ind, i, ispin, j, k, kind_ind, &
563  n_donor_states, n_kinds, nao, &
564  nat_of_kind, natom, nex_atoms, &
565  nex_kinds, nmatch, nspins
566  INTEGER, DIMENSION(2) :: homo, n_mo, n_moloc
567  INTEGER, DIMENSION(:), POINTER :: ind_of_kind
568  LOGICAL :: do_os, do_uks, unique
569  REAL(dp) :: fact
570  REAL(dp), DIMENSION(:), POINTER :: mo_evals
571  TYPE(admm_type), POINTER :: admm_env
572  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: at_kind_set
573  TYPE(cell_type), POINTER :: cell
574  TYPE(cp_fm_type), POINTER :: mo_coeff
575  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
576  TYPE(dbcsr_type) :: matrix_tmp
577  TYPE(dft_control_type), POINTER :: dft_control
578  TYPE(gto_basis_set_type), POINTER :: tmp_basis
579  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
580  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
581  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
582  TYPE(qs_loc_env_type), POINTER :: qs_loc_env
583  TYPE(qs_rho_type), POINTER :: rho
584  TYPE(section_type), POINTER :: dummy_section
585  TYPE(section_vals_type), POINTER :: loc_section, xas_tdp_section
586 
587  NULLIFY (xas_tdp_section, at_kind_set, ind_of_kind, dft_control, qs_kind_set, tmp_basis)
588  NULLIFY (qs_loc_env, loc_section, mos, particle_set, rho, rho_ao, mo_evals, cell)
589  NULLIFY (mo_coeff, matrix_ks, admm_env, dummy_section)
590 
591 ! XAS TDP control type initialization
592  xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
593 
594  CALL get_qs_env(qs_env, dft_control=dft_control)
595  CALL xas_tdp_control_create(xas_tdp_control)
596  CALL read_xas_tdp_control(xas_tdp_control, xas_tdp_section)
597 
598 ! Check the qs_env for a LSD/ROKS calculation
599  IF (dft_control%uks) xas_tdp_control%do_uks = .true.
600  IF (dft_control%roks) xas_tdp_control%do_roks = .true.
601  do_uks = xas_tdp_control%do_uks
602  do_os = do_uks .OR. xas_tdp_control%do_roks
603 
604 ! XAS TDP environment type initialization
605  CALL xas_tdp_env_create(xas_tdp_env)
606 
607 ! Retrieving the excited atoms indices and correspondig state types
608  IF (xas_tdp_control%define_excited == xas_tdp_by_index) THEN
609 
610 ! simply copy indices from xas_tdp_control
611  nex_atoms = SIZE(xas_tdp_control%list_ex_atoms)
612  CALL set_xas_tdp_env(xas_tdp_env, nex_atoms=nex_atoms)
613  ALLOCATE (xas_tdp_env%ex_atom_indices(nex_atoms))
614  ALLOCATE (xas_tdp_env%state_types(SIZE(xas_tdp_control%state_types, 1), nex_atoms))
615  xas_tdp_env%ex_atom_indices = xas_tdp_control%list_ex_atoms
616  xas_tdp_env%state_types = xas_tdp_control%state_types
617 
618 ! Test that these indices are within the range of available atoms
619  CALL get_qs_env(qs_env=qs_env, natom=natom)
620  IF (any(xas_tdp_env%ex_atom_indices > natom)) THEN
621  cpabort("Invalid index for the ATOM_LIST keyword.")
622  END IF
623 
624 ! Check atom kinds and fill corresponding array
625  ALLOCATE (xas_tdp_env%ex_kind_indices(nex_atoms))
626  xas_tdp_env%ex_kind_indices = 0
627  k = 0
628  CALL get_qs_env(qs_env, particle_set=particle_set)
629  DO i = 1, nex_atoms
630  at_ind = xas_tdp_env%ex_atom_indices(i)
631  CALL get_atomic_kind(particle_set(at_ind)%atomic_kind, kind_number=j)
632  IF (all(abs(xas_tdp_env%ex_kind_indices - j) .NE. 0)) THEN
633  k = k + 1
634  xas_tdp_env%ex_kind_indices(k) = j
635  END IF
636  END DO
637  nex_kinds = k
638  CALL set_xas_tdp_env(xas_tdp_env, nex_kinds=nex_kinds)
639  CALL reallocate(xas_tdp_env%ex_kind_indices, 1, nex_kinds)
640 
641  ELSE IF (xas_tdp_control%define_excited == xas_tdp_by_kind) THEN
642 
643 ! need to find out which atom of which kind is excited
644  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=at_kind_set)
645  n_kinds = SIZE(at_kind_set)
646  nex_atoms = 0
647 
648  nex_kinds = SIZE(xas_tdp_control%list_ex_kinds)
649  ALLOCATE (xas_tdp_env%ex_kind_indices(nex_kinds))
650  k = 0
651 
652  DO i = 1, n_kinds
653  CALL get_atomic_kind(atomic_kind=at_kind_set(i), name=kind_name, &
654  natom=nat_of_kind, kind_number=kind_ind)
655  IF (any(xas_tdp_control%list_ex_kinds == kind_name)) THEN
656  nex_atoms = nex_atoms + nat_of_kind
657  k = k + 1
658  xas_tdp_env%ex_kind_indices(k) = kind_ind
659  END IF
660  END DO
661 
662  ALLOCATE (xas_tdp_env%ex_atom_indices(nex_atoms))
663  ALLOCATE (xas_tdp_env%state_types(SIZE(xas_tdp_control%state_types, 1), nex_atoms))
664  nex_atoms = 0
665  nmatch = 0
666 
667  DO i = 1, n_kinds
668  CALL get_atomic_kind(atomic_kind=at_kind_set(i), name=kind_name, &
669  natom=nat_of_kind, atom_list=ind_of_kind)
670  DO j = 1, nex_kinds
671  IF (xas_tdp_control%list_ex_kinds(j) == kind_name) THEN
672  xas_tdp_env%ex_atom_indices(nex_atoms + 1:nex_atoms + nat_of_kind) = ind_of_kind
673  DO k = 1, SIZE(xas_tdp_control%state_types, 1)
674  xas_tdp_env%state_types(k, nex_atoms + 1:nex_atoms + nat_of_kind) = &
675  xas_tdp_control%state_types(k, j)
676  END DO
677  nex_atoms = nex_atoms + nat_of_kind
678  nmatch = nmatch + 1
679  END IF
680  END DO
681  END DO
682 
683  CALL set_xas_tdp_env(xas_tdp_env, nex_atoms=nex_atoms, nex_kinds=nex_kinds)
684 
685 ! Verifying that the input was valid
686  IF (nmatch .NE. SIZE(xas_tdp_control%list_ex_kinds)) THEN
687  cpabort("Invalid kind(s) for the KIND_LIST keyword.")
688  END IF
689 
690  END IF
691 
692 ! Sort the excited atoms indices (for convinience and use of locate function)
693  CALL sort_unique(xas_tdp_env%ex_atom_indices, unique)
694  IF (.NOT. unique) THEN
695  cpabort("Excited atoms not uniquely defined.")
696  END IF
697 
698 ! Check for periodicity
699  CALL get_qs_env(qs_env, cell=cell)
700  IF (all(cell%perd == 0)) THEN
701  xas_tdp_control%is_periodic = .false.
702  ELSE IF (all(cell%perd == 1)) THEN
703  xas_tdp_control%is_periodic = .true.
704  ELSE
705  cpabort("XAS TDP only implemented for full PBCs or non-PBCs")
706  END IF
707 
708 ! Allocating memory for the array of donor states
709  n_donor_states = count(xas_tdp_env%state_types /= xas_not_excited)
710  ALLOCATE (xas_tdp_env%donor_states(n_donor_states))
711  DO i = 1, n_donor_states
712  CALL donor_state_create(xas_tdp_env%donor_states(i))
713  END DO
714 
715 ! In case of ADMM, for the whole duration of the XAS_TDP, we need the total KS matrix
716  IF (dft_control%do_admm) THEN
717  CALL get_qs_env(qs_env, admm_env=admm_env, matrix_ks=matrix_ks)
718 
719  DO ispin = 1, dft_control%nspins
720  CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks(ispin)%matrix)
721  END DO
722  END IF
723 
724 ! In case of externally imposed EPS_PGF_XAS, need to update ORB and RI_XAS interaction radii
725  IF (xas_tdp_control%eps_pgf > 0.0_dp) THEN
726  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
727 
728  DO i = 1, SIZE(qs_kind_set)
729  CALL get_qs_kind(qs_kind_set(i), basis_set=tmp_basis, basis_type="ORB")
730  CALL init_interaction_radii_orb_basis(tmp_basis, eps_pgf_orb=xas_tdp_control%eps_pgf)
731  CALL get_qs_kind(qs_kind_set(i), basis_set=tmp_basis, basis_type="RI_XAS")
732  CALL init_interaction_radii_orb_basis(tmp_basis, eps_pgf_orb=xas_tdp_control%eps_pgf)
733  END DO
734  END IF
735 
736 ! In case of ground state OT optimization, compute the MO eigenvalues and get canonical MOs
737  IF (qs_env%scf_env%method == ot_method_nr) THEN
738 
739  CALL get_qs_env(qs_env, mos=mos, matrix_ks=matrix_ks)
740  nspins = 1; IF (do_uks) nspins = 2
741 
742  DO ispin = 1, nspins
743  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_evals)
744  CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, evals_arg=mo_evals)
745  END DO
746  END IF
747 
748 ! Initializing the qs_loc_env from the LOCALIZE subsection of XAS_TDP (largely inpired by MI's XAS)
749 ! We create the LOCALIZE subsection here, since it is completely overwritten anyways
750  CALL create_localize_section(dummy_section)
751  CALL section_vals_create(xas_tdp_control%loc_subsection, dummy_section)
752  CALL section_vals_val_set(xas_tdp_control%loc_subsection, "_SECTION_PARAMETERS_", &
753  l_val=xas_tdp_control%do_loc)
754  CALL section_release(dummy_section)
755  xas_tdp_control%print_loc_subsection => section_vals_get_subs_vals( &
756  xas_tdp_control%loc_subsection, "PRINT")
757 
758  ALLOCATE (xas_tdp_env%qs_loc_env)
759  CALL qs_loc_env_create(xas_tdp_env%qs_loc_env)
760  qs_loc_env => xas_tdp_env%qs_loc_env
761  loc_section => xas_tdp_control%loc_subsection
762 ! getting the number of MOs
763  CALL get_qs_env(qs_env, mos=mos)
764  CALL get_mo_set(mos(1), nmo=n_mo(1), homo=homo(1), nao=nao)
765  n_mo(2) = n_mo(1)
766  homo(2) = homo(1)
767  nspins = 1
768  IF (do_os) CALL get_mo_set(mos(2), nmo=n_mo(2), homo=homo(2))
769  IF (do_uks) nspins = 2 !in roks, same MOs for both spins
770 
771  ! by default, all (doubly occupied) homo are localized
772  IF (xas_tdp_control%n_search < 0 .OR. xas_tdp_control%n_search > minval(homo)) &
773  xas_tdp_control%n_search = minval(homo)
774  CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.true., do_xas=.true., &
775  nloc_xas=xas_tdp_control%n_search, spin_xas=1)
776 
777  ! do_xas argument above only prepares spin-alpha localization
778  IF (do_uks) THEN
779  qs_loc_env%localized_wfn_control%nloc_states(2) = xas_tdp_control%n_search
780  qs_loc_env%localized_wfn_control%lu_bound_states(1, 2) = 1
781  qs_loc_env%localized_wfn_control%lu_bound_states(2, 2) = xas_tdp_control%n_search
782  END IF
783 
784 ! final qs_loc_env initialization. Impose Berry operator
785  qs_loc_env%localized_wfn_control%operator_type = op_loc_berry
786  qs_loc_env%localized_wfn_control%max_iter = 25000
787  IF (.NOT. xas_tdp_control%do_loc) THEN
788  qs_loc_env%localized_wfn_control%localization_method = do_loc_none
789  ELSE
790  n_moloc = qs_loc_env%localized_wfn_control%nloc_states
791  CALL set_loc_centers(qs_loc_env%localized_wfn_control, n_moloc, nspins)
792  IF (do_uks) THEN
793  CALL qs_loc_env_init(qs_loc_env, qs_loc_env%localized_wfn_control, &
794  qs_env, do_localize=.true.)
795  ELSE
796  CALL qs_loc_env_init(qs_loc_env, qs_loc_env%localized_wfn_control, &
797  qs_env, do_localize=.true., myspin=1)
798  END IF
799  END IF
800 
801 ! Allocating memory for the array of excited atoms MOs. Worst case senario, all searched MOs are
802 ! associated to the same atom
803  ALLOCATE (xas_tdp_env%mos_of_ex_atoms(xas_tdp_control%n_search, nex_atoms, nspins))
804 
805 ! Compute the projector on the unoccupied, unperturbed ground state: Q = 1 - SP, sor each spin
806  IF (do_os) nspins = 2
807  CALL get_qs_env(qs_env, rho=rho, matrix_s=matrix_s)
808  CALL qs_rho_get(rho, rho_ao=rho_ao)
809 
810  ALLOCATE (xas_tdp_env%q_projector(nspins))
811  ALLOCATE (xas_tdp_env%q_projector(1)%matrix)
812  CALL dbcsr_create(xas_tdp_env%q_projector(1)%matrix, name="Q PROJECTOR ALPHA", &
813  template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
814  IF (do_os) THEN
815  ALLOCATE (xas_tdp_env%q_projector(2)%matrix)
816  CALL dbcsr_create(xas_tdp_env%q_projector(2)%matrix, name="Q PROJECTOR BETA", &
817  template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
818  END IF
819 
820 ! In the case of spin-restricted calculations, rho_ao includes double occupency => 0.5 prefactor
821  fact = -0.5_dp; IF (do_os) fact = -1.0_dp
822  CALL dbcsr_multiply('N', 'N', fact, matrix_s(1)%matrix, rho_ao(1)%matrix, 0.0_dp, &
823  xas_tdp_env%q_projector(1)%matrix, filter_eps=xas_tdp_control%eps_filter)
824  CALL dbcsr_add_on_diag(xas_tdp_env%q_projector(1)%matrix, 1.0_dp)
825  CALL dbcsr_finalize(xas_tdp_env%q_projector(1)%matrix)
826 
827  IF (do_os) THEN
828  CALL dbcsr_multiply('N', 'N', fact, matrix_s(1)%matrix, rho_ao(2)%matrix, 0.0_dp, &
829  xas_tdp_env%q_projector(2)%matrix, filter_eps=xas_tdp_control%eps_filter)
830  CALL dbcsr_add_on_diag(xas_tdp_env%q_projector(2)%matrix, 1.0_dp)
831  CALL dbcsr_finalize(xas_tdp_env%q_projector(2)%matrix)
832  END IF
833 
834 ! Create the structure for the dipole in the AO basis
835  ALLOCATE (xas_tdp_env%dipmat(3))
836  DO i = 1, 3
837  ALLOCATE (xas_tdp_env%dipmat(i)%matrix)
838  CALL dbcsr_copy(matrix_tmp, matrix_s(1)%matrix, name="XAS TDP dipole matrix")
839  IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
840  CALL dbcsr_create(xas_tdp_env%dipmat(i)%matrix, template=matrix_s(1)%matrix, &
841  matrix_type=dbcsr_type_antisymmetric)
842  CALL dbcsr_complete_redistribute(matrix_tmp, xas_tdp_env%dipmat(i)%matrix)
843  ELSE
844  CALL dbcsr_create(xas_tdp_env%dipmat(i)%matrix, template=matrix_s(1)%matrix, &
845  matrix_type=dbcsr_type_symmetric)
846  CALL dbcsr_copy(xas_tdp_env%dipmat(i)%matrix, matrix_tmp)
847  END IF
848  CALL dbcsr_set(xas_tdp_env%dipmat(i)%matrix, 0.0_dp)
849  CALL dbcsr_release(matrix_tmp)
850  END DO
851 
852 ! Create the structure for the electric quadrupole in the AO basis, if required
853  IF (xas_tdp_control%do_quad) THEN
854  ALLOCATE (xas_tdp_env%quadmat(6))
855  DO i = 1, 6
856  ALLOCATE (xas_tdp_env%quadmat(i)%matrix)
857  CALL dbcsr_copy(xas_tdp_env%quadmat(i)%matrix, matrix_s(1)%matrix, name="XAS TDP quadrupole matrix")
858  CALL dbcsr_set(xas_tdp_env%quadmat(i)%matrix, 0.0_dp)
859  END DO
860  END IF
861 
862 ! Precompute it in the velocity representation, if so chosen
863  IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
864  !enforce minimum image to avoid any PBCs related issues. Ok because very localized densities
865  CALL p_xyz_ao(xas_tdp_env%dipmat, qs_env, minimum_image=.true.)
866  END IF
867 
868 ! Allocate SOC in AO basis matrices
869  IF (xas_tdp_control%do_soc .OR. xas_tdp_control%do_gw2x) THEN
870  ALLOCATE (xas_tdp_env%orb_soc(3))
871  DO i = 1, 3
872  ALLOCATE (xas_tdp_env%orb_soc(i)%matrix)
873  END DO
874  END IF
875 
876 ! Check that everything is allowed
877  CALL safety_check(xas_tdp_control, qs_env)
878 
879 ! Initialize libint for the 3-center integrals
880  CALL cp_libint_static_init()
881 
882 ! Compute LUMOs as guess for OT solver and/or for GW2X correction
883  IF (xas_tdp_control%do_ot .OR. xas_tdp_control%do_gw2x) THEN
884  CALL make_lumo_guess(xas_tdp_env, xas_tdp_control, qs_env)
885  END IF
886 
887  END SUBROUTINE xas_tdp_init
888 
889 ! **************************************************************************************************
890 !> \brief splits the excited atoms of a kind into batches for RI 3c integrals load balance
891 !> \param ex_atoms_of_kind the excited atoms for the current kind, randomly shuffled
892 !> \param nbatch number of batches to loop over
893 !> \param batch_size standard size of a batch
894 !> \param atoms_of_kind number of atoms for the current kind (excited or not)
895 !> \param xas_tdp_env ...
896 ! **************************************************************************************************
897  SUBROUTINE get_ri_3c_batches(ex_atoms_of_kind, nbatch, batch_size, atoms_of_kind, xas_tdp_env)
898 
899  INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: ex_atoms_of_kind
900  INTEGER, INTENT(OUT) :: nbatch
901  INTEGER, INTENT(IN) :: batch_size
902  INTEGER, DIMENSION(:), INTENT(IN) :: atoms_of_kind
903  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
904 
905  INTEGER :: iat, iatom, nex_atom
906  TYPE(rng_stream_type), ALLOCATABLE :: rng_stream
907 
908  !Get the atoms from atoms_of_kind that are excited
909  nex_atom = 0
910  DO iat = 1, SIZE(atoms_of_kind)
911  iatom = atoms_of_kind(iat)
912  IF (.NOT. any(xas_tdp_env%ex_atom_indices == iatom)) cycle
913  nex_atom = nex_atom + 1
914  END DO
915 
916  ALLOCATE (ex_atoms_of_kind(nex_atom))
917  nex_atom = 0
918  DO iat = 1, SIZE(atoms_of_kind)
919  iatom = atoms_of_kind(iat)
920  IF (.NOT. any(xas_tdp_env%ex_atom_indices == iatom)) cycle
921  nex_atom = nex_atom + 1
922  ex_atoms_of_kind(nex_atom) = iatom
923  END DO
924 
925  !We shuffle those atoms to spread them
926  rng_stream = rng_stream_type(name="uniform_rng", distribution_type=uniform)
927  CALL rng_stream%shuffle(ex_atoms_of_kind(1:nex_atom))
928 
929  nbatch = nex_atom/batch_size
930  IF (nbatch*batch_size .NE. nex_atom) nbatch = nbatch + 1
931 
932  END SUBROUTINE get_ri_3c_batches
933 
934 ! **************************************************************************************************
935 !> \brief Checks for forbidden keywords combinations
936 !> \param xas_tdp_control ...
937 !> \param qs_env ...
938 ! **************************************************************************************************
939  SUBROUTINE safety_check(xas_tdp_control, qs_env)
940 
941  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
942  TYPE(qs_environment_type), POINTER :: qs_env
943 
944  TYPE(dft_control_type), POINTER :: dft_control
945 
946  !PB only available without exact exchange
947  IF (xas_tdp_control%is_periodic .AND. xas_tdp_control%do_hfx &
948  .AND. xas_tdp_control%x_potential%potential_type == do_potential_coulomb) THEN
949  cpabort("XAS TDP with Coulomb operator for exact exchange only supports non-periodic BCs")
950  END IF
951 
952  !open-shell/closed-shell tests
953  IF (xas_tdp_control%do_roks .OR. xas_tdp_control%do_uks) THEN
954 
955  IF (.NOT. (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_spin_flip)) THEN
956  cpabort("Need spin-conserving and/or spin-flip excitations for open-shell systems")
957  END IF
958 
959  IF (xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet) THEN
960  cpabort("Singlet/triplet excitations only for restricted closed-shell systems")
961  END IF
962 
963  IF (xas_tdp_control%do_soc .AND. .NOT. &
964  (xas_tdp_control%do_spin_flip .AND. xas_tdp_control%do_spin_cons)) THEN
965 
966  cpabort("Both spin-conserving and spin-flip excitations are required for SOC")
967  END IF
968  ELSE
969 
970  IF (.NOT. (xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet)) THEN
971  cpabort("Need singlet and/or triplet excitations for closed-shell systems")
972  END IF
973 
974  IF (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_spin_flip) THEN
975  cpabort("Spin-conserving/spin-flip excitations only for open-shell systems")
976  END IF
977 
978  IF (xas_tdp_control%do_soc .AND. .NOT. &
979  (xas_tdp_control%do_singlet .AND. xas_tdp_control%do_triplet)) THEN
980 
981  cpabort("Both singlet and triplet excitations are needed for SOC")
982  END IF
983  END IF
984 
985  !Warn against using E_RANGE with SOC
986  IF (xas_tdp_control%do_soc .AND. xas_tdp_control%e_range > 0.0_dp) THEN
987  cpwarn("Using E_RANGE and SOC together may lead to crashes, use N_EXCITED for safety.")
988  END IF
989 
990  !TDA, full-TDDFT and diagonalization
991  IF (.NOT. xas_tdp_control%tamm_dancoff) THEN
992 
993  IF (xas_tdp_control%do_spin_flip) THEN
994  cpabort("Spin-flip kernel only implemented for Tamm-Dancoff approximation")
995  END IF
996 
997  IF (xas_tdp_control%do_ot) THEN
998  cpabort("OT diagonalization only available within the Tamm-Dancoff approximation")
999  END IF
1000  END IF
1001 
1002  !GW2X, need hfx kernel and LOCALIZE
1003  IF (xas_tdp_control%do_gw2x) THEN
1004  IF (.NOT. xas_tdp_control%do_hfx) THEN
1005  cpabort("GW2x requires the definition of the EXACT_EXCHANGE kernel")
1006  END IF
1007  IF (.NOT. xas_tdp_control%do_loc) THEN
1008  cpabort("GW2X requires the LOCALIZE keyword in DONOR_STATES")
1009  END IF
1010  END IF
1011 
1012  !Only allow ADMM schemes that correct for eigenvalues
1013  CALL get_qs_env(qs_env, dft_control=dft_control)
1014  IF (dft_control%do_admm) THEN
1015  IF ((.NOT. qs_env%admm_env%purification_method == do_admm_purify_none) .AND. &
1016  (.NOT. qs_env%admm_env%purification_method == do_admm_purify_cauchy_subspace) .AND. &
1017  (.NOT. qs_env%admm_env%purification_method == do_admm_purify_mo_diag)) THEN
1018 
1019  cpabort("XAS_TDP only compatible with ADMM purification NONE, CAUCHY_SUBSPACE and MO_DIAG")
1020 
1021  END IF
1022  END IF
1023 
1024  END SUBROUTINE safety_check
1025 
1026 ! **************************************************************************************************
1027 !> \brief Prints some basic info about the chosen parameters
1028 !> \param ou the output unis
1029 !> \param xas_tdp_control ...
1030 !> \param qs_env ...
1031 ! **************************************************************************************************
1032  SUBROUTINE print_info(ou, xas_tdp_control, qs_env)
1033 
1034  INTEGER, INTENT(IN) :: ou
1035  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1036  TYPE(qs_environment_type), POINTER :: qs_env
1037 
1038  INTEGER :: i
1039  REAL(dp) :: occ
1040  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1041  TYPE(dft_control_type), POINTER :: dft_control
1042  TYPE(section_vals_type), POINTER :: input, kernel_section
1043 
1044  NULLIFY (input, kernel_section, dft_control, matrix_s)
1045 
1046  CALL get_qs_env(qs_env, input=input, dft_control=dft_control, matrix_s=matrix_s)
1047 
1048  !Overlap matrix sparsity
1049  occ = dbcsr_get_occupation(matrix_s(1)%matrix)
1050 
1051  IF (ou .LE. 0) RETURN
1052 
1053  !Reference calculation
1054  IF (xas_tdp_control%do_uks) THEN
1055  WRITE (unit=ou, fmt="(/,T3,A)") &
1056  "XAS_TDP| Reference calculation: Unrestricted Kohn-Sham"
1057  ELSE IF (xas_tdp_control%do_roks) THEN
1058  WRITE (unit=ou, fmt="(/,T3,A)") &
1059  "XAS_TDP| Reference calculation: Restricted Open-Shell Kohn-Sham"
1060  ELSE
1061  WRITE (unit=ou, fmt="(/,T3,A)") &
1062  "XAS_TDP| Reference calculation: Restricted Closed-Shell Kohn-Sham"
1063  END IF
1064 
1065  !TDA
1066  IF (xas_tdp_control%tamm_dancoff) THEN
1067  WRITE (unit=ou, fmt="(T3,A)") &
1068  "XAS_TDP| Tamm-Dancoff Approximation (TDA): On"
1069  ELSE
1070  WRITE (unit=ou, fmt="(T3,A)") &
1071  "XAS_TDP| Tamm-Dancoff Approximation (TDA): Off"
1072  END IF
1073 
1074  !Dipole form
1075  IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
1076  WRITE (unit=ou, fmt="(T3,A)") &
1077  "XAS_TDP| Transition Dipole Representation: VELOCITY"
1078  ELSE
1079  WRITE (unit=ou, fmt="(T3,A)") &
1080  "XAS_TDP| Transition Dipole Representation: LENGTH"
1081  END IF
1082 
1083  !Quadrupole
1084  IF (xas_tdp_control%do_quad) THEN
1085  WRITE (unit=ou, fmt="(T3,A)") &
1086  "XAS_TDP| Transition Quadrupole: On"
1087  END IF
1088 
1089  !EPS_PGF
1090  IF (xas_tdp_control%eps_pgf > 0.0_dp) THEN
1091  WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1092  "XAS_TDP| EPS_PGF_XAS: ", xas_tdp_control%eps_pgf
1093  ELSE
1094  WRITE (unit=ou, fmt="(T3,A,ES7.1,A)") &
1095  "XAS_TDP| EPS_PGF_XAS: ", dft_control%qs_control%eps_pgf_orb, " (= EPS_PGF_ORB)"
1096  END IF
1097 
1098  !EPS_FILTER
1099  WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1100  "XAS_TDP| EPS_FILTER: ", xas_tdp_control%eps_filter
1101 
1102  !Grid info
1103  IF (xas_tdp_control%do_xc) THEN
1104  WRITE (unit=ou, fmt="(T3,A)") &
1105  "XAS_TDP| Radial Grid(s) Info: Kind, na, nr"
1106  DO i = 1, SIZE(xas_tdp_control%grid_info, 1)
1107  WRITE (unit=ou, fmt="(T3,A,A6,A,A,A,A)") &
1108  " ", trim(xas_tdp_control%grid_info(i, 1)), ", ", &
1109  trim(xas_tdp_control%grid_info(i, 2)), ", ", trim(xas_tdp_control%grid_info(i, 3))
1110  END DO
1111  END IF
1112 
1113  !No kernel
1114  IF (.NOT. xas_tdp_control%do_coulomb) THEN
1115  WRITE (unit=ou, fmt="(/,T3,A)") &
1116  "XAS_TDP| No kernel (standard DFT)"
1117  END IF
1118 
1119  !XC kernel
1120  IF (xas_tdp_control%do_xc) THEN
1121 
1122  WRITE (unit=ou, fmt="(/,T3,A,F5.2,A)") &
1123  "XAS_TDP| RI Region's Radius: ", xas_tdp_control%ri_radius*angstrom, " Ang"
1124 
1125  WRITE (unit=ou, fmt="(T3,A,/)") &
1126  "XAS_TDP| XC Kernel Functional(s) used for the kernel:"
1127 
1128  kernel_section => section_vals_get_subs_vals(input, "DFT%XAS_TDP%KERNEL")
1129  CALL xc_write(ou, kernel_section, lsd=.true.)
1130  END IF
1131 
1132  !HFX kernel
1133  IF (xas_tdp_control%do_hfx) THEN
1134  WRITE (unit=ou, fmt="(/,T3,A,/,/,T3,A,F5.3)") &
1135  "XAS_TDP| Exact Exchange Kernel: Yes ", &
1136  "EXACT_EXCHANGE| Scale: ", xas_tdp_control%sx
1137  IF (xas_tdp_control%x_potential%potential_type == do_potential_coulomb) THEN
1138  WRITE (unit=ou, fmt="(T3,A)") &
1139  "EXACT_EXCHANGE| Potential : Coulomb"
1140  ELSE IF (xas_tdp_control%x_potential%potential_type == do_potential_truncated) THEN
1141  WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,A)") &
1142  "EXACT_EXCHANGE| Potential: Truncated Coulomb", &
1143  "EXACT_EXCHANGE| Range: ", xas_tdp_control%x_potential%cutoff_radius*angstrom, ", (Ang)", &
1144  "EXACT_EXCHANGE| T_C_G_DATA: ", trim(xas_tdp_control%x_potential%filename)
1145  ELSE IF (xas_tdp_control%x_potential%potential_type == do_potential_short) THEN
1146  WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,F5.2,A,/,T3,A,ES7.1)") &
1147  "EXACT_EXCHANGE| Potential: Short Range", &
1148  "EXACT_EXCHANGE| Omega: ", xas_tdp_control%x_potential%omega, ", (1/a0)", &
1149  "EXACT_EXCHANGE| Effective Range: ", xas_tdp_control%x_potential%cutoff_radius*angstrom, ", (Ang)", &
1150  "EXACT_EXCHANGE| EPS_RANGE: ", xas_tdp_control%eps_range
1151  END IF
1152  IF (xas_tdp_control%eps_screen > 1.0e-16) THEN
1153  WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1154  "EXACT_EXCHANGE| EPS_SCREENING: ", xas_tdp_control%eps_screen
1155  END IF
1156 
1157  !RI metric
1158  IF (xas_tdp_control%do_ri_metric) THEN
1159 
1160  WRITE (unit=ou, fmt="(/,T3,A)") &
1161  "EXACT_EXCHANGE| Using a RI metric"
1162  IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_id) THEN
1163  WRITE (unit=ou, fmt="(T3,A)") &
1164  "EXACT_EXCHANGE RI_METRIC| Potential : Overlap"
1165  ELSE IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_truncated) THEN
1166  WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,A)") &
1167  "EXACT_EXCHANGE RI_METRIC| Potential: Truncated Coulomb", &
1168  "EXACT_EXCHANGE RI_METRIC| Range: ", xas_tdp_control%ri_m_potential%cutoff_radius &
1169  *angstrom, ", (Ang)", &
1170  "EXACT_EXCHANGE RI_METRIC| T_C_G_DATA: ", trim(xas_tdp_control%ri_m_potential%filename)
1171  ELSE IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_short) THEN
1172  WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,F5.2,A,/,T3,A,ES7.1)") &
1173  "EXACT_EXCHANGE RI_METRIC| Potential: Short Range", &
1174  "EXACT_EXCHANGE RI_METRIC| Omega: ", xas_tdp_control%ri_m_potential%omega, ", (1/a0)", &
1175  "EXACT_EXCHANGE RI_METRIC| Effective Range: ", &
1176  xas_tdp_control%ri_m_potential%cutoff_radius*angstrom, ", (Ang)", &
1177  "EXACT_EXCHANGE RI_METRIC| EPS_RANGE: ", xas_tdp_control%eps_range
1178  END IF
1179  END IF
1180  ELSE
1181  WRITE (unit=ou, fmt="(/,T3,A,/)") &
1182  "XAS_TDP| Exact Exchange Kernel: No "
1183  END IF
1184 
1185  !overlap mtrix occupation
1186  WRITE (unit=ou, fmt="(/,T3,A,F5.2)") &
1187  "XAS_TDP| Overlap matrix occupation: ", occ
1188 
1189  !GW2X parameter
1190  IF (xas_tdp_control%do_gw2x) THEN
1191  WRITE (unit=ou, fmt="(T3,A,/)") &
1192  "XAS_TDP| GW2X correction enabled"
1193 
1194  IF (xas_tdp_control%xps_only) THEN
1195  WRITE (unit=ou, fmt="(T3,A)") &
1196  "GW2X| Only computing ionizations potentials for XPS"
1197  END IF
1198 
1199  IF (xas_tdp_control%pseudo_canonical) THEN
1200  WRITE (unit=ou, fmt="(T3,A)") &
1201  "GW2X| Using the pseudo-canonical scheme"
1202  ELSE
1203  WRITE (unit=ou, fmt="(T3,A)") &
1204  "GW2X| Using the GW2X* scheme"
1205  END IF
1206 
1207  WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1208  "GW2X| EPS_GW2X: ", xas_tdp_control%gw2x_eps
1209 
1210  WRITE (unit=ou, fmt="(T3,A,I5)") &
1211  "GW2X| contraction batch size: ", xas_tdp_control%batch_size
1212 
1213  IF ((int(xas_tdp_control%c_os) .NE. 1) .OR. (int(xas_tdp_control%c_ss) .NE. 1)) THEN
1214  WRITE (unit=ou, fmt="(T3,A,F7.4,/,T3,A,F7.4)") &
1215  "GW2X| Same-spin scaling factor: ", xas_tdp_control%c_ss, &
1216  "GW2X| Opposite-spin scaling factor: ", xas_tdp_control%c_os
1217  END IF
1218 
1219  END IF
1220 
1221  END SUBROUTINE print_info
1222 
1223 ! **************************************************************************************************
1224 !> \brief Assosciate (possibly localized) lowest energy MOs to each excited atoms. The procedure
1225 !> looks for MOs "centered" on the excited atoms by comparing distances. It
1226 !> then fills the mos_of_ex_atoms arrays of the xas_tdp_env. Only the xas_tdp_control%n_search
1227 !> lowest energy MOs are considered. Largely inspired by MI's implementation of XAS
1228 !> It is assumed that the Berry phase is used to compute centers.
1229 !> \param xas_tdp_env ...
1230 !> \param xas_tdp_control ...
1231 !> \param qs_env ...
1232 !> \note Whether localization took place or not, the procedure is the same as centers are stored in
1233 !> xas_tdp_env%qs_loc_env%localized_wfn_control%centers_set
1234 !> Assumes that find_mo_centers has been run previously
1235 ! **************************************************************************************************
1236  SUBROUTINE assign_mos_to_ex_atoms(xas_tdp_env, xas_tdp_control, qs_env)
1237 
1238  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1239  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1240  TYPE(qs_environment_type), POINTER :: qs_env
1241 
1242  INTEGER :: at_index, iat, iat_memo, imo, ispin, &
1243  n_atoms, n_search, nex_atoms, nspins
1244  INTEGER, DIMENSION(3) :: perd_init
1245  INTEGER, DIMENSION(:, :, :), POINTER :: mos_of_ex_atoms
1246  REAL(dp) :: dist, dist_min
1247  REAL(dp), DIMENSION(3) :: at_pos, r_ac, wfn_center
1248  TYPE(cell_type), POINTER :: cell
1249  TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1250  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1251 
1252  NULLIFY (localized_wfn_control, mos_of_ex_atoms, cell, particle_set)
1253 
1254 ! Initialization. mos_of_ex_atoms filled with -1, meaning no assigned state
1255  mos_of_ex_atoms => xas_tdp_env%mos_of_ex_atoms
1256  mos_of_ex_atoms(:, :, :) = -1
1257  n_search = xas_tdp_control%n_search
1258  nex_atoms = xas_tdp_env%nex_atoms
1259  localized_wfn_control => xas_tdp_env%qs_loc_env%localized_wfn_control
1260  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
1261  n_atoms = SIZE(particle_set)
1262  nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1263 
1264 ! Temporarly impose periodic BCs because of Berry's phase operator used for localization
1265  perd_init = cell%perd
1266  cell%perd = 1
1267 
1268 ! Loop over n_search lowest energy MOs and all atoms, for each spin
1269  DO ispin = 1, nspins
1270  DO imo = 1, n_search
1271 ! retrieve MO wave function center coordinates.
1272  wfn_center(1:3) = localized_wfn_control%centers_set(ispin)%array(1:3, imo)
1273  iat_memo = 0
1274 
1275 ! a large enough value to avoid bad surprises
1276  dist_min = 10000.0_dp
1277  DO iat = 1, n_atoms
1278  at_pos = particle_set(iat)%r
1279  r_ac = pbc(at_pos, wfn_center, cell)
1280  dist = norm2(r_ac)
1281 
1282 ! keep memory of which atom is the closest to the wave function center
1283  IF (dist < dist_min) THEN
1284  iat_memo = iat
1285  dist_min = dist
1286  END IF
1287  END DO
1288 
1289 ! Verify that the closest atom is actually excited and assign the MO if so
1290  IF (any(xas_tdp_env%ex_atom_indices == iat_memo)) THEN
1291  at_index = locate(xas_tdp_env%ex_atom_indices, iat_memo)
1292  mos_of_ex_atoms(imo, at_index, ispin) = 1
1293  END IF
1294  END DO !imo
1295  END DO !ispin
1296 
1297 ! Go back to initial BCs
1298  cell%perd = perd_init
1299 
1300  END SUBROUTINE assign_mos_to_ex_atoms
1301 
1302 ! **************************************************************************************************
1303 !> \brief Re-initialize the qs_loc_env to the current MOs.
1304 !> \param qs_loc_env the env to re-initialize
1305 !> \param n_loc_states the number of states to include
1306 !> \param do_uks in cas of spin unrestricted calculation, initialize for both spins
1307 !> \param qs_env ...
1308 !> \note Useful when one needs to make use of qs_loc features and it is either with canonical MOs
1309 !> or the localized MOs have been modified. do_localize is overwritten.
1310 !> Same loc range for both spins
1311 ! **************************************************************************************************
1312  SUBROUTINE reinit_qs_loc_env(qs_loc_env, n_loc_states, do_uks, qs_env)
1313 
1314  TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1315  INTEGER, INTENT(IN) :: n_loc_states
1316  LOGICAL, INTENT(IN) :: do_uks
1317  TYPE(qs_environment_type), POINTER :: qs_env
1318 
1319  INTEGER :: i, nspins
1320  TYPE(localized_wfn_control_type), POINTER :: loc_wfn_control
1321 
1322 ! First, release the old env
1323  CALL qs_loc_env_release(qs_loc_env)
1324 
1325 ! Re-create it
1326  CALL qs_loc_env_create(qs_loc_env)
1327  CALL localized_wfn_control_create(qs_loc_env%localized_wfn_control)
1328  loc_wfn_control => qs_loc_env%localized_wfn_control
1329 
1330 ! Initialize it
1331  loc_wfn_control%localization_method = do_loc_none
1332  loc_wfn_control%operator_type = op_loc_berry
1333  loc_wfn_control%nloc_states(:) = n_loc_states
1334  loc_wfn_control%eps_occ = 0.0_dp
1335  loc_wfn_control%lu_bound_states(1, :) = 1
1336  loc_wfn_control%lu_bound_states(2, :) = n_loc_states
1337  loc_wfn_control%set_of_states = state_loc_list
1338  loc_wfn_control%do_homo = .true.
1339  ALLOCATE (loc_wfn_control%loc_states(n_loc_states, 2))
1340  DO i = 1, n_loc_states
1341  loc_wfn_control%loc_states(i, :) = i
1342  END DO
1343 
1344  nspins = 1; IF (do_uks) nspins = 2
1345  CALL set_loc_centers(loc_wfn_control, loc_wfn_control%nloc_states, nspins=nspins)
1346  ! need to set do_localize=.TRUE. because otherwise no routine works
1347  IF (do_uks) THEN
1348  CALL qs_loc_env_init(qs_loc_env, loc_wfn_control, qs_env, do_localize=.true.)
1349  ELSE
1350  CALL qs_loc_env_init(qs_loc_env, loc_wfn_control, qs_env, myspin=1, do_localize=.true.)
1351  END IF
1352 
1353  END SUBROUTINE reinit_qs_loc_env
1354 
1355 ! *************************************************************************************************
1356 !> \brief Diagonalize the subset of previously localized MOs that are associated to each excited
1357 !> atoms. Updates the MO coeffs accordingly.
1358 !> \param xas_tdp_env ...
1359 !> \param xas_tdp_control ...
1360 !> \param qs_env ...
1361 !> \note Needed because after localization, the MOs loose their identity (1s, 2s , 2p, etc)
1362 ! **************************************************************************************************
1363  SUBROUTINE diagonalize_assigned_mo_subset(xas_tdp_env, xas_tdp_control, qs_env)
1364 
1365  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1366  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1367  TYPE(qs_environment_type), POINTER :: qs_env
1368 
1369  INTEGER :: i, iat, ilmo, ispin, nao, nlmo, nspins
1370  REAL(dp), ALLOCATABLE, DIMENSION(:) :: evals
1371  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1372  TYPE(cp_fm_struct_type), POINTER :: ks_struct, lmo_struct
1373  TYPE(cp_fm_type) :: evecs, ks_fm, lmo_fm, work
1374  TYPE(cp_fm_type), POINTER :: mo_coeff
1375  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1376  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1377  TYPE(mp_para_env_type), POINTER :: para_env
1378 
1379  NULLIFY (mos, mo_coeff, matrix_ks, para_env, blacs_env, lmo_struct, ks_struct)
1380 
1381  ! Get what we need from qs_env
1382  CALL get_qs_env(qs_env, mos=mos, matrix_ks=matrix_ks, para_env=para_env, blacs_env=blacs_env)
1383 
1384  nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1385 
1386  ! Loop over the excited atoms and spin
1387  DO ispin = 1, nspins
1388  DO iat = 1, xas_tdp_env%nex_atoms
1389 
1390  ! get the MOs
1391  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao)
1392 
1393  ! count how many MOs are associated to this atom and create a fm/struct
1394  nlmo = count(xas_tdp_env%mos_of_ex_atoms(:, iat, ispin) == 1)
1395  CALL cp_fm_struct_create(lmo_struct, nrow_global=nao, ncol_global=nlmo, &
1396  para_env=para_env, context=blacs_env)
1397  CALL cp_fm_create(lmo_fm, lmo_struct)
1398  CALL cp_fm_create(work, lmo_struct)
1399 
1400  CALL cp_fm_struct_create(ks_struct, nrow_global=nlmo, ncol_global=nlmo, &
1401  para_env=para_env, context=blacs_env)
1402  CALL cp_fm_create(ks_fm, ks_struct)
1403  CALL cp_fm_create(evecs, ks_struct)
1404 
1405  ! Loop over the localized MOs associated to this atom
1406  i = 0
1407  DO ilmo = 1, xas_tdp_control%n_search
1408  IF (xas_tdp_env%mos_of_ex_atoms(ilmo, iat, ispin) == -1) cycle
1409 
1410  i = i + 1
1411  ! put the coeff in our atom-restricted lmo_fm
1412  CALL cp_fm_to_fm_submat(mo_coeff, lmo_fm, nrow=nao, ncol=1, s_firstrow=1, &
1413  s_firstcol=ilmo, t_firstrow=1, t_firstcol=i)
1414 
1415  END DO !ilmo
1416 
1417  ! Computing the KS matrix in the subset of MOs
1418  CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, lmo_fm, work, ncol=nlmo)
1419  CALL parallel_gemm('T', 'N', nlmo, nlmo, nao, 1.0_dp, lmo_fm, work, 0.0_dp, ks_fm)
1420 
1421  ! Diagonalizing the KS matrix in the subset of MOs
1422  ALLOCATE (evals(nlmo))
1423  CALL cp_fm_syevd(ks_fm, evecs, evals)
1424  DEALLOCATE (evals)
1425 
1426  ! Express the MOs in the basis that diagonalizes KS
1427  CALL parallel_gemm('N', 'N', nao, nlmo, nlmo, 1.0_dp, lmo_fm, evecs, 0.0_dp, work)
1428 
1429  ! Replacing the new MOs back in the MO coeffs
1430  i = 0
1431  DO ilmo = 1, xas_tdp_control%n_search
1432  IF (xas_tdp_env%mos_of_ex_atoms(ilmo, iat, ispin) == -1) cycle
1433 
1434  i = i + 1
1435  CALL cp_fm_to_fm_submat(work, mo_coeff, nrow=nao, ncol=1, s_firstrow=1, &
1436  s_firstcol=i, t_firstrow=1, t_firstcol=ilmo)
1437 
1438  END DO
1439 
1440  ! Excited atom clean-up
1441  CALL cp_fm_release(lmo_fm)
1442  CALL cp_fm_release(work)
1443  CALL cp_fm_struct_release(lmo_struct)
1444  CALL cp_fm_release(ks_fm)
1445  CALL cp_fm_release(evecs)
1446  CALL cp_fm_struct_release(ks_struct)
1447  END DO !iat
1448  END DO !ispin
1449 
1450  END SUBROUTINE diagonalize_assigned_mo_subset
1451 
1452 ! **************************************************************************************************
1453 !> \brief Assign core MO(s) to a given donor_state, taking the type (1S, 2S, etc) into account.
1454 !> The projection on a representative Slater-type orbital basis is used as a indicator.
1455 !> It is assumed that MOs are already assigned to excited atoms based on their center
1456 !> \param donor_state the donor_state to which a MO must be assigned
1457 !> \param xas_tdp_env ...
1458 !> \param xas_tdp_control ...
1459 !> \param qs_env ...
1460 ! **************************************************************************************************
1461  SUBROUTINE assign_mos_to_donor_state(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
1462 
1463  TYPE(donor_state_type), POINTER :: donor_state
1464  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1465  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1466  TYPE(qs_environment_type), POINTER :: qs_env
1467 
1468  INTEGER :: at_index, i, iat, imo, ispin, l, my_mo, &
1469  n_search, n_states, nao, ndo_so, nj, &
1470  nsgf_kind, nsgf_sto, nspins, &
1471  output_unit, zval
1472  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: my_mos
1473  INTEGER, DIMENSION(2) :: next_best_overlap_ind
1474  INTEGER, DIMENSION(4, 7) :: ne
1475  INTEGER, DIMENSION(:), POINTER :: first_sgf, lq, nq
1476  INTEGER, DIMENSION(:, :, :), POINTER :: mos_of_ex_atoms
1477  LOGICAL :: unique
1478  REAL(dp) :: zeff
1479  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, overlap, sto_overlap
1480  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: max_overlap
1481  REAL(dp), DIMENSION(2) :: next_best_overlap
1482  REAL(dp), DIMENSION(:), POINTER :: mo_evals, zeta
1483  REAL(dp), DIMENSION(:, :), POINTER :: overlap_matrix, tmp_coeff
1484  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1485  TYPE(cp_fm_struct_type), POINTER :: eval_mat_struct, gs_struct
1486  TYPE(cp_fm_type) :: eval_mat, work_mat
1487  TYPE(cp_fm_type), POINTER :: gs_coeffs, mo_coeff
1488  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1489  TYPE(gto_basis_set_type), POINTER :: kind_basis_set, sto_to_gto_basis_set
1490  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1491  TYPE(mp_para_env_type), POINTER :: para_env
1492  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1493  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1494  TYPE(sto_basis_set_type), POINTER :: sto_basis_set
1495 
1496  NULLIFY (sto_basis_set, sto_to_gto_basis_set, qs_kind_set, kind_basis_set, lq, nq, zeta)
1497  NULLIFY (overlap_matrix, mos, mo_coeff, mos_of_ex_atoms, tmp_coeff, first_sgf, particle_set)
1498  NULLIFY (mo_evals, matrix_ks, para_env, blacs_env)
1499  NULLIFY (eval_mat_struct, gs_struct, gs_coeffs)
1500 
1501  output_unit = cp_logger_get_default_io_unit()
1502 
1503  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, mos=mos, particle_set=particle_set, &
1504  matrix_ks=matrix_ks, para_env=para_env, blacs_env=blacs_env)
1505 
1506  nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1507 
1508 ! Construction of a STO that fits the type of orbital we look for
1509  ALLOCATE (zeta(1))
1510  ALLOCATE (lq(1))
1511  ALLOCATE (nq(1))
1512 ! Retrieving quantum numbers
1513  IF (donor_state%state_type == xas_1s_type) THEN
1514  nq(1) = 1
1515  lq(1) = 0
1516  n_states = 1
1517  ELSE IF (donor_state%state_type == xas_2s_type) THEN
1518  nq(1) = 2
1519  lq(1) = 0
1520  n_states = 1
1521  ELSE IF (donor_state%state_type == xas_2p_type) THEN
1522  nq(1) = 2
1523  lq(1) = 1
1524  n_states = 3
1525  ELSE
1526  cpabort("Procedure for required type not implemented")
1527  END IF
1528  ALLOCATE (my_mos(n_states, nspins))
1529  ALLOCATE (max_overlap(n_states, nspins))
1530 
1531 ! Getting the atomic number
1532  CALL get_qs_kind(qs_kind_set(donor_state%kind_index), zeff=zeff)
1533  zval = int(zeff)
1534 
1535 ! Electronic configuration (copied from MI's XAS)
1536  ne = 0
1537  DO l = 1, 4
1538  nj = 2*(l - 1) + 1
1539  DO i = l, 7
1540  ne(l, i) = ptable(zval)%e_conv(l - 1) - 2*nj*(i - l)
1541  ne(l, i) = max(ne(l, i), 0)
1542  ne(l, i) = min(ne(l, i), 2*nj)
1543  END DO
1544  END DO
1545 
1546 ! computing zeta with the Slater sum rules
1547  zeta(1) = srules(zval, ne, nq(1), lq(1))
1548 
1549 ! Allocating memory and initiate STO
1550  CALL allocate_sto_basis_set(sto_basis_set)
1551  CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, lq=lq, zet=zeta)
1552 
1553 ! Some clean-up
1554  DEALLOCATE (nq, lq, zeta)
1555 
1556 ! Expanding the STO into (normalized) GTOs for later calculations, use standard 3 gaussians
1557  CALL create_gto_from_sto_basis(sto_basis_set=sto_basis_set, &
1558  gto_basis_set=sto_to_gto_basis_set, &
1559  ngauss=3)
1560  sto_to_gto_basis_set%norm_type = 2
1561  CALL init_orb_basis_set(sto_to_gto_basis_set)
1562 
1563 ! Retrieving the atomic kind related GTO in which MOs are expanded
1564  CALL get_qs_kind(qs_kind_set(donor_state%kind_index), basis_set=kind_basis_set)
1565 
1566 ! Allocating and computing the overlap between the two basis (they share the same center)
1567  CALL get_gto_basis_set(gto_basis_set=kind_basis_set, nsgf=nsgf_kind)
1568  CALL get_gto_basis_set(gto_basis_set=sto_to_gto_basis_set, nsgf=nsgf_sto)
1569  ALLOCATE (overlap_matrix(nsgf_sto, nsgf_kind))
1570 
1571 ! Making use of MI's subroutine
1572  CALL calc_stogto_overlap(sto_to_gto_basis_set, kind_basis_set, overlap_matrix)
1573 
1574 ! Some clean-up
1575  CALL deallocate_sto_basis_set(sto_basis_set)
1576  CALL deallocate_gto_basis_set(sto_to_gto_basis_set)
1577 
1578 ! Looping over the potential donor states to compute overlap with STO basis
1579  mos_of_ex_atoms => xas_tdp_env%mos_of_ex_atoms
1580  n_search = xas_tdp_control%n_search
1581  at_index = donor_state%at_index
1582  iat = locate(xas_tdp_env%ex_atom_indices, at_index)
1583  ALLOCATE (first_sgf(SIZE(particle_set))) !probably do not need that
1584  CALL get_particle_set(particle_set=particle_set, qs_kind_set=qs_kind_set, first_sgf=first_sgf)
1585  ALLOCATE (tmp_coeff(nsgf_kind, 1))
1586  ALLOCATE (sto_overlap(nsgf_kind))
1587  ALLOCATE (overlap(n_search))
1588 
1589  next_best_overlap = 0.0_dp
1590  max_overlap = 0.0_dp
1591 
1592  DO ispin = 1, nspins
1593 
1594  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao)
1595  overlap = 0.0_dp
1596 
1597  my_mo = 0
1598  DO imo = 1, n_search
1599  IF (mos_of_ex_atoms(imo, iat, ispin) > 0) THEN
1600 
1601  sto_overlap = 0.0_dp
1602  tmp_coeff = 0.0_dp
1603 
1604 ! Getting the relevant coefficients for the candidate state
1605  CALL cp_fm_get_submatrix(fm=mo_coeff, target_m=tmp_coeff, start_row=first_sgf(at_index), &
1606  start_col=imo, n_rows=nsgf_kind, n_cols=1, transpose=.false.)
1607 
1608 ! Computing the product overlap_matrix*coeffs
1609  CALL dgemm('N', 'N', nsgf_sto, 1, nsgf_kind, 1.0_dp, overlap_matrix, nsgf_sto, &
1610  tmp_coeff, nsgf_kind, 0.0_dp, sto_overlap, nsgf_sto)
1611 
1612 ! Each element of column vector sto_overlap is the overlap of a basis element of the
1613 ! generated STO basis with the kind specific orbital basis. Take the sum of the absolute
1614 ! values so that rotation (of the px, py, pz for example) does not hinder our search
1615  overlap(imo) = sum(abs(sto_overlap))
1616 
1617  END IF
1618  END DO
1619 
1620 ! Finding the best overlap(s)
1621  DO i = 1, n_states
1622  my_mo = maxloc(overlap, 1)
1623  my_mos(i, ispin) = my_mo
1624  max_overlap(i, ispin) = maxval(overlap, 1)
1625  overlap(my_mo) = 0.0_dp
1626  END DO
1627 ! Getting the next best overlap (for validation purposes)
1628  next_best_overlap(ispin) = maxval(overlap, 1)
1629  next_best_overlap_ind(ispin) = maxloc(overlap, 1)
1630 
1631 ! Sort MO indices
1632  CALL sort_unique(my_mos(:, ispin), unique)
1633 
1634  END DO !ispin
1635 
1636 ! Some clean-up
1637  DEALLOCATE (overlap_matrix, tmp_coeff)
1638 
1639 ! Dealing with the result
1640  IF (all(my_mos > 0) .AND. all(my_mos <= n_search)) THEN
1641 ! Assigning the MO indices to the donor_state
1642  ALLOCATE (donor_state%mo_indices(n_states, nspins))
1643  donor_state%mo_indices = my_mos
1644  donor_state%ndo_mo = n_states
1645 
1646 ! Storing the MOs in the donor_state, as vectors column: first columns alpha spin, then beta
1647  CALL cp_fm_struct_create(gs_struct, nrow_global=nao, ncol_global=n_states*nspins, &
1648  para_env=para_env, context=blacs_env)
1649  ALLOCATE (donor_state%gs_coeffs)
1650  CALL cp_fm_create(donor_state%gs_coeffs, gs_struct)
1651 
1652  DO ispin = 1, nspins
1653  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1654  DO i = 1, n_states
1655  CALL cp_fm_to_fm_submat(msource=mo_coeff, mtarget=donor_state%gs_coeffs, nrow=nao, &
1656  ncol=1, s_firstrow=1, s_firstcol=my_mos(i, ispin), &
1657  t_firstrow=1, t_firstcol=(ispin - 1)*n_states + i)
1658  END DO
1659  END DO
1660  gs_coeffs => donor_state%gs_coeffs
1661 
1662  !Keep the subset of the coeffs centered on the excited atom as global array (used a lot)
1663  ALLOCATE (donor_state%contract_coeffs(nsgf_kind, n_states*nspins))
1664  CALL cp_fm_get_submatrix(gs_coeffs, donor_state%contract_coeffs, start_row=first_sgf(at_index), &
1665  start_col=1, n_rows=nsgf_kind, n_cols=n_states*nspins)
1666 
1667 ! Assigning corresponding energy eigenvalues and writing some info in standard input file
1668 
1669  !standard eigenvalues as gotten from the KS diagonalization in the ground state
1670  IF (.NOT. xas_tdp_control%do_loc .AND. .NOT. xas_tdp_control%do_roks) THEN
1671  IF (output_unit > 0) THEN
1672  WRITE (unit=output_unit, fmt="(T5,A,/,T5,A,/,T5,A)") &
1673  "The following canonical MO(s) have been associated with the donor state(s)", &
1674  "based on the overlap with the components of a minimal STO basis: ", &
1675  " Spin MO index overlap(sum)"
1676  END IF
1677 
1678  ALLOCATE (donor_state%energy_evals(n_states, nspins))
1679  donor_state%energy_evals = 0.0_dp
1680 
1681 ! Canonical MO, no change in eigenvalues, only diagonal elements
1682  DO ispin = 1, nspins
1683  CALL get_mo_set(mos(ispin), eigenvalues=mo_evals)
1684  DO i = 1, n_states
1685  donor_state%energy_evals(i, ispin) = mo_evals(my_mos(i, ispin))
1686 
1687  IF (output_unit > 0) THEN
1688  WRITE (unit=output_unit, fmt="(T46,I4,I11,F17.5)") &
1689  ispin, my_mos(i, ispin), max_overlap(i, ispin)
1690  END IF
1691  END DO
1692  END DO
1693 
1694  !either localization of MOs or ROKS, in both cases the MO eigenvalues from the KS
1695  !digonalization mat have changed
1696  ELSE
1697  IF (output_unit > 0) THEN
1698  WRITE (unit=output_unit, fmt="(T5,A,/,T5,A,/,T5,A)") &
1699  "The following localized MO(s) have been associated with the donor state(s)", &
1700  "based on the overlap with the components of a minimal STO basis: ", &
1701  " Spin MO index overlap(sum)"
1702  END IF
1703 
1704 ! Loop over the donor states and print
1705  DO ispin = 1, nspins
1706  DO i = 1, n_states
1707 
1708 ! Print info
1709  IF (output_unit > 0) THEN
1710  WRITE (unit=output_unit, fmt="(T46,I4,I11,F17.5)") &
1711  ispin, my_mos(i, ispin), max_overlap(i, ispin)
1712  END IF
1713  END DO
1714  END DO
1715 
1716 ! MO have been rotated or non-physical ROKS MO eigrenvalues:
1717 ! => need epsilon_ij = <psi_i|F|psi_j> = sum_{pq} c_{qi}c_{pj} F_{pq}
1718 ! Note: only have digonal elements by construction
1719  ndo_so = nspins*n_states
1720  CALL cp_fm_create(work_mat, gs_struct)
1721  CALL cp_fm_struct_create(eval_mat_struct, nrow_global=ndo_so, ncol_global=ndo_so, &
1722  para_env=para_env, context=blacs_env)
1723  CALL cp_fm_create(eval_mat, eval_mat_struct)
1724  ALLOCATE (diag(ndo_so))
1725 
1726  IF (.NOT. xas_tdp_control%do_roks) THEN
1727 
1728  ALLOCATE (donor_state%energy_evals(n_states, nspins))
1729  donor_state%energy_evals = 0.0_dp
1730 
1731 ! Compute gs_coeff^T * matrix_ks * gs_coeff to get the epsilon_ij matrix
1732  DO ispin = 1, nspins
1733  CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, gs_coeffs, work_mat, ncol=ndo_so)
1734  CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, work_mat, 0.0_dp, eval_mat)
1735 
1736 ! Put the epsilon_ii into the donor_state. No off-diagonal element because of subset diag
1737  CALL cp_fm_get_diag(eval_mat, diag)
1738  donor_state%energy_evals(:, ispin) = diag((ispin - 1)*n_states + 1:ispin*n_states)
1739 
1740  END DO
1741 
1742  ELSE
1743  ! If ROKS, slightly different procedure => 2 KS matrices but one type of MOs
1744  ALLOCATE (donor_state%energy_evals(n_states, 2))
1745  donor_state%energy_evals = 0.0_dp
1746 
1747 ! Compute gs_coeff^T * matrix_ks * gs_coeff to get the epsilon_ij matrix
1748  DO ispin = 1, 2
1749  CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, gs_coeffs, work_mat, ncol=ndo_so)
1750  CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, work_mat, 0.0_dp, eval_mat)
1751 
1752  CALL cp_fm_get_diag(eval_mat, diag)
1753  donor_state%energy_evals(:, ispin) = diag(:)
1754 
1755  END DO
1756 
1757  DEALLOCATE (diag)
1758  END IF
1759 
1760 ! Clean-up
1761  CALL cp_fm_release(work_mat)
1762  CALL cp_fm_release(eval_mat)
1763  CALL cp_fm_struct_release(eval_mat_struct)
1764 
1765  END IF ! do_localize and/or ROKS
1766 
1767 ! Allocate and initialize GW2X corrected IPs as energy_evals
1768  ALLOCATE (donor_state%gw2x_evals(SIZE(donor_state%energy_evals, 1), SIZE(donor_state%energy_evals, 2)))
1769  donor_state%gw2x_evals(:, :) = donor_state%energy_evals(:, :)
1770 
1771 ! Clean-up
1772  CALL cp_fm_struct_release(gs_struct)
1773  DEALLOCATE (first_sgf)
1774 
1775  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T5,A)") " "
1776 
1777  DO ispin = 1, nspins
1778  IF (output_unit > 0) THEN
1779  WRITE (unit=output_unit, fmt="(T5,A,I1,A,F7.5,A,I4)") &
1780  "The next best overlap for spin ", ispin, " is ", next_best_overlap(ispin), &
1781  " for MO with index ", next_best_overlap_ind(ispin)
1782  END IF
1783  END DO
1784  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T5,A)") " "
1785 
1786  ELSE
1787  cpabort("A core donor state could not be assigned MO(s). Increasing NSEARCH might help.")
1788  END IF
1789 
1790  END SUBROUTINE assign_mos_to_donor_state
1791 
1792 ! **************************************************************************************************
1793 !> \brief Compute the centers and spreads of (core) MOs using the Berry phase operator
1794 !> \param xas_tdp_env ...
1795 !> \param xas_tdp_control ...
1796 !> \param qs_env ...
1797 !> \note xas_tdp_env%qs_loc_env is used and modified. OK since no localization done after this
1798 !> subroutine is used.
1799 ! **************************************************************************************************
1800  SUBROUTINE find_mo_centers(xas_tdp_env, xas_tdp_control, qs_env)
1801 
1802  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1803  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1804  TYPE(qs_environment_type), POINTER :: qs_env
1805 
1806  INTEGER :: dim_op, i, ispin, j, n_centers, nao, &
1807  nspins
1808  REAL(dp), DIMENSION(6) :: weights
1809  TYPE(cell_type), POINTER :: cell
1810  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1811  TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1812  TYPE(cp_fm_type) :: opvec
1813  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij_fm_set
1814  TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
1815  TYPE(cp_fm_type), POINTER :: vectors
1816  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
1817  TYPE(mp_para_env_type), POINTER :: para_env
1818  TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1819  TYPE(section_vals_type), POINTER :: print_loc_section, prog_run_info
1820 
1821  NULLIFY (qs_loc_env, cell, print_loc_section, op_sm_set, moloc_coeff, vectors)
1822  NULLIFY (tmp_fm_struct, para_env, blacs_env, prog_run_info)
1823 
1824 ! Initialization
1825  print_loc_section => xas_tdp_control%print_loc_subsection
1826  n_centers = xas_tdp_control%n_search
1827  CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env, cell=cell)
1828 
1829 ! Set print option to debug to keep clean output file
1830  prog_run_info => section_vals_get_subs_vals(print_loc_section, "PROGRAM_RUN_INFO")
1831  CALL section_vals_val_set(prog_run_info, keyword_name="_SECTION_PARAMETERS_", &
1832  i_val=debug_print_level)
1833 
1834 ! Re-initialize the qs_loc_env to get the current MOs. Use force_loc because needed for centers
1835  CALL reinit_qs_loc_env(xas_tdp_env%qs_loc_env, n_centers, xas_tdp_control%do_uks, qs_env)
1836  qs_loc_env => xas_tdp_env%qs_loc_env
1837 
1838 ! Get what we need from the qs_lovc_env
1839  CALL get_qs_loc_env(qs_loc_env=qs_loc_env, weights=weights, op_sm_set=op_sm_set, &
1840  moloc_coeff=moloc_coeff)
1841 
1842 ! Prepare for zij
1843  vectors => moloc_coeff(1)
1844  CALL cp_fm_get_info(vectors, nrow_global=nao)
1845  CALL cp_fm_create(opvec, vectors%matrix_struct)
1846 
1847  CALL cp_fm_struct_create(tmp_fm_struct, para_env=para_env, context=blacs_env, &
1848  ncol_global=n_centers, nrow_global=n_centers)
1849 
1850  IF (cell%orthorhombic) THEN
1851  dim_op = 3
1852  ELSE
1853  dim_op = 6
1854  END IF
1855  ALLOCATE (zij_fm_set(2, dim_op))
1856  DO i = 1, dim_op
1857  DO j = 1, 2
1858  CALL cp_fm_create(zij_fm_set(j, i), tmp_fm_struct)
1859  END DO
1860  END DO
1861 
1862  ! If spin-unrestricted, need to go spin by spin
1863  nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1864 
1865  DO ispin = 1, nspins
1866 ! zij computation, copied from qs_loc_methods:optimize_loc_berry
1867  vectors => moloc_coeff(ispin)
1868  DO i = 1, dim_op
1869  DO j = 1, 2
1870  CALL cp_fm_set_all(zij_fm_set(j, i), 0.0_dp)
1871  CALL cp_dbcsr_sm_fm_multiply(op_sm_set(j, i)%matrix, vectors, opvec, ncol=n_centers)
1872  CALL parallel_gemm("T", "N", n_centers, n_centers, nao, 1.0_dp, vectors, opvec, 0.0_dp, &
1873  zij_fm_set(j, i))
1874  END DO
1875  END DO
1876 
1877 ! Compute centers (and spread)
1878  CALL centers_spreads_berry(qs_loc_env=qs_loc_env, zij=zij_fm_set, nmoloc=n_centers, &
1879  cell=cell, weights=weights, ispin=ispin, &
1880  print_loc_section=print_loc_section, only_initial_out=.true.)
1881  END DO !ispins
1882 
1883 ! Clean-up
1884  CALL cp_fm_release(opvec)
1885  CALL cp_fm_struct_release(tmp_fm_struct)
1886  CALL cp_fm_release(zij_fm_set)
1887 
1888 ! Make sure we leave with the correct do_loc value
1889  qs_loc_env%do_localize = xas_tdp_control%do_loc
1890 
1891  END SUBROUTINE find_mo_centers
1892 
1893 ! **************************************************************************************************
1894 !> \brief Prints the MO to donor_state assocaition with overlap and Mulliken population analysis
1895 !> \param xas_tdp_env ...
1896 !> \param xas_tdp_control ...
1897 !> \param qs_env ...
1898 !> \note Called only in case of CHECK_ONLY run
1899 ! **************************************************************************************************
1900  SUBROUTINE print_checks(xas_tdp_env, xas_tdp_control, qs_env)
1901 
1902  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1903  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1904  TYPE(qs_environment_type), POINTER :: qs_env
1905 
1906  CHARACTER(LEN=default_string_length) :: kind_name
1907  INTEGER :: current_state_index, iat, iatom, ikind, &
1908  istate, output_unit, tmp_index
1909  INTEGER, DIMENSION(:), POINTER :: atoms_of_kind
1910  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1911  TYPE(donor_state_type), POINTER :: current_state
1912 
1913  NULLIFY (atomic_kind_set, atoms_of_kind, current_state)
1914 
1915  output_unit = cp_logger_get_default_io_unit()
1916 
1917  IF (output_unit > 0) THEN
1918  WRITE (output_unit, "(/,T3,A,/,T3,A,/,T3,A)") &
1919  "# Check the donor states for their quality. They need to have a well defined type ", &
1920  " (1s, 2s, etc) which is indicated by the overlap. They also need to be localized, ", &
1921  " for which the Mulliken population analysis is one indicator (must be close to 1.0)"
1922  END IF
1923 
1924 ! Loop over the donor states (as in the main xas_tdp loop)
1925  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1926  current_state_index = 1
1927 
1928  !loop over atomic kinds
1929  DO ikind = 1, SIZE(atomic_kind_set)
1930 
1931  CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), name=kind_name, &
1932  atom_list=atoms_of_kind)
1933 
1934  IF (.NOT. any(xas_tdp_env%ex_kind_indices == ikind)) cycle
1935 
1936  !loop over atoms of kind
1937  DO iat = 1, SIZE(atoms_of_kind)
1938  iatom = atoms_of_kind(iat)
1939 
1940  IF (.NOT. any(xas_tdp_env%ex_atom_indices == iatom)) cycle
1941  tmp_index = locate(xas_tdp_env%ex_atom_indices, iatom)
1942 
1943  !loop over states of excited atom
1944  DO istate = 1, SIZE(xas_tdp_env%state_types, 1)
1945 
1946  IF (xas_tdp_env%state_types(istate, tmp_index) == xas_not_excited) cycle
1947 
1948  current_state => xas_tdp_env%donor_states(current_state_index)
1949  CALL set_donor_state(current_state, at_index=iatom, &
1950  at_symbol=kind_name, kind_index=ikind, &
1951  state_type=xas_tdp_env%state_types(istate, tmp_index))
1952 
1953  IF (output_unit > 0) THEN
1954  WRITE (output_unit, "(/,T4,A,A2,A,I4,A,A,A)") &
1955  "-Donor state of type ", xas_tdp_env%state_type_char(current_state%state_type), &
1956  " for atom", current_state%at_index, " of kind ", trim(current_state%at_symbol), ":"
1957  END IF
1958 
1959  !Assign the MOs and perform Mulliken
1960  CALL assign_mos_to_donor_state(current_state, xas_tdp_env, xas_tdp_control, qs_env)
1961  CALL perform_mulliken_on_donor_state(current_state, qs_env)
1962 
1963  current_state_index = current_state_index + 1
1964  NULLIFY (current_state)
1965 
1966  END DO !istate
1967  END DO !iat
1968  END DO !ikind
1969 
1970  IF (output_unit > 0) THEN
1971  WRITE (output_unit, "(/,T5,A)") &
1972  "Use LOCALIZE and/or increase N_SEARCH for better results, if so required."
1973  END IF
1974 
1975  END SUBROUTINE print_checks
1976 
1977 ! **************************************************************************************************
1978 !> \brief Computes the required multipole moment in the length representation for a given atom
1979 !> \param iatom index of the given atom
1980 !> \param xas_tdp_env ...
1981 !> \param xas_tdp_control ...
1982 !> \param qs_env ...
1983 !> \note Assumes that wither dipole or quadrupole in length rep is required
1984 ! **************************************************************************************************
1985  SUBROUTINE compute_lenrep_multipole(iatom, xas_tdp_env, xas_tdp_control, qs_env)
1986 
1987  INTEGER, INTENT(IN) :: iatom
1988  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1989  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1990  TYPE(qs_environment_type), POINTER :: qs_env
1991 
1992  INTEGER :: i, order
1993  REAL(dp), DIMENSION(3) :: rc
1994  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: work
1995  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1996 
1997  NULLIFY (work, particle_set)
1998 
1999  CALL get_qs_env(qs_env, particle_set=particle_set)
2000  rc = particle_set(iatom)%r
2001 
2002  ALLOCATE (work(9))
2003  IF (xas_tdp_control%dipole_form == xas_dip_len) THEN
2004  DO i = 1, 3
2005  CALL dbcsr_set(xas_tdp_env%dipmat(i)%matrix, 0.0_dp)
2006  work(i)%matrix => xas_tdp_env%dipmat(i)%matrix
2007  END DO
2008  order = 1
2009  END IF
2010  IF (xas_tdp_control%do_quad) THEN
2011  DO i = 1, 6
2012  CALL dbcsr_set(xas_tdp_env%quadmat(i)%matrix, 0.0_dp)
2013  work(3 + i)%matrix => xas_tdp_env%quadmat(i)%matrix
2014  END DO
2015  order = 2
2016  IF (xas_tdp_control%dipole_form == xas_dip_vel) order = -2
2017  END IF
2018 
2019  !enforce minimum image to avoid PBCs related issues, ok because localized densities
2020  CALL rrc_xyz_ao(work, qs_env, rc, order=order, minimum_image=.true.)
2021  DEALLOCATE (work)
2022 
2023  END SUBROUTINE compute_lenrep_multipole
2024 
2025 ! **************************************************************************************************
2026 !> \brief Computes the oscillator strength based on the dipole moment (velocity or length rep) for
2027 !> all available excitation energies and store the results in the donor_state. There is no
2028 !> triplet dipole in the spin-restricted ground state.
2029 !> \param donor_state the donor state which is excited
2030 !> \param xas_tdp_control ...
2031 !> \param xas_tdp_env ...
2032 !> \note The oscillator strength is a scalar: osc_str = 2/(3*omega)*(dipole_v)^2 in the velocity rep
2033 !> or : osc_str = 2/3*omega*(dipole_r)^2 in the length representation
2034 !> The formulae for the dipoles come from the trace of the dipole operator with the transition
2035 !> densities, i.e. what we get from solving the xas_tdp problem. Same procedure with or wo TDA
2036 ! **************************************************************************************************
2037  SUBROUTINE compute_dipole_fosc(donor_state, xas_tdp_control, xas_tdp_env)
2038 
2039  TYPE(donor_state_type), POINTER :: donor_state
2040  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2041  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2042 
2043  CHARACTER(len=*), PARAMETER :: routinen = 'compute_dipole_fosc'
2044 
2045  INTEGER :: handle, iosc, j, nao, ndo_mo, ndo_so, &
2046  ngs, nosc, nspins
2047  LOGICAL :: do_sc, do_sg
2048  REAL(dp) :: osc_xyz, pref
2049  REAL(dp), ALLOCATABLE, DIMENSION(:) :: tot_contr
2050  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dip_block
2051  REAL(dp), DIMENSION(:), POINTER :: lr_evals
2052  REAL(dp), DIMENSION(:, :), POINTER :: osc_str
2053  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2054  TYPE(cp_fm_struct_type), POINTER :: col_struct, mat_struct
2055  TYPE(cp_fm_type) :: col_work, mat_work
2056  TYPE(cp_fm_type), POINTER :: lr_coeffs
2057  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat
2058  TYPE(mp_para_env_type), POINTER :: para_env
2059 
2060  NULLIFY (dipmat, col_struct, mat_struct, para_env, blacs_env, lr_coeffs)
2061  NULLIFY (lr_evals, osc_str)
2062 
2063  CALL timeset(routinen, handle)
2064 
2065 ! Initialization
2066  do_sc = xas_tdp_control%do_spin_cons
2067  do_sg = xas_tdp_control%do_singlet
2068  IF (do_sc) THEN
2069  nspins = 2
2070  lr_evals => donor_state%sc_evals
2071  lr_coeffs => donor_state%sc_coeffs
2072  ELSE IF (do_sg) THEN
2073  nspins = 1
2074  lr_evals => donor_state%sg_evals
2075  lr_coeffs => donor_state%sg_coeffs
2076  ELSE
2077  cpabort("Dipole oscilaltor strength only for singlets and spin-conserving excitations.")
2078  END IF
2079  ndo_mo = donor_state%ndo_mo
2080  ndo_so = ndo_mo*nspins
2081  ngs = ndo_so; IF (xas_tdp_control%do_roks) ngs = ndo_mo !in ROKS, same gs coeffs
2082  nosc = SIZE(lr_evals)
2083  ALLOCATE (donor_state%osc_str(nosc, 4))
2084  osc_str => donor_state%osc_str
2085  osc_str = 0.0_dp
2086  dipmat => xas_tdp_env%dipmat
2087 
2088  ! do some work matrix initialization
2089  CALL cp_fm_get_info(donor_state%gs_coeffs, matrix_struct=col_struct, para_env=para_env, &
2090  context=blacs_env, nrow_global=nao)
2091  CALL cp_fm_struct_create(mat_struct, para_env=para_env, context=blacs_env, &
2092  nrow_global=ndo_so*nosc, ncol_global=ngs)
2093  CALL cp_fm_create(mat_work, mat_struct)
2094  CALL cp_fm_create(col_work, col_struct)
2095 
2096  ALLOCATE (tot_contr(ndo_mo), dip_block(ndo_so, ngs))
2097  pref = 2.0_dp; IF (do_sc) pref = 1.0_dp !because of singlet definition u = 1/sqrt(2)(c_a+c_b)
2098 
2099 ! Looping over cartesian coord
2100  DO j = 1, 3
2101 
2102  !Compute dip*gs_coeffs
2103  CALL cp_dbcsr_sm_fm_multiply(dipmat(j)%matrix, donor_state%gs_coeffs, col_work, ncol=ngs)
2104  !compute lr_coeffs*dip*gs_coeffs
2105  CALL parallel_gemm('T', 'N', ndo_so*nosc, ngs, nao, 1.0_dp, lr_coeffs, col_work, 0.0_dp, mat_work)
2106 
2107  !Loop over the excited states
2108  DO iosc = 1, nosc
2109 
2110  tot_contr = 0.0_dp
2111  CALL cp_fm_get_submatrix(fm=mat_work, target_m=dip_block, start_row=(iosc - 1)*ndo_so + 1, &
2112  start_col=1, n_rows=ndo_so, n_cols=ngs)
2113  IF (do_sg) THEN
2114  tot_contr(:) = get_diag(dip_block)
2115  ELSE IF (do_sc .AND. xas_tdp_control%do_uks) THEN
2116  tot_contr(:) = get_diag(dip_block(1:ndo_mo, 1:ndo_mo)) !alpha
2117  tot_contr(:) = tot_contr(:) + get_diag(dip_block(ndo_mo + 1:ndo_so, ndo_mo + 1:ndo_so)) !beta
2118  ELSE
2119  !roks
2120  tot_contr(:) = get_diag(dip_block(1:ndo_mo, :)) !alpha
2121  tot_contr(:) = tot_contr(:) + get_diag(dip_block(ndo_mo + 1:ndo_so, :)) !beta
2122  END IF
2123 
2124  osc_xyz = sum(tot_contr)**2
2125  osc_str(iosc, 4) = osc_str(iosc, 4) + osc_xyz
2126  osc_str(iosc, j) = osc_xyz
2127 
2128  END DO !iosc
2129  END DO !j
2130 
2131  !compute the prefactor
2132  DO j = 1, 4
2133  IF (xas_tdp_control%dipole_form == xas_dip_len) THEN
2134  osc_str(:, j) = pref*2.0_dp/3.0_dp*lr_evals(:)*osc_str(:, j)
2135  ELSE
2136  osc_str(:, j) = pref*2.0_dp/3.0_dp/lr_evals(:)*osc_str(:, j)
2137  END IF
2138  END DO
2139 
2140  !clean-up
2141  CALL cp_fm_release(mat_work)
2142  CALL cp_fm_release(col_work)
2143  CALL cp_fm_struct_release(mat_struct)
2144 
2145  CALL timestop(handle)
2146 
2147  END SUBROUTINE compute_dipole_fosc
2148 
2149 ! **************************************************************************************************
2150 !> \brief Computes the oscillator strength due to the electric quadrupole moment and store it in
2151 !> the donor_state (for singlet or spin-conserving)
2152 !> \param donor_state the donor state which is excited
2153 !> \param xas_tdp_control ...
2154 !> \param xas_tdp_env ...
2155 !> \note Formula: 1/20*a_fine^2*omega^3 * sum_ab (sum_i r_ia*r_ib - 1/3*ri^2*delta_ab)
2156 ! **************************************************************************************************
2157  SUBROUTINE compute_quadrupole_fosc(donor_state, xas_tdp_control, xas_tdp_env)
2158 
2159  TYPE(donor_state_type), POINTER :: donor_state
2160  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2161  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2162 
2163  CHARACTER(len=*), PARAMETER :: routinen = 'compute_quadrupole_fosc'
2164 
2165  INTEGER :: handle, iosc, j, nao, ndo_mo, ndo_so, &
2166  ngs, nosc, nspins
2167  LOGICAL :: do_sc, do_sg
2168  REAL(dp) :: pref
2169  REAL(dp), ALLOCATABLE, DIMENSION(:) :: tot_contr, trace
2170  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: quad_block
2171  REAL(dp), DIMENSION(:), POINTER :: lr_evals, osc_str
2172  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2173  TYPE(cp_fm_struct_type), POINTER :: col_struct, mat_struct
2174  TYPE(cp_fm_type) :: col_work, mat_work
2175  TYPE(cp_fm_type), POINTER :: lr_coeffs
2176  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: quadmat
2177  TYPE(mp_para_env_type), POINTER :: para_env
2178 
2179  NULLIFY (lr_evals, osc_str, lr_coeffs, col_struct, mat_struct, para_env)
2180  NULLIFY (blacs_env)
2181 
2182  CALL timeset(routinen, handle)
2183 
2184  ! Initialization
2185  do_sc = xas_tdp_control%do_spin_cons
2186  do_sg = xas_tdp_control%do_singlet
2187  IF (do_sc) THEN
2188  nspins = 2
2189  lr_evals => donor_state%sc_evals
2190  lr_coeffs => donor_state%sc_coeffs
2191  ELSE IF (do_sg) THEN
2192  nspins = 1
2193  lr_evals => donor_state%sg_evals
2194  lr_coeffs => donor_state%sg_coeffs
2195  ELSE
2196  cpabort("Quadrupole oscillator strengths only for singlet and spin-conserving excitations")
2197  END IF
2198  ndo_mo = donor_state%ndo_mo
2199  ndo_so = ndo_mo*nspins
2200  ngs = ndo_so; IF (xas_tdp_control%do_roks) ngs = ndo_mo !only alpha do_mo in ROKS
2201  nosc = SIZE(lr_evals)
2202  ALLOCATE (donor_state%quad_osc_str(nosc))
2203  osc_str => donor_state%quad_osc_str
2204  osc_str = 0.0_dp
2205  quadmat => xas_tdp_env%quadmat
2206 
2207  !work matrices init
2208  CALL cp_fm_get_info(donor_state%gs_coeffs, matrix_struct=col_struct, para_env=para_env, &
2209  context=blacs_env, nrow_global=nao)
2210  CALL cp_fm_struct_create(mat_struct, para_env=para_env, context=blacs_env, &
2211  nrow_global=ndo_so*nosc, ncol_global=ngs)
2212  CALL cp_fm_create(mat_work, mat_struct)
2213  CALL cp_fm_create(col_work, col_struct)
2214 
2215  ALLOCATE (quad_block(ndo_so, ngs), tot_contr(ndo_mo))
2216  pref = 2.0_dp; IF (do_sc) pref = 1.0_dp !because of singlet definition u = 1/sqrt(2)*...
2217  ALLOCATE (trace(nosc))
2218  trace = 0.0_dp
2219 
2220  !Loop over the cartesioan coord :x2, xy, xz, y2, yz, z2
2221  DO j = 1, 6
2222 
2223  !Compute quad*gs_coeffs
2224  CALL cp_dbcsr_sm_fm_multiply(quadmat(j)%matrix, donor_state%gs_coeffs, col_work, ncol=ngs)
2225  !compute lr_coeffs*quadmat*gs_coeffs
2226  CALL parallel_gemm('T', 'N', ndo_so*nosc, ngs, nao, 1.0_dp, lr_coeffs, col_work, 0.0_dp, mat_work)
2227 
2228  !Loop over the excited states
2229  DO iosc = 1, nosc
2230 
2231  tot_contr = 0.0_dp
2232  CALL cp_fm_get_submatrix(fm=mat_work, target_m=quad_block, start_row=(iosc - 1)*ndo_so + 1, &
2233  start_col=1, n_rows=ndo_so, n_cols=ngs)
2234 
2235  IF (do_sg) THEN
2236  tot_contr(:) = get_diag(quad_block)
2237  ELSE IF (do_sc .AND. xas_tdp_control%do_uks) THEN
2238  tot_contr(:) = get_diag(quad_block(1:ndo_mo, 1:ndo_mo)) !alpha
2239  tot_contr(:) = tot_contr(:) + get_diag(quad_block(ndo_mo + 1:ndo_so, ndo_mo + 1:ndo_so)) !beta
2240  ELSE
2241  !roks
2242  tot_contr(:) = get_diag(quad_block(1:ndo_mo, :)) !alpha
2243  tot_contr(:) = tot_contr(:) + get_diag(quad_block(ndo_mo + 1:ndo_so, :)) !beta
2244  END IF
2245 
2246  !if x2, y2, or z2 direction, need to update the trace (for later)
2247  IF (j == 1 .OR. j == 4 .OR. j == 6) THEN
2248  osc_str(iosc) = osc_str(iosc) + sum(tot_contr)**2
2249  trace(iosc) = trace(iosc) + sum(tot_contr)
2250 
2251  !if xy, xz or yz, need to count twice the contribution (for yx, zx and zy)
2252  ELSE
2253  osc_str(iosc) = osc_str(iosc) + 2.0_dp*sum(tot_contr)**2
2254  END IF
2255 
2256  END DO !iosc
2257  END DO !j
2258 
2259  !compute the prefactor, and remove 1/3*trace^2
2260  osc_str(:) = pref*1._dp/20._dp*a_fine**2*lr_evals(:)**3*(osc_str(:) - 1._dp/3._dp*trace(:)**2)
2261 
2262  !clean-up
2263  CALL cp_fm_release(mat_work)
2264  CALL cp_fm_release(col_work)
2265  CALL cp_fm_struct_release(mat_struct)
2266 
2267  CALL timestop(handle)
2268 
2269  END SUBROUTINE compute_quadrupole_fosc
2270 
2271 ! **************************************************************************************************
2272 !> \brief Writes the core MOs to excited atoms associations in the main output file
2273 !> \param xas_tdp_env ...
2274 !> \param xas_tdp_control ...
2275 !> \param qs_env ...
2276 !> \note Look at alpha spin MOs, as we are dealing with core states and alpha/beta MOs are the same
2277 ! **************************************************************************************************
2278  SUBROUTINE write_mos_to_ex_atoms_association(xas_tdp_env, xas_tdp_control, qs_env)
2279 
2280  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2281  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2282  TYPE(qs_environment_type), POINTER :: qs_env
2283 
2284  CHARACTER(LEN=default_string_length) :: kind_name
2285  INTEGER :: at_index, imo, ispin, nmo, nspins, &
2286  output_unit, tmp_index
2287  INTEGER, DIMENSION(3) :: perd_init
2288  INTEGER, DIMENSION(:), POINTER :: ex_atom_indices
2289  INTEGER, DIMENSION(:, :, :), POINTER :: mos_of_ex_atoms
2290  REAL(dp) :: dist, mo_spread
2291  REAL(dp), DIMENSION(3) :: at_pos, r_ac, wfn_center
2292  TYPE(cell_type), POINTER :: cell
2293  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2294 
2295  NULLIFY (cell, particle_set, mos_of_ex_atoms, ex_atom_indices)
2296 
2297  output_unit = cp_logger_get_default_io_unit()
2298 
2299  IF (output_unit > 0) THEN
2300  WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,/,T3,A)") &
2301  " Associated Associated Distance to MO spread (Ang^2)", &
2302  "Spin MO index atom index atom kind MO center (Ang) -w_i ln(|z_ij|^2)", &
2303  "---------------------------------------------------------------------------------"
2304  END IF
2305 
2306 ! Initialization
2307  nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
2308  mos_of_ex_atoms => xas_tdp_env%mos_of_ex_atoms
2309  ex_atom_indices => xas_tdp_env%ex_atom_indices
2310  nmo = xas_tdp_control%n_search
2311  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
2312 
2313 ! because the use of Berry's phase operator implies PBCs
2314  perd_init = cell%perd
2315  cell%perd = 1
2316 
2317 ! Retrieving all the info for each MO and spin
2318  DO imo = 1, nmo
2319  DO ispin = 1, nspins
2320 
2321 ! each Mo is associated to at most one atom (only 1 in array of -1)
2322  IF (any(mos_of_ex_atoms(imo, :, ispin) == 1)) THEN
2323  tmp_index = maxloc(mos_of_ex_atoms(imo, :, ispin), 1)
2324  at_index = ex_atom_indices(tmp_index)
2325  kind_name = particle_set(at_index)%atomic_kind%name
2326 
2327  at_pos = particle_set(at_index)%r
2328  wfn_center = xas_tdp_env%qs_loc_env%localized_wfn_control%centers_set(ispin)%array(1:3, imo)
2329  r_ac = pbc(at_pos, wfn_center, cell)
2330  dist = norm2(r_ac)
2331 ! convert distance from a.u. to Angstrom
2332  dist = dist*angstrom
2333 
2334  mo_spread = xas_tdp_env%qs_loc_env%localized_wfn_control%centers_set(ispin)%array(4, imo)
2335  mo_spread = mo_spread*angstrom*angstrom
2336 
2337  IF (output_unit > 0) THEN
2338  WRITE (unit=output_unit, fmt="(T3,I4,I10,I14,A14,ES19.3,ES20.3)") &
2339  ispin, imo, at_index, trim(kind_name), dist, mo_spread
2340  END IF
2341 
2342  END IF
2343  END DO !ispin
2344  END DO !imo
2345 
2346  IF (output_unit > 0) THEN
2347  WRITE (unit=output_unit, fmt="(T3,A,/)") &
2348  "---------------------------------------------------------------------------------"
2349  END IF
2350 
2351 ! Go back to initial BCs
2352  cell%perd = perd_init
2353 
2354  END SUBROUTINE write_mos_to_ex_atoms_association
2355 
2356 ! **************************************************************************************************
2357 !> \brief Performs Mulliken population analysis for the MO(s) of a donor_state_type so that user
2358 !> can verify it is indeed a core state
2359 !> \param donor_state ...
2360 !> \param qs_env ...
2361 !> \note This is a specific case of Mulliken analysis. In general one computes sum_i (SP)_ii, where
2362 !> i labels the basis function centered on the atom of interest. For a specific MO with index
2363 !> j, one need to compute sum_{ik} c_{ij} S_{ik} c_{kj}, k = 1,nao
2364 ! **************************************************************************************************
2365  SUBROUTINE perform_mulliken_on_donor_state(donor_state, qs_env)
2366  TYPE(donor_state_type), POINTER :: donor_state
2367  TYPE(qs_environment_type), POINTER :: qs_env
2368 
2369  INTEGER :: at_index, i, ispin, nao, natom, ndo_mo, &
2370  ndo_so, nsgf, nspins, output_unit
2371  INTEGER, DIMENSION(:), POINTER :: first_sgf, last_sgf
2372  INTEGER, DIMENSION(:, :), POINTER :: mo_indices
2373  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: mul_pop, pop_mat
2374  REAL(dp), DIMENSION(:, :), POINTER :: work_array
2375  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2376  TYPE(cp_fm_struct_type), POINTER :: col_vect_struct
2377  TYPE(cp_fm_type) :: work_vect
2378  TYPE(cp_fm_type), POINTER :: gs_coeffs
2379  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2380  TYPE(mp_para_env_type), POINTER :: para_env
2381  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2382  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2383 
2384  NULLIFY (mo_indices, qs_kind_set, particle_set, first_sgf, work_array)
2385  NULLIFY (matrix_s, para_env, blacs_env, col_vect_struct, last_sgf)
2386 
2387 ! Initialization
2388  at_index = donor_state%at_index
2389  mo_indices => donor_state%mo_indices
2390  ndo_mo = donor_state%ndo_mo
2391  gs_coeffs => donor_state%gs_coeffs
2392  output_unit = cp_logger_get_default_io_unit()
2393  nspins = 1; IF (SIZE(mo_indices, 2) == 2) nspins = 2
2394  ndo_so = ndo_mo*nspins
2395  ALLOCATE (mul_pop(ndo_mo, nspins))
2396  mul_pop = 0.0_dp
2397 
2398  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set, &
2399  para_env=para_env, blacs_env=blacs_env, matrix_s=matrix_s)
2400  CALL cp_fm_get_info(gs_coeffs, nrow_global=nao, matrix_struct=col_vect_struct)
2401 
2402  natom = SIZE(particle_set, 1)
2403  ALLOCATE (first_sgf(natom))
2404  ALLOCATE (last_sgf(natom))
2405 
2406  CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
2407  nsgf = last_sgf(at_index) - first_sgf(at_index) + 1
2408 
2409  CALL cp_fm_create(work_vect, col_vect_struct)
2410 
2411 ! Take the product of S*coeffs
2412  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, gs_coeffs, work_vect, ncol=ndo_so)
2413 
2414 ! Only consider the product coeffs^T * S * coeffs on the atom of interest
2415  ALLOCATE (work_array(nsgf, ndo_so))
2416  ALLOCATE (pop_mat(ndo_so, ndo_so))
2417 
2418  CALL cp_fm_get_submatrix(fm=work_vect, target_m=work_array, start_row=first_sgf(at_index), &
2419  start_col=1, n_rows=nsgf, n_cols=ndo_so, transpose=.false.)
2420 
2421  CALL dgemm('T', 'N', ndo_so, ndo_so, nsgf, 1.0_dp, donor_state%contract_coeffs, nsgf, &
2422  work_array, nsgf, 0.0_dp, pop_mat, ndo_so)
2423 
2424 ! The Mullikan population for the MOs in on the diagonal.
2425  DO ispin = 1, nspins
2426  DO i = 1, ndo_mo
2427  mul_pop(i, ispin) = pop_mat((ispin - 1)*ndo_mo + i, (ispin - 1)*ndo_mo + i)
2428  END DO
2429  END DO
2430 
2431 ! Printing in main output file
2432  IF (output_unit > 0) THEN
2433  WRITE (unit=output_unit, fmt="(T5,A,/,T5,A)") &
2434  "Mulliken population analysis retricted to the associated MO(s) yields: ", &
2435  " Spin MO index charge"
2436  DO ispin = 1, nspins
2437  DO i = 1, ndo_mo
2438  WRITE (unit=output_unit, fmt="(T51,I4,I10,F11.3)") &
2439  ispin, mo_indices(i, ispin), mul_pop(i, ispin)
2440  END DO
2441  END DO
2442  END IF
2443 
2444 ! Clean-up
2445  DEALLOCATE (first_sgf, last_sgf, work_array)
2446  CALL cp_fm_release(work_vect)
2447 
2448  END SUBROUTINE perform_mulliken_on_donor_state
2449 
2450 ! **************************************************************************************************
2451 !> \brief write the PDOS wrt the LR-orbitals for the current donor_state and/or the CUBES files
2452 !> \param ex_type the excitation type: singlet, triplet, spin-conserving, etc
2453 !> \param donor_state ...
2454 !> \param xas_tdp_env ...
2455 !> \param xas_tdp_section ...
2456 !> \param qs_env ...
2457 ! **************************************************************************************************
2458  SUBROUTINE xas_tdp_post(ex_type, donor_state, xas_tdp_env, xas_tdp_section, qs_env)
2459 
2460  INTEGER, INTENT(IN) :: ex_type
2461  TYPE(donor_state_type), POINTER :: donor_state
2462  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2463  TYPE(section_vals_type), POINTER :: xas_tdp_section
2464  TYPE(qs_environment_type), POINTER :: qs_env
2465 
2466  CHARACTER(len=*), PARAMETER :: routinen = 'xas_tdp_post'
2467 
2468  CHARACTER(len=default_string_length) :: domo, domon, excite, pos, xas_mittle
2469  INTEGER :: ex_state_idx, handle, ic, ido_mo, imo, irep, ispin, n_dependent, n_rep, nao, &
2470  ncubes, ndo_mo, ndo_so, nlumo, nmo, nspins, output_unit
2471  INTEGER, DIMENSION(:), POINTER :: bounds, list, state_list
2472  LOGICAL :: append_cube, do_cubes, do_pdos, &
2473  do_wfn_restart
2474  REAL(dp), DIMENSION(:), POINTER :: lr_evals
2475  REAL(dp), DIMENSION(:, :), POINTER :: centers
2476  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2477  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2478  TYPE(cp_fm_struct_type), POINTER :: fm_struct, mo_struct
2479  TYPE(cp_fm_type) :: mo_coeff, work_fm
2480  TYPE(cp_fm_type), POINTER :: lr_coeffs
2481  TYPE(cp_logger_type), POINTER :: logger
2482  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2483  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2484  TYPE(mo_set_type), POINTER :: mo_set
2485  TYPE(mp_para_env_type), POINTER :: para_env
2486  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2487  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2488  TYPE(section_vals_type), POINTER :: print_key
2489 
2490  NULLIFY (atomic_kind_set, particle_set, qs_kind_set, mo_set, lr_evals, lr_coeffs)
2491  NULLIFY (mo_struct, para_env, blacs_env, fm_struct, matrix_s, print_key, logger)
2492  NULLIFY (bounds, state_list, list, mos)
2493 
2494  !Tests on what to do
2495  logger => cp_get_default_logger()
2496  do_pdos = .false.; do_cubes = .false.; do_wfn_restart = .false.
2497 
2498  IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
2499  "PRINT%PDOS"), cp_p_file)) do_pdos = .true.
2500 
2501  IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
2502  "PRINT%CUBES"), cp_p_file)) do_cubes = .true.
2503 
2504  IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
2505  "PRINT%RESTART_WFN"), cp_p_file)) do_wfn_restart = .true.
2506 
2507  IF (.NOT. (do_pdos .OR. do_cubes .OR. do_wfn_restart)) RETURN
2508 
2509  CALL timeset(routinen, handle)
2510 
2511  !Initialization
2512  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
2513  qs_kind_set=qs_kind_set, para_env=para_env, blacs_env=blacs_env, &
2514  matrix_s=matrix_s, mos=mos)
2515 
2516  SELECT CASE (ex_type)
2517  CASE (tddfpt_spin_cons)
2518  lr_evals => donor_state%sc_evals
2519  lr_coeffs => donor_state%sc_coeffs
2520  nspins = 2
2521  excite = "spincons"
2522  CASE (tddfpt_spin_flip)
2523  lr_evals => donor_state%sf_evals
2524  lr_coeffs => donor_state%sf_coeffs
2525  nspins = 2
2526  excite = "spinflip"
2527  CASE (tddfpt_singlet)
2528  lr_evals => donor_state%sg_evals
2529  lr_coeffs => donor_state%sg_coeffs
2530  nspins = 1
2531  excite = "singlet"
2532  CASE (tddfpt_triplet)
2533  lr_evals => donor_state%tp_evals
2534  lr_coeffs => donor_state%tp_coeffs
2535  nspins = 1
2536  excite = "triplet"
2537  END SELECT
2538 
2539  SELECT CASE (donor_state%state_type)
2540  CASE (xas_1s_type)
2541  domo = "1s"
2542  CASE (xas_2s_type)
2543  domo = "2s"
2544  CASE (xas_2p_type)
2545  domo = "2p"
2546  END SELECT
2547 
2548  ndo_mo = donor_state%ndo_mo
2549  ndo_so = ndo_mo*nspins
2550  nmo = SIZE(lr_evals)
2551  CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
2552 
2553  CALL cp_fm_struct_create(mo_struct, context=blacs_env, para_env=para_env, &
2554  nrow_global=nao, ncol_global=nmo)
2555  CALL cp_fm_create(mo_coeff, mo_struct)
2556 
2557  !Dump the TDDFT excited state AMEW wavefunction into a file for restart in RTP
2558  IF (do_wfn_restart) THEN
2559  block
2560  TYPE(mo_set_type), DIMENSION(2) :: restart_mos
2561  IF (.NOT. (nspins == 1 .AND. donor_state%state_type == xas_1s_type)) THEN
2562  cpabort("RESTART.wfn file only available for RKS K-edge XAS spectroscopy")
2563  END IF
2564 
2565  CALL section_vals_val_get(xas_tdp_section, "PRINT%RESTART_WFN%EXCITED_STATE_INDEX", n_rep_val=n_rep)
2566 
2567  DO irep = 1, n_rep
2568  CALL section_vals_val_get(xas_tdp_section, "PRINT%RESTART_WFN%EXCITED_STATE_INDEX", &
2569  i_rep_val=irep, i_val=ex_state_idx)
2570  cpassert(ex_state_idx <= SIZE(lr_evals))
2571 
2572  DO ispin = 1, 2
2573  CALL duplicate_mo_set(restart_mos(ispin), mos(1))
2574  ! Set the new occupation number in the case of spin-independent based calculaltion since the restart is spin-depedent
2575  IF (SIZE(mos) == 1) &
2576  restart_mos(ispin)%occupation_numbers = mos(1)%occupation_numbers/2
2577  END DO
2578 
2579  CALL cp_fm_to_fm_submat(msource=lr_coeffs, mtarget=restart_mos(1)%mo_coeff, nrow=nao, &
2580  ncol=1, s_firstrow=1, s_firstcol=ex_state_idx, t_firstrow=1, &
2581  t_firstcol=donor_state%mo_indices(1, 1))
2582 
2583  xas_mittle = 'xasat'//trim(adjustl(cp_to_string(donor_state%at_index)))//'_'//trim(domo)// &
2584  '_'//trim(excite)//'_idx'//trim(adjustl(cp_to_string(ex_state_idx)))
2585  output_unit = cp_print_key_unit_nr(logger, xas_tdp_section, "PRINT%RESTART_WFN", &
2586  extension=".wfn", file_status="REPLACE", &
2587  file_action="WRITE", file_form="UNFORMATTED", &
2588  middle_name=xas_mittle)
2589 
2590  CALL write_mo_set_low(restart_mos, particle_set=particle_set, &
2591  qs_kind_set=qs_kind_set, ires=output_unit)
2592 
2593  CALL cp_print_key_finished_output(output_unit, logger, xas_tdp_section, "PRINT%RESTART_WFN")
2594 
2595  DO ispin = 1, 2
2596  CALL deallocate_mo_set(restart_mos(ispin))
2597  END DO
2598  END DO
2599  END block
2600  END IF
2601 
2602  !PDOS related stuff
2603  IF (do_pdos) THEN
2604 
2605  !If S^0.5 not yet stored, compute it once and for all
2606  IF (.NOT. ASSOCIATED(xas_tdp_env%matrix_shalf) .AND. do_pdos) THEN
2607  CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
2608  nrow_global=nao, ncol_global=nao)
2609  ALLOCATE (xas_tdp_env%matrix_shalf)
2610  CALL cp_fm_create(xas_tdp_env%matrix_shalf, fm_struct)
2611  CALL cp_fm_create(work_fm, fm_struct)
2612 
2613  CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, xas_tdp_env%matrix_shalf)
2614  CALL cp_fm_power(xas_tdp_env%matrix_shalf, work_fm, 0.5_dp, epsilon(0.0_dp), n_dependent)
2615 
2616  CALL cp_fm_release(work_fm)
2617  CALL cp_fm_struct_release(fm_struct)
2618  END IF
2619 
2620  !Giving some PDOS info
2621  output_unit = cp_logger_get_default_io_unit()
2622  IF (output_unit > 0) THEN
2623  WRITE (unit=output_unit, fmt="(/,T5,A,/,T5,A,/,T5,A)") &
2624  "Computing the PDOS of linear-response orbitals for spectral features analysis", &
2625  "Note: using standard PDOS routines => ignore mentions of KS states and MO ", &
2626  " occupation numbers. Eigenvalues in *.pdos files are excitations energies."
2627  END IF
2628 
2629  !Check on NLUMO
2630  CALL section_vals_val_get(xas_tdp_section, "PRINT%PDOS%NLUMO", i_val=nlumo)
2631  IF (nlumo .NE. 0) THEN
2632  cpwarn("NLUMO is irrelevant for XAS_TDP PDOS. It was overwritten to 0.")
2633  END IF
2634  CALL section_vals_val_set(xas_tdp_section, "PRINT%PDOS%NLUMO", i_val=0)
2635  END IF
2636 
2637  !CUBES related stuff
2638  IF (do_cubes) THEN
2639 
2640  print_key => section_vals_get_subs_vals(xas_tdp_section, "PRINT%CUBES")
2641 
2642  CALL section_vals_val_get(print_key, "CUBES_LU_BOUNDS", i_vals=bounds)
2643  ncubes = bounds(2) - bounds(1) + 1
2644  IF (ncubes > 0) THEN
2645  ALLOCATE (state_list(ncubes))
2646  DO ic = 1, ncubes
2647  state_list(ic) = bounds(1) + ic - 1
2648  END DO
2649  END IF
2650 
2651  IF (.NOT. ASSOCIATED(state_list)) THEN
2652  CALL section_vals_val_get(print_key, "CUBES_LIST", n_rep_val=n_rep)
2653 
2654  ncubes = 0
2655  DO irep = 1, n_rep
2656  NULLIFY (list)
2657  CALL section_vals_val_get(print_key, "CUBES_LIST", i_rep_val=irep, i_vals=list)
2658  IF (ASSOCIATED(list)) THEN
2659  CALL reallocate(state_list, 1, ncubes + SIZE(list))
2660  DO ic = 1, SIZE(list)
2661  state_list(ncubes + ic) = list(ic)
2662  END DO
2663  ncubes = ncubes + SIZE(list)
2664  END IF
2665  END DO
2666  END IF
2667 
2668  IF (.NOT. ASSOCIATED(state_list)) THEN
2669  ncubes = 1
2670  ALLOCATE (state_list(1))
2671  state_list(1) = 1
2672  END IF
2673 
2674  CALL section_vals_val_get(print_key, "APPEND", l_val=append_cube)
2675  pos = "REWIND"
2676  IF (append_cube) pos = "APPEND"
2677 
2678  ALLOCATE (centers(6, ncubes))
2679  centers = 0.0_dp
2680 
2681  END IF
2682 
2683  !Loop over MOs and spin, one PDOS/CUBE for each
2684  DO ido_mo = 1, ndo_mo
2685  DO ispin = 1, nspins
2686 
2687  !need to create a mo set for the LR-orbitals
2688  ALLOCATE (mo_set)
2689  CALL allocate_mo_set(mo_set, nao=nao, nmo=nmo, nelectron=nmo, n_el_f=real(nmo, dp), &
2690  maxocc=1.0_dp, flexible_electron_count=0.0_dp)
2691  CALL init_mo_set(mo_set, fm_ref=mo_coeff, name="PDOS XAS_TDP MOs")
2692  mo_set%eigenvalues(:) = lr_evals(:)
2693 
2694  !get the actual coeff => most common case: closed-shell K-edge, can directly take lr_coeffs
2695  IF (nspins == 1 .AND. ndo_mo == 1) THEN
2696  CALL cp_fm_to_fm(lr_coeffs, mo_set%mo_coeff)
2697  ELSE
2698  DO imo = 1, nmo
2699  CALL cp_fm_to_fm_submat(msource=lr_coeffs, mtarget=mo_set%mo_coeff, &
2700  nrow=nao, ncol=1, s_firstrow=1, &
2701  s_firstcol=(imo - 1)*ndo_so + (ispin - 1)*ndo_mo + ido_mo, &
2702  t_firstrow=1, t_firstcol=imo)
2703  END DO
2704  END IF
2705 
2706  !naming the output
2707  domon = domo
2708  IF (donor_state%state_type == xas_2p_type) domon = trim(domo)//trim(adjustl(cp_to_string(ido_mo)))
2709  xas_mittle = 'xasat'//trim(adjustl(cp_to_string(donor_state%at_index)))//'_'// &
2710  trim(domon)//'_'//trim(excite)
2711 
2712  IF (do_pdos) THEN
2713  CALL calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
2714  qs_env, xas_tdp_section, ispin, xas_mittle, &
2715  external_matrix_shalf=xas_tdp_env%matrix_shalf)
2716  END IF
2717 
2718  IF (do_cubes) THEN
2719  CALL qs_print_cubes(qs_env, mo_set%mo_coeff, ncubes, state_list, centers, &
2720  print_key=print_key, root=xas_mittle, ispin=ispin, &
2721  file_position=pos)
2722  END IF
2723 
2724  !clean-up
2725  CALL deallocate_mo_set(mo_set)
2726  DEALLOCATE (mo_set)
2727 
2728  END DO
2729  END DO
2730 
2731  !clean-up
2732  CALL cp_fm_release(mo_coeff)
2733  CALL cp_fm_struct_release(mo_struct)
2734  IF (do_cubes) DEALLOCATE (centers, state_list)
2735 
2736  CALL timestop(handle)
2737 
2738  END SUBROUTINE xas_tdp_post
2739 
2740 ! **************************************************************************************************
2741 !> \brief Computed the LUMOs for the OT eigensolver guesses
2742 !> \param xas_tdp_env ...
2743 !> \param xas_tdp_control ...
2744 !> \param qs_env ...
2745 !> \note Uses stendard diagonalization. Do not use the stendard make_lumo subroutine as it uses
2746 !> the OT eigensolver and there is no guarantee that it will converge fast
2747 ! **************************************************************************************************
2748  SUBROUTINE make_lumo_guess(xas_tdp_env, xas_tdp_control, qs_env)
2749 
2750  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2751  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2752  TYPE(qs_environment_type), POINTER :: qs_env
2753 
2754  CHARACTER(len=*), PARAMETER :: routinen = 'make_lumo_guess'
2755 
2756  INTEGER :: handle, ispin, nao, nelec_spin(2), &
2757  nlumo(2), nocc(2), nspins
2758  LOGICAL :: do_os
2759  REAL(dp), ALLOCATABLE, DIMENSION(:) :: evals
2760  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2761  TYPE(cp_fm_struct_type), POINTER :: fm_struct, lumo_struct
2762  TYPE(cp_fm_type) :: amatrix, bmatrix, evecs, work_fm
2763  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
2764  TYPE(mp_para_env_type), POINTER :: para_env
2765 
2766  NULLIFY (matrix_ks, matrix_s, para_env, blacs_env)
2767  NULLIFY (lumo_struct, fm_struct)
2768 
2769  CALL timeset(routinen, handle)
2770 
2771  do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
2772  nspins = 1; IF (do_os) nspins = 2
2773  ALLOCATE (xas_tdp_env%lumo_evecs(nspins))
2774  ALLOCATE (xas_tdp_env%lumo_evals(nspins))
2775  CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s, nelectron_spin=nelec_spin, &
2776  para_env=para_env, blacs_env=blacs_env)
2777  CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
2778 
2779  IF (do_os) THEN
2780  nlumo = nao - nelec_spin
2781  nocc = nelec_spin
2782  ELSE
2783  nlumo = nao - nelec_spin(1)/2
2784  nocc = nelec_spin(1)/2
2785  END IF
2786 
2787  ALLOCATE (xas_tdp_env%ot_prec(nspins))
2788 
2789  DO ispin = 1, nspins
2790 
2791  !Going through fm to diagonalize
2792  CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2793  nrow_global=nao, ncol_global=nao)
2794  CALL cp_fm_create(amatrix, fm_struct)
2795  CALL cp_fm_create(bmatrix, fm_struct)
2796  CALL cp_fm_create(evecs, fm_struct)
2797  CALL cp_fm_create(work_fm, fm_struct)
2798  ALLOCATE (evals(nao))
2799  ALLOCATE (xas_tdp_env%lumo_evals(ispin)%array(nlumo(ispin)))
2800 
2801  CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, amatrix)
2802  CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, bmatrix)
2803 
2804  !The actual diagonalization through Cholesky decomposition
2805  CALL cp_fm_geeig(amatrix, bmatrix, evecs, evals, work_fm)
2806 
2807  !Storing results
2808  CALL cp_fm_struct_create(lumo_struct, para_env=para_env, context=blacs_env, &
2809  nrow_global=nao, ncol_global=nlumo(ispin))
2810  CALL cp_fm_create(xas_tdp_env%lumo_evecs(ispin), lumo_struct)
2811 
2812  CALL cp_fm_to_fm_submat(evecs, xas_tdp_env%lumo_evecs(ispin), nrow=nao, &
2813  ncol=nlumo(ispin), s_firstrow=1, s_firstcol=nocc(ispin) + 1, &
2814  t_firstrow=1, t_firstcol=1)
2815 
2816  xas_tdp_env%lumo_evals(ispin)%array(1:nlumo(ispin)) = evals(nocc(ispin) + 1:nao)
2817 
2818  CALL build_ot_spin_prec(evecs, evals, ispin, xas_tdp_env, xas_tdp_control, qs_env)
2819 
2820  !clean-up
2821  CALL cp_fm_release(amatrix)
2822  CALL cp_fm_release(bmatrix)
2823  CALL cp_fm_release(evecs)
2824  CALL cp_fm_release(work_fm)
2825  CALL cp_fm_struct_release(fm_struct)
2826  CALL cp_fm_struct_release(lumo_struct)
2827  DEALLOCATE (evals)
2828  END DO
2829 
2830  CALL timestop(handle)
2831 
2832  END SUBROUTINE make_lumo_guess
2833 
2834 ! **************************************************************************************************
2835 !> \brief Builds a preconditioner for the OT eigensolver, based on some heurstics that prioritize
2836 !> LUMOs with lower eigenvalues
2837 !> \param evecs all the ground state eigenvectors
2838 !> \param evals all the ground state eigenvalues
2839 !> \param ispin ...
2840 !> \param xas_tdp_env ...
2841 !> \param xas_tdp_control ...
2842 !> \param qs_env ...
2843 !> \note assumes that the preconditioner matrix array is allocated
2844 ! **************************************************************************************************
2845  SUBROUTINE build_ot_spin_prec(evecs, evals, ispin, xas_tdp_env, xas_tdp_control, qs_env)
2846 
2847  TYPE(cp_fm_type), INTENT(IN) :: evecs
2848  REAL(dp), DIMENSION(:) :: evals
2849  INTEGER :: ispin
2850  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2851  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2852  TYPE(qs_environment_type), POINTER :: qs_env
2853 
2854  CHARACTER(len=*), PARAMETER :: routinen = 'build_ot_spin_prec'
2855 
2856  INTEGER :: handle, nao, nelec_spin(2), nguess, &
2857  nocc, nspins
2858  LOGICAL :: do_os
2859  REAL(dp) :: shift
2860  REAL(dp), ALLOCATABLE, DIMENSION(:) :: scaling
2861  TYPE(cp_fm_struct_type), POINTER :: fm_struct
2862  TYPE(cp_fm_type) :: fm_prec, work_fm
2863  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2864  TYPE(mp_para_env_type), POINTER :: para_env
2865 
2866  NULLIFY (fm_struct, para_env, matrix_s)
2867 
2868  CALL timeset(routinen, handle)
2869 
2870  do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
2871  CALL get_qs_env(qs_env, para_env=para_env, nelectron_spin=nelec_spin, matrix_s=matrix_s)
2872  CALL cp_fm_get_info(evecs, nrow_global=nao, matrix_struct=fm_struct)
2873  CALL cp_fm_create(fm_prec, fm_struct)
2874  ALLOCATE (scaling(nao))
2875  nocc = nelec_spin(1)/2
2876  nspins = 1
2877  IF (do_os) THEN
2878  nocc = nelec_spin(ispin)
2879  nspins = 2
2880  END IF
2881 
2882  !rough estimate of the number of required evals
2883  nguess = nao - nocc
2884  IF (xas_tdp_control%n_excited > 0 .AND. xas_tdp_control%n_excited < nguess) THEN
2885  nguess = xas_tdp_control%n_excited/nspins
2886  ELSE IF (xas_tdp_control%e_range > 0.0_dp) THEN
2887  nguess = count(evals(nocc + 1:nao) - evals(nocc + 1) .LE. xas_tdp_control%e_range)
2888  END IF
2889 
2890  !Give max weight to the first LUMOs
2891  scaling(nocc + 1:nocc + nguess) = 100.0_dp
2892  !Then gradually decrease weight
2893  shift = evals(nocc + 1) - 0.01_dp
2894  scaling(nocc + nguess:nao) = 1.0_dp/(evals(nocc + nguess:nao) - shift)
2895  !HOMOs do not matter, but need well behaved matrix
2896  scaling(1:nocc) = 1.0_dp
2897 
2898  !Building the precond as an fm
2899  CALL cp_fm_create(work_fm, fm_struct)
2900 
2901  CALL cp_fm_copy_general(evecs, work_fm, para_env)
2902  CALL cp_fm_column_scale(work_fm, scaling)
2903 
2904  CALL parallel_gemm('N', 'T', nao, nao, nao, 1.0_dp, work_fm, evecs, 0.0_dp, fm_prec)
2905 
2906  !Copy into dbcsr format
2907  ALLOCATE (xas_tdp_env%ot_prec(ispin)%matrix)
2908  CALL dbcsr_create(xas_tdp_env%ot_prec(ispin)%matrix, template=matrix_s(1)%matrix, name="OT_PREC")
2909  CALL copy_fm_to_dbcsr(fm_prec, xas_tdp_env%ot_prec(ispin)%matrix)
2910  CALL dbcsr_filter(xas_tdp_env%ot_prec(ispin)%matrix, xas_tdp_control%eps_filter)
2911 
2912  CALL cp_fm_release(work_fm)
2913  CALL cp_fm_release(fm_prec)
2914 
2915  CALL timestop(handle)
2916 
2917  END SUBROUTINE build_ot_spin_prec
2918 
2919 ! **************************************************************************************************
2920 !> \brief Prints GW2X corrected ionization potentials to main output file, including SOC splitting
2921 !> \param donor_state ...
2922 !> \param xas_tdp_env ...
2923 !> \param xas_tdp_control ...
2924 !> \param qs_env ...
2925 ! **************************************************************************************************
2926  SUBROUTINE print_xps(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
2927 
2928  TYPE(donor_state_type), POINTER :: donor_state
2929  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2930  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2931  TYPE(qs_environment_type), POINTER :: qs_env
2932 
2933  INTEGER :: ido_mo, ispin, nspins, output_unit
2934  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ips, soc_shifts
2935 
2936  output_unit = cp_logger_get_default_io_unit()
2937 
2938  nspins = 1; IF (xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks) nspins = 2
2939 
2940  ALLOCATE (ips(SIZE(donor_state%gw2x_evals, 1), SIZE(donor_state%gw2x_evals, 2)))
2941  ips(:, :) = donor_state%gw2x_evals(:, :)
2942 
2943  !IPs in PBCs cannot be trusted because of a lack of a potential reference
2944  IF (.NOT. xas_tdp_control%is_periodic) THEN
2945 
2946  !Apply SOC splitting
2947  IF (donor_state%ndo_mo > 1) THEN
2948  CALL get_soc_splitting(soc_shifts, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
2949  ips(:, :) = ips(:, :) + soc_shifts
2950 
2951  IF (output_unit > 0) THEN
2952  WRITE (output_unit, fmt="(/,T5,A,F23.6)") &
2953  "Ionization potentials for XPS (GW2X + SOC): ", -ips(1, 1)*evolt
2954 
2955  DO ispin = 1, nspins
2956  DO ido_mo = 1, donor_state%ndo_mo
2957 
2958  IF (ispin == 1 .AND. ido_mo == 1) cycle
2959 
2960  WRITE (output_unit, fmt="(T5,A,F23.6)") &
2961  " ", -ips(ido_mo, ispin)*evolt
2962 
2963  END DO
2964  END DO
2965  END IF
2966 
2967  ELSE
2968 
2969  ! No SOC, only 1 donor MO per spin
2970  IF (output_unit > 0) THEN
2971  WRITE (output_unit, fmt="(/,T5,A,F29.6)") &
2972  "Ionization potentials for XPS (GW2X): ", -ips(1, 1)*evolt
2973 
2974  IF (nspins == 2) THEN
2975  WRITE (output_unit, fmt="(T5,A,F29.6)") &
2976  " ", -ips(1, 2)*evolt
2977  END IF
2978  END IF
2979 
2980  END IF
2981  END IF
2982 
2983  END SUBROUTINE print_xps
2984 
2985 ! **************************************************************************************************
2986 !> \brief Prints the excitation energies and the oscillator strengths for a given donor_state in a file
2987 !> \param donor_state the donor_state to print
2988 !> \param xas_tdp_env ...
2989 !> \param xas_tdp_control ...
2990 !> \param xas_tdp_section ...
2991 ! **************************************************************************************************
2992  SUBROUTINE print_xas_tdp_to_file(donor_state, xas_tdp_env, xas_tdp_control, xas_tdp_section)
2993 
2994  TYPE(donor_state_type), POINTER :: donor_state
2995  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2996  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2997  TYPE(section_vals_type), POINTER :: xas_tdp_section
2998 
2999  INTEGER :: i, output_unit, xas_tdp_unit
3000  TYPE(cp_logger_type), POINTER :: logger
3001 
3002  NULLIFY (logger)
3003  logger => cp_get_default_logger()
3004 
3005  xas_tdp_unit = cp_print_key_unit_nr(logger, xas_tdp_section, "PRINT%SPECTRUM", &
3006  extension=".spectrum", file_position="APPEND", &
3007  file_action="WRITE", file_form="FORMATTED")
3008 
3009  output_unit = cp_logger_get_default_io_unit()
3010 
3011  IF (output_unit > 0) THEN
3012  WRITE (output_unit, fmt="(/,T5,A,/)") &
3013  "Calculations done: "
3014  END IF
3015 
3016  IF (xas_tdp_control%do_spin_cons) THEN
3017  IF (xas_tdp_unit > 0) THEN
3018 
3019 ! Printing the general donor state information
3020  WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3021  "==================================================================================", &
3022  "XAS TDP open-shell spin-conserving (no SOC) excitations for DONOR STATE: ", &
3023  xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3024  "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3025  donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3026  "=================================================================================="
3027 
3028 ! Simply dump the excitation energies/ oscillator strength as they come
3029 
3030  IF (xas_tdp_control%do_quad) THEN
3031  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3032  " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3033  DO i = 1, SIZE(donor_state%sc_evals)
3034  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3035  i, donor_state%sc_evals(i)*evolt, donor_state%osc_str(i, 4), &
3036  donor_state%quad_osc_str(i)
3037  END DO
3038  ELSE IF (xas_tdp_control%xyz_dip) THEN
3039  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3040  " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3041  DO i = 1, SIZE(donor_state%sc_evals)
3042  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3043  i, donor_state%sc_evals(i)*evolt, donor_state%osc_str(i, 4), &
3044  donor_state%osc_str(i, 1), donor_state%osc_str(i, 2), donor_state%osc_str(i, 3)
3045  END DO
3046  ELSE
3047  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3048  " Index Excitation energy (eV) fosc dipole (a.u.)"
3049  DO i = 1, SIZE(donor_state%sc_evals)
3050  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3051  i, donor_state%sc_evals(i)*evolt, donor_state%osc_str(i, 4)
3052  END DO
3053  END IF
3054 
3055  WRITE (xas_tdp_unit, fmt="(A,/)") " "
3056  END IF !xas_tdp_unit > 0
3057 
3058  IF (output_unit > 0) THEN
3059  WRITE (output_unit, fmt="(T5,A,F17.6)") &
3060  "First spin-conserving XAS excitation energy (eV): ", donor_state%sc_evals(1)*evolt
3061  END IF
3062 
3063  END IF ! do_spin_cons
3064 
3065  IF (xas_tdp_control%do_spin_flip) THEN
3066  IF (xas_tdp_unit > 0) THEN
3067 
3068 ! Printing the general donor state information
3069  WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3070  "==================================================================================", &
3071  "XAS TDP open-shell spin-flip (no SOC) excitations for DONOR STATE: ", &
3072  xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3073  "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3074  donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3075  "=================================================================================="
3076 
3077 ! Simply dump the excitation energies/ oscillator strength as they come
3078 
3079  IF (xas_tdp_control%do_quad) THEN
3080  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3081  " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3082  DO i = 1, SIZE(donor_state%sf_evals)
3083  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3084  i, donor_state%sf_evals(i)*evolt, 0.0_dp, 0.0_dp !spin-forbidden !
3085  END DO
3086  ELSE IF (xas_tdp_control%xyz_dip) THEN
3087  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3088  " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3089  DO i = 1, SIZE(donor_state%sf_evals)
3090  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3091  i, donor_state%sf_evals(i)*evolt, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp
3092  END DO
3093  ELSE
3094  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3095  " Index Excitation energy (eV) fosc dipole (a.u.)"
3096  DO i = 1, SIZE(donor_state%sf_evals)
3097  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3098  i, donor_state%sf_evals(i)*evolt, 0.0_dp
3099  END DO
3100  END IF
3101 
3102  WRITE (xas_tdp_unit, fmt="(A,/)") " "
3103  END IF !xas_tdp_unit
3104 
3105  IF (output_unit > 0) THEN
3106  WRITE (output_unit, fmt="(T5,A,F23.6)") &
3107  "First spin-flip XAS excitation energy (eV): ", donor_state%sf_evals(1)*evolt
3108  END IF
3109  END IF ! do_spin_flip
3110 
3111  IF (xas_tdp_control%do_singlet) THEN
3112  IF (xas_tdp_unit > 0) THEN
3113 
3114 ! Printing the general donor state information
3115  WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3116  "==================================================================================", &
3117  "XAS TDP singlet excitations (no SOC) for DONOR STATE: ", &
3118  xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3119  "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3120  donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3121  "=================================================================================="
3122 
3123 ! Simply dump the excitation energies/ oscillator strength as they come
3124 
3125  IF (xas_tdp_control%do_quad) THEN
3126  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3127  " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3128  DO i = 1, SIZE(donor_state%sg_evals)
3129  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3130  i, donor_state%sg_evals(i)*evolt, donor_state%osc_str(i, 4), &
3131  donor_state%quad_osc_str(i)
3132  END DO
3133  ELSE IF (xas_tdp_control%xyz_dip) THEN
3134  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3135  " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3136  DO i = 1, SIZE(donor_state%sg_evals)
3137  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3138  i, donor_state%sg_evals(i)*evolt, donor_state%osc_str(i, 4), &
3139  donor_state%osc_str(i, 1), donor_state%osc_str(i, 2), donor_state%osc_str(i, 3)
3140  END DO
3141  ELSE
3142  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3143  " Index Excitation energy (eV) fosc dipole (a.u.)"
3144  DO i = 1, SIZE(donor_state%sg_evals)
3145  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3146  i, donor_state%sg_evals(i)*evolt, donor_state%osc_str(i, 4)
3147  END DO
3148  END IF
3149 
3150  WRITE (xas_tdp_unit, fmt="(A,/)") " "
3151  END IF !xas_tdp_unit
3152 
3153  IF (output_unit > 0) THEN
3154  WRITE (output_unit, fmt="(T5,A,F25.6)") &
3155  "First singlet XAS excitation energy (eV): ", donor_state%sg_evals(1)*evolt
3156  END IF
3157  END IF ! do_singlet
3158 
3159  IF (xas_tdp_control%do_triplet) THEN
3160  IF (xas_tdp_unit > 0) THEN
3161 
3162 ! Printing the general donor state information
3163  WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3164  "==================================================================================", &
3165  "XAS TDP triplet excitations (no SOC) for DONOR STATE: ", &
3166  xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3167  "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3168  donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3169  "=================================================================================="
3170 
3171 ! Simply dump the excitation energies/ oscillator strength as they come
3172 
3173  IF (xas_tdp_control%do_quad) THEN
3174  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3175  " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3176  DO i = 1, SIZE(donor_state%tp_evals)
3177  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3178  i, donor_state%tp_evals(i)*evolt, 0.0_dp, 0.0_dp !spin-forbidden !
3179  END DO
3180  ELSE IF (xas_tdp_control%xyz_dip) THEN
3181  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3182  " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3183  DO i = 1, SIZE(donor_state%tp_evals)
3184  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3185  i, donor_state%tp_evals(i)*evolt, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp
3186  END DO
3187  ELSE
3188  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3189  " Index Excitation energy (eV) fosc dipole (a.u.)"
3190  DO i = 1, SIZE(donor_state%tp_evals)
3191  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3192  i, donor_state%tp_evals(i)*evolt, 0.0_dp
3193  END DO
3194  END IF
3195 
3196  WRITE (xas_tdp_unit, fmt="(A,/)") " "
3197  END IF !xas_tdp_unit
3198 
3199  IF (output_unit > 0) THEN
3200  WRITE (output_unit, fmt="(T5,A,F25.6)") &
3201  "First triplet XAS excitation energy (eV): ", donor_state%tp_evals(1)*evolt
3202  END IF
3203  END IF ! do_triplet
3204 
3205  IF (xas_tdp_control%do_soc .AND. donor_state%state_type == xas_2p_type) THEN
3206  IF (xas_tdp_unit > 0) THEN
3207 
3208 ! Printing the general donor state information
3209  WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3210  "==================================================================================", &
3211  "XAS TDP excitations after spin-orbit coupling for DONOR STATE: ", &
3212  xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3213  "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3214  donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3215  "=================================================================================="
3216 
3217 ! Simply dump the excitation energies/ oscillator strength as they come
3218  IF (xas_tdp_control%do_quad) THEN
3219  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3220  " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3221  DO i = 1, SIZE(donor_state%soc_evals)
3222  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3223  i, donor_state%soc_evals(i)*evolt, donor_state%soc_osc_str(i, 4), &
3224  donor_state%soc_quad_osc_str(i)
3225  END DO
3226  ELSE IF (xas_tdp_control%xyz_dip) THEN
3227  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3228  " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3229  DO i = 1, SIZE(donor_state%soc_evals)
3230  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3231  i, donor_state%soc_evals(i)*evolt, donor_state%soc_osc_str(i, 4), &
3232  donor_state%soc_osc_str(i, 1), donor_state%soc_osc_str(i, 2), donor_state%soc_osc_str(i, 3)
3233  END DO
3234  ELSE
3235  WRITE (xas_tdp_unit, fmt="(T3,A)") &
3236  " Index Excitation energy (eV) fosc dipole (a.u.)"
3237  DO i = 1, SIZE(donor_state%soc_evals)
3238  WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3239  i, donor_state%soc_evals(i)*evolt, donor_state%soc_osc_str(i, 4)
3240  END DO
3241  END IF
3242 
3243  WRITE (xas_tdp_unit, fmt="(A,/)") " "
3244  END IF !xas_tdp_unit
3245 
3246  IF (output_unit > 0) THEN
3247  WRITE (output_unit, fmt="(T5,A,F29.6)") &
3248  "First SOC XAS excitation energy (eV): ", donor_state%soc_evals(1)*evolt
3249  END IF
3250  END IF !do_soc
3251 
3252  CALL cp_print_key_finished_output(xas_tdp_unit, logger, xas_tdp_section, "PRINT%SPECTRUM")
3253 
3254  END SUBROUTINE print_xas_tdp_to_file
3255 
3256 ! **************************************************************************************************
3257 !> \brief Prints the donor_state and excitation_type info into a RESTART file for cheap PDOS and/or
3258 !> CUBE printing without expensive computation
3259 !> \param ex_type singlet, triplet, etc.
3260 !> \param donor_state ...
3261 !> \param xas_tdp_section ...
3262 !> \param qs_env ...
3263 ! **************************************************************************************************
3264  SUBROUTINE write_donor_state_restart(ex_type, donor_state, xas_tdp_section, qs_env)
3265 
3266  INTEGER, INTENT(IN) :: ex_type
3267  TYPE(donor_state_type), POINTER :: donor_state
3268  TYPE(section_vals_type), POINTER :: xas_tdp_section
3269  TYPE(qs_environment_type), POINTER :: qs_env
3270 
3271  CHARACTER(len=*), PARAMETER :: routinen = 'write_donor_state_restart'
3272 
3273  CHARACTER(len=default_path_length) :: filename
3274  CHARACTER(len=default_string_length) :: domo, excite, my_middle
3275  INTEGER :: ex_atom, handle, ispin, nao, ndo_mo, &
3276  nex, nspins, output_unit, rst_unit, &
3277  state_type
3278  INTEGER, DIMENSION(:, :), POINTER :: mo_indices
3279  LOGICAL :: do_print
3280  REAL(dp), DIMENSION(:), POINTER :: lr_evals
3281  TYPE(cp_fm_type), POINTER :: lr_coeffs
3282  TYPE(cp_logger_type), POINTER :: logger
3283  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3284  TYPE(section_vals_type), POINTER :: print_key
3285 
3286  NULLIFY (logger, lr_coeffs, lr_evals, print_key, mos)
3287 
3288  !Initialization
3289  logger => cp_get_default_logger()
3290  do_print = .false.
3291  IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
3292  "PRINT%RESTART", used_print_key=print_key), cp_p_file)) do_print = .true.
3293 
3294  IF (.NOT. do_print) RETURN
3295 
3296  CALL timeset(routinen, handle)
3297 
3298  output_unit = cp_logger_get_default_io_unit()
3299 
3300  !Get general info
3301  SELECT CASE (ex_type)
3302  CASE (tddfpt_spin_cons)
3303  lr_evals => donor_state%sc_evals
3304  lr_coeffs => donor_state%sc_coeffs
3305  excite = "spincons"
3306  nspins = 2
3307  CASE (tddfpt_spin_flip)
3308  lr_evals => donor_state%sf_evals
3309  lr_coeffs => donor_state%sf_coeffs
3310  excite = "spinflip"
3311  nspins = 2
3312  CASE (tddfpt_singlet)
3313  lr_evals => donor_state%sg_evals
3314  lr_coeffs => donor_state%sg_coeffs
3315  excite = "singlet"
3316  nspins = 1
3317  CASE (tddfpt_triplet)
3318  lr_evals => donor_state%tp_evals
3319  lr_coeffs => donor_state%tp_coeffs
3320  excite = "triplet"
3321  nspins = 1
3322  END SELECT
3323 
3324  SELECT CASE (donor_state%state_type)
3325  CASE (xas_1s_type)
3326  domo = "1s"
3327  CASE (xas_2s_type)
3328  domo = "2s"
3329  CASE (xas_2p_type)
3330  domo = "2p"
3331  END SELECT
3332 
3333  ndo_mo = donor_state%ndo_mo
3334  nex = SIZE(lr_evals)
3335  CALL cp_fm_get_info(lr_coeffs, nrow_global=nao)
3336  state_type = donor_state%state_type
3337  ex_atom = donor_state%at_index
3338  mo_indices => donor_state%mo_indices
3339 
3340  !Opening restart file
3341  rst_unit = -1
3342  my_middle = 'xasat'//trim(adjustl(cp_to_string(ex_atom)))//'_'//trim(domo)//'_'//trim(excite)
3343  rst_unit = cp_print_key_unit_nr(logger, xas_tdp_section, "PRINT%RESTART", extension=".rst", &
3344  file_status="REPLACE", file_action="WRITE", &
3345  file_form="UNFORMATTED", middle_name=trim(my_middle))
3346 
3347  filename = cp_print_key_generate_filename(logger, print_key, middle_name=trim(my_middle), &
3348  extension=".rst", my_local=.false.)
3349 
3350  IF (output_unit > 0) THEN
3351  WRITE (unit=output_unit, fmt="(/,T5,A,/T5,A,A,A)") &
3352  "Linear-response orbitals and excitation energies are written in: ", &
3353  '"', trim(filename), '"'
3354  END IF
3355 
3356  !Writing
3357  IF (rst_unit > 0) THEN
3358  WRITE (rst_unit) ex_atom, state_type, ndo_mo, ex_type
3359  WRITE (rst_unit) nao, nex, nspins
3360  WRITE (rst_unit) mo_indices(:, :)
3361  WRITE (rst_unit) lr_evals(:)
3362  END IF
3363  CALL cp_fm_write_unformatted(lr_coeffs, rst_unit)
3364 
3365  !The MOs as well (because the may have been localized)
3366  CALL get_qs_env(qs_env, mos=mos)
3367  DO ispin = 1, nspins
3368  CALL cp_fm_write_unformatted(mos(ispin)%mo_coeff, rst_unit)
3369  END DO
3370 
3371  !closing
3372  CALL cp_print_key_finished_output(rst_unit, logger, xas_tdp_section, "PRINT%RESTART")
3373 
3374  CALL timestop(handle)
3375 
3376  END SUBROUTINE write_donor_state_restart
3377 
3378 ! **************************************************************************************************
3379 !> \brief Reads donor_state info from a restart file
3380 !> \param donor_state the pre-allocated donor_state
3381 !> \param ex_type the excitations stored in this specific file
3382 !> \param filename the restart file to read from
3383 !> \param qs_env ...
3384 ! **************************************************************************************************
3385  SUBROUTINE read_donor_state_restart(donor_state, ex_type, filename, qs_env)
3386 
3387  TYPE(donor_state_type), POINTER :: donor_state
3388  INTEGER, INTENT(OUT) :: ex_type
3389  CHARACTER(len=*), INTENT(IN) :: filename
3390  TYPE(qs_environment_type), POINTER :: qs_env
3391 
3392  CHARACTER(len=*), PARAMETER :: routinen = 'read_donor_state_restart'
3393 
3394  INTEGER :: handle, ispin, nao, nex, nspins, &
3395  output_unit, read_params(7), rst_unit
3396  INTEGER, DIMENSION(:, :), POINTER :: mo_indices
3397  LOGICAL :: file_exists
3398  REAL(dp), DIMENSION(:), POINTER :: lr_evals
3399  TYPE(cp_blacs_env_type), POINTER :: blacs_env
3400  TYPE(cp_fm_struct_type), POINTER :: fm_struct
3401  TYPE(cp_fm_type), POINTER :: lr_coeffs
3402  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3403  TYPE(mp_comm_type) :: group
3404  TYPE(mp_para_env_type), POINTER :: para_env
3405 
3406  NULLIFY (lr_evals, lr_coeffs, para_env, fm_struct, blacs_env, mos)
3407 
3408  CALL timeset(routinen, handle)
3409 
3410  output_unit = cp_logger_get_default_io_unit()
3411  cpassert(ASSOCIATED(donor_state))
3412  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
3413  group = para_env
3414 
3415  file_exists = .false.
3416  rst_unit = -1
3417 
3418  IF (para_env%is_source()) THEN
3419 
3420  INQUIRE (file=filename, exist=file_exists)
3421  IF (.NOT. file_exists) cpabort("Trying to read non-existing XAS_TDP restart file")
3422 
3423  CALL open_file(file_name=trim(filename), file_action="READ", file_form="UNFORMATTED", &
3424  file_position="REWIND", file_status="OLD", unit_number=rst_unit)
3425  END IF
3426 
3427  IF (output_unit > 0) THEN
3428  WRITE (unit=output_unit, fmt="(/,T5,A,/,T5,A,A,A)") &
3429  "Reading linear-response orbitals and excitation energies from file: ", &
3430  '"', filename, '"'
3431  END IF
3432 
3433  !read general params
3434  IF (rst_unit > 0) THEN
3435  READ (rst_unit) read_params(1:4)
3436  READ (rst_unit) read_params(5:7)
3437  END IF
3438  CALL group%bcast(read_params)
3439  donor_state%at_index = read_params(1)
3440  donor_state%state_type = read_params(2)
3441  donor_state%ndo_mo = read_params(3)
3442  ex_type = read_params(4)
3443  nao = read_params(5)
3444  nex = read_params(6)
3445  nspins = read_params(7)
3446 
3447  ALLOCATE (mo_indices(donor_state%ndo_mo, nspins))
3448  IF (rst_unit > 0) THEN
3449  READ (rst_unit) mo_indices(1:donor_state%ndo_mo, 1:nspins)
3450  END IF
3451  CALL group%bcast(mo_indices)
3452  donor_state%mo_indices => mo_indices
3453 
3454  !read evals
3455  ALLOCATE (lr_evals(nex))
3456  IF (rst_unit > 0) READ (rst_unit) lr_evals(1:nex)
3457  CALL group%bcast(lr_evals)
3458 
3459  !read evecs
3460  CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
3461  nrow_global=nao, ncol_global=nex*donor_state%ndo_mo*nspins)
3462  ALLOCATE (lr_coeffs)
3463  CALL cp_fm_create(lr_coeffs, fm_struct)
3464  CALL cp_fm_read_unformatted(lr_coeffs, rst_unit)
3465  CALL cp_fm_struct_release(fm_struct)
3466 
3467  !read MO coeffs and replace in qs_env
3468  CALL get_qs_env(qs_env, mos=mos)
3469  DO ispin = 1, nspins
3470  CALL cp_fm_read_unformatted(mos(ispin)%mo_coeff, rst_unit)
3471  END DO
3472 
3473  !closing file
3474  IF (para_env%is_source()) THEN
3475  CALL close_file(unit_number=rst_unit)
3476  END IF
3477 
3478  !case study on excitation type
3479  SELECT CASE (ex_type)
3480  CASE (tddfpt_spin_cons)
3481  donor_state%sc_evals => lr_evals
3482  donor_state%sc_coeffs => lr_coeffs
3483  CASE (tddfpt_spin_flip)
3484  donor_state%sf_evals => lr_evals
3485  donor_state%sf_coeffs => lr_coeffs
3486  CASE (tddfpt_singlet)
3487  donor_state%sg_evals => lr_evals
3488  donor_state%sg_coeffs => lr_coeffs
3489  CASE (tddfpt_triplet)
3490  donor_state%tp_evals => lr_evals
3491  donor_state%tp_coeffs => lr_coeffs
3492  END SELECT
3493 
3494  CALL timestop(handle)
3495 
3496  END SUBROUTINE read_donor_state_restart
3497 
3498 ! **************************************************************************************************
3499 !> \brief Checks whether this is a restart calculation and runs it if so
3500 !> \param rst_filename the file to read for restart
3501 !> \param xas_tdp_section ...
3502 !> \param qs_env ...
3503 ! **************************************************************************************************
3504  SUBROUTINE restart_calculation(rst_filename, xas_tdp_section, qs_env)
3505 
3506  CHARACTER(len=*), INTENT(IN) :: rst_filename
3507  TYPE(section_vals_type), POINTER :: xas_tdp_section
3508  TYPE(qs_environment_type), POINTER :: qs_env
3509 
3510  INTEGER :: ex_type
3511  TYPE(donor_state_type), POINTER :: donor_state
3512  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
3513 
3514  NULLIFY (xas_tdp_env, donor_state)
3515 
3516  !create a donor_state that we fill with the information we read
3517  ALLOCATE (donor_state)
3518  CALL donor_state_create(donor_state)
3519  CALL read_donor_state_restart(donor_state, ex_type, rst_filename, qs_env)
3520 
3521  !create a dummy xas_tdp_env and compute the post XAS_TDP stuff
3522  CALL xas_tdp_env_create(xas_tdp_env)
3523  CALL xas_tdp_post(ex_type, donor_state, xas_tdp_env, xas_tdp_section, qs_env)
3524 
3525  !clean-up
3526  CALL xas_tdp_env_release(xas_tdp_env)
3527  CALL free_ds_memory(donor_state)
3528  DEALLOCATE (donor_state%mo_indices)
3529  DEALLOCATE (donor_state)
3530 
3531  END SUBROUTINE restart_calculation
3532 
3533 END MODULE xas_tdp_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
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
admm_utils
Contains methods used in the context of density fitting.
Definition: admm_utils.F:15
admm_utils::admm_uncorrect_for_eigenvalues
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition: admm_utils.F:127
admm_utils::admm_correct_for_eigenvalues
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition: admm_utils.F:53
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::deallocate_gto_basis_set
subroutine, public deallocate_gto_basis_set(gto_basis_set)
...
Definition: basis_set_types.F:168
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
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::bussy2021a
integer, save, public bussy2021a
Definition: bibliography.F:43
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
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_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_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
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_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1882
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition: cp_fm_diag.F:896
cp_fm_diag::cp_fm_syevd
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
Definition: cp_fm_diag.F:413
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition: cp_fm_diag.F:1274
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_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition: cp_fm_types.F:1538
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition: cp_fm_types.F:570
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition: cp_fm_types.F:2147
cp_fm_types::cp_fm_read_unformatted
subroutine, public cp_fm_read_unformatted(fm, unit)
...
Definition: cp_fm_types.F:2379
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition: cp_fm_types.F:1473
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::debug_print_level
integer, parameter, public debug_print_level
Definition: cp_output_handling.F:73
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_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition: cp_output_handling.F:704
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::do_admm_purify_none
integer, parameter, public do_admm_purify_none
Definition: input_constants.F:848
input_constants::do_loc_none
integer, parameter, public do_loc_none
Definition: input_constants.F:452
input_constants::op_loc_berry
integer, parameter, public op_loc_berry
Definition: input_constants.F:467
input_constants::xas_not_excited
integer, parameter, public xas_not_excited
Definition: input_constants.F:671
input_constants::tddfpt_singlet
integer, parameter, public tddfpt_singlet
Definition: input_constants.F:575
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::tddfpt_triplet
integer, parameter, public tddfpt_triplet
Definition: input_constants.F:575
input_constants::xas_2s_type
integer, parameter, public xas_2s_type
Definition: input_constants.F:671
input_constants::tddfpt_spin_flip
integer, parameter, public tddfpt_spin_flip
Definition: input_constants.F:575
input_constants::xas_tdp_by_kind
integer, parameter, public xas_tdp_by_kind
Definition: input_constants.F:688
input_constants::do_admm_purify_cauchy_subspace
integer, parameter, public do_admm_purify_cauchy_subspace
Definition: input_constants.F:848
input_constants::state_loc_list
integer, parameter, public state_loc_list
Definition: input_constants.F:467
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition: input_constants.F:801
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition: input_constants.F:801
input_constants::xas_tdp_by_index
integer, parameter, public xas_tdp_by_index
Definition: input_constants.F:688
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition: input_constants.F:801
input_constants::do_admm_purify_mo_diag
integer, parameter, public do_admm_purify_mo_diag
Definition: input_constants.F:848
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition: input_constants.F:801
input_constants::tddfpt_spin_cons
integer, parameter, public tddfpt_spin_cons
Definition: input_constants.F:575
input_cp2k_loc
Definition: input_cp2k_loc.F:7
input_cp2k_loc::create_localize_section
subroutine, public create_localize_section(section)
parameters fo the localization of wavefunctions
Definition: input_cp2k_loc.F:53
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_create
recursive subroutine, public section_vals_create(section_vals, section)
creates a object where to store the values of a section
Definition: input_section_types.F:601
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition: input_section_types.F:1216
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_release
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
Definition: input_section_types.F:231
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
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
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition: libint_wrapper.F:15
libint_wrapper::cp_libint_static_init
subroutine, public cp_libint_static_init()
Definition: libint_wrapper.F:1236
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
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::get_diag
pure real(kind=dp) function, dimension(min(size(a, 1), size(a, 2))), public get_diag(a)
Return the diagonal elements of matrix a as a vector.
Definition: mathlib.F:493
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
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition: parallel_rng_types.F:51
parallel_rng_types::uniform
integer, parameter, public uniform
Definition: parallel_rng_types.F:73
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::a_fine
real(kind=dp), parameter, public a_fine
Definition: physcon.F:119
physcon::evolt
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
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_interactions
Calculate the interaction radii for the operator matrix calculation.
Definition: qs_interactions.F:17
qs_interactions::init_interaction_radii_orb_basis
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Definition: qs_interactions.F:420
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::centers_spreads_berry
subroutine, public centers_spreads_berry(qs_loc_env, zij, nmoloc, cell, weights, ispin, print_loc_section, only_initial_out)
...
Definition: qs_loc_methods.F:576
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::localized_wfn_control_create
subroutine, public localized_wfn_control_create(localized_wfn_control)
create the localized_wfn_control_type
Definition: qs_loc_types.F:234
qs_loc_types::qs_loc_env_release
subroutine, public qs_loc_env_release(qs_loc_env)
...
Definition: qs_loc_types.F:192
qs_loc_types::get_qs_loc_env
subroutine, public get_qs_loc_env(qs_loc_env, cell, local_molecules, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set, para_env, particle_set, weights, dim_op)
...
Definition: qs_loc_types.F:317
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::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_mo_io
Definition and initialisation of the mo data type.
Definition: qs_mo_io.F:21
qs_mo_io::write_mo_set_low
subroutine, public write_mo_set_low(mo_array, qs_kind_set, particle_set, ires, rt_mos)
...
Definition: qs_mo_io.F:285
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::duplicate_mo_set
subroutine, public duplicate_mo_set(mo_set_new, mo_set_old)
allocate a new mo_set, and copy the old data
Definition: qs_mo_types.F:149
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition: qs_mo_types.F:206
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
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_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition: qs_mo_types.F:245
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_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
qs_scf_types::ot_method_nr
integer, parameter, public ot_method_nr
Definition: qs_scf_types.F:51
util
All kind of helpful little routines.
Definition: util.F:14
util::locate
pure integer function, public locate(array, x)
Purpose: Given an array array(1:n), and given a value x, a value x_index is returned which is the ind...
Definition: util.F:61
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333
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_tdp_atom
This module deals with all the integrals done on local atomic grids in xas_tdp. This is mostly used t...
Definition: xas_tdp_atom.F:15
xas_tdp_atom::init_xas_atom_env
subroutine, public init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env, ltddfpt)
Initializes a xas_atom_env type given the qs_enxas_atom_env, qs_envv.
Definition: xas_tdp_atom.F:159
xas_tdp_atom::integrate_soc_atoms
subroutine, public integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env, soc_atom_env)
Computes the SOC matrix elements with respect to the ORB basis set for each atomic kind and put them ...
Definition: xas_tdp_atom.F:3369
xas_tdp_atom::integrate_fxc_atoms
subroutine, public integrate_fxc_atoms(int_fxc, xas_atom_env, xas_tdp_control, qs_env)
Integrate the xc kernel as a function of r on the atomic grids for the RI_XAS basis.
Definition: xas_tdp_atom.F:1875
xas_tdp_correction
Second order perturbation correction to XAS_TDP spectra (i.e. shift)
Definition: xas_tdp_correction.F:13
xas_tdp_correction::gw2x_shift
subroutine, public gw2x_shift(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Computes the ionization potential using the GW2X method of Shigeta et. al. The result cam be used for...
Definition: xas_tdp_correction.F:103
xas_tdp_correction::get_soc_splitting
subroutine, public get_soc_splitting(soc_shifts, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
We try to compute the spin-orbit splitting via perturbation theory. We keep it \ cheap by only inculd...
Definition: xas_tdp_correction.F:1734
xas_tdp_integrals
3-center integrals machinery for the XAS_TDP method
Definition: xas_tdp_integrals.F:13
xas_tdp_integrals::compute_ri_coulomb2_int
subroutine, public compute_ri_coulomb2_int(ex_kind, xas_tdp_env, xas_tdp_control, qs_env)
Computes the two-center Coulomb integral needed for the RI in kernel calculation. Stores the integral...
Definition: xas_tdp_integrals.F:1400
xas_tdp_integrals::compute_ri_3c_coulomb
subroutine, public compute_ri_3c_coulomb(xas_tdp_env, qs_env)
Computes the RI Coulomb 3-center integrals (ab|c), where c is from the RI_XAS basis and centered on t...
Definition: xas_tdp_integrals.F:1146
xas_tdp_integrals::compute_ri_exchange2_int
subroutine, public compute_ri_exchange2_int(ex_kind, xas_tdp_env, xas_tdp_control, qs_env)
Computes the two-center Exchange integral needed for the RI in kernel calculation....
Definition: xas_tdp_integrals.F:1241
xas_tdp_integrals::compute_ri_3c_exchange
subroutine, public compute_ri_3c_exchange(ex_atoms, xas_tdp_env, xas_tdp_control, qs_env)
Computes the RI exchange 3-center integrals (ab|c), where c is from the RI_XAS basis and centered on ...
Definition: xas_tdp_integrals.F:1059
xas_tdp_methods
Methods for X-Ray absorption spectroscopy (XAS) using TDDFPT.
Definition: xas_tdp_methods.F:13
xas_tdp_methods::xas_tdp_init
subroutine, public xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
Overall control and environment types initialization.
Definition: xas_tdp_methods.F:556
xas_tdp_methods::xas_tdp
subroutine, public xas_tdp(qs_env)
Driver for XAS TDDFT calculations.
Definition: xas_tdp_methods.F:164
xas_tdp_types
Define XAS TDP control type and associated create, release, etc subroutines, as well as XAS TDP envir...
Definition: xas_tdp_types.F:13
xas_tdp_types::xas_tdp_env_create
subroutine, public xas_tdp_env_create(xas_tdp_env)
Creates a TDP XAS environment type.
Definition: xas_tdp_types.F:723
xas_tdp_types::xas_tdp_env_release
subroutine, public xas_tdp_env_release(xas_tdp_env)
Releases the TDP XAS environment type.
Definition: xas_tdp_types.F:768
xas_tdp_types::donor_state_create
subroutine, public donor_state_create(donor_state)
Creates a donor_state.
Definition: xas_tdp_types.F:909
xas_tdp_types::xas_tdp_control_release
subroutine, public xas_tdp_control_release(xas_tdp_control)
Releases the xas_tdp_control_type.
Definition: xas_tdp_types.F:445
xas_tdp_types::xas_atom_env_create
subroutine, public xas_atom_env_create(xas_atom_env)
Creates a xas_atom_env type.
Definition: xas_tdp_types.F:1201
xas_tdp_types::set_xas_tdp_env
subroutine, public set_xas_tdp_env(xas_tdp_env, nex_atoms, nex_kinds)
Sets values of selected variables within the TDP XAS environment type.
Definition: xas_tdp_types.F:893
xas_tdp_types::free_ds_memory
subroutine, public free_ds_memory(donor_state)
Deallocate a donor_state's heavy attributes.
Definition: xas_tdp_types.F:1108
xas_tdp_types::set_donor_state
subroutine, public set_donor_state(donor_state, at_index, at_symbol, kind_index, state_type)
sets specified values of the donor state type
Definition: xas_tdp_types.F:950
xas_tdp_types::xas_tdp_control_create
subroutine, public xas_tdp_control_create(xas_tdp_control)
Creates and initializes the xas_tdp_control_type.
Definition: xas_tdp_types.F:380
xas_tdp_types::free_exat_memory
subroutine, public free_exat_memory(xas_tdp_env, atom, end_of_batch)
Releases the memory heavy attribute of xas_tdp_env that are specific to the current excited atom.
Definition: xas_tdp_types.F:1365
xas_tdp_types::xas_atom_env_release
subroutine, public xas_atom_env_release(xas_atom_env)
Releases the xas_atom_env type.
Definition: xas_tdp_types.F:1230
xas_tdp_types::read_xas_tdp_control
subroutine, public read_xas_tdp_control(xas_tdp_control, xas_tdp_section)
Reads the inputs and stores in xas_tdp_control_type.
Definition: xas_tdp_types.F:479
xas_tdp_utils
Utilities for X-ray absorption spectroscopy using TDDFPT.
Definition: xas_tdp_utils.F:13
xas_tdp_utils::include_rcs_soc
subroutine, public include_rcs_soc(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Includes the SOC effects on the precomputed restricted closed-shell singlet and triplet excitations....
Definition: xas_tdp_utils.F:1939
xas_tdp_utils::setup_xas_tdp_prob
subroutine, public setup_xas_tdp_prob(donor_state, qs_env, xas_tdp_env, xas_tdp_control)
Builds the matrix that defines the XAS TDDFPT generalized eigenvalue problem to be solved for excitat...
Definition: xas_tdp_utils.F:168
xas_tdp_utils::solve_xas_tdp_prob
subroutine, public solve_xas_tdp_prob(donor_state, xas_tdp_control, xas_tdp_env, qs_env, ex_type)
Solves the XAS TDP generalized eigenvalue problem omega*C = matrix_tdp*C using standard full diagonal...
Definition: xas_tdp_utils.F:458
xas_tdp_utils::include_os_soc
subroutine, public include_os_soc(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Includes the SOC effects on the precomputed spin-conserving and spin-flip excitations from an open-sh...
Definition: xas_tdp_utils.F:1502
xc_write_output
Writes information on XC functionals to output.
Definition: xc_write_output.F:11
xc_write_output::xc_write
subroutine, public xc_write(iounit, xc_section, lsd)
...
Definition: xc_write_output.F:42