(git:374b731)
Loading...
Searching...
No Matches
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
13MODULE 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
155CONTAINS
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
3533END MODULE xas_tdp_methods
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
cell_types::pbc
Definition cell_types.F:92
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:66
util::sort_unique
Definition util.F:36
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
mathlib::diag
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition mathlib.F:1507
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
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
basis_set_types::sto_basis_set_type
Definition basis_set_types.F:100
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_loc_types::localized_wfn_control_type
A type that holds controlling information for the calculation of the spread of wfn and the optimizati...
Definition qs_loc_types.F:121
qs_loc_types::qs_loc_env_type
contains all the info needed by quickstep to calculate the spread of a selected set of orbitals and i...
Definition qs_loc_types.F:80
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
xas_tdp_types::donor_state_type
Type containing informations about a single donor state.
Definition xas_tdp_types.F:269
xas_tdp_types::xas_atom_env_type
a environment type that contains all the info needed for XAS_TDP atomic grid calculations
Definition xas_tdp_types.F:343
xas_tdp_types::xas_tdp_control_type
Type containing control information for TDP XAS calculations.
Definition xas_tdp_types.F:102
xas_tdp_types::xas_tdp_env_type
Type containing informations such as inputs and results for TDP XAS calculations.
Definition xas_tdp_types.F:194