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