(git:f56c6e3)
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-2025 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 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
36 copy_fm_to_dbcsr,&
37 cp_dbcsr_sm_fm_multiply
38 USE cp_files, ONLY: close_file,&
39 open_file
40 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
41 USE cp_fm_diag, ONLY: cp_fm_geeig,&
42 cp_fm_power,&
43 cp_fm_syevd
44 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
45 cp_fm_struct_release,&
46 cp_fm_struct_type
47 USE cp_fm_types, ONLY: &
48 cp_fm_copy_general, cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, &
49 cp_fm_read_unformatted, cp_fm_release, cp_fm_set_all, cp_fm_to_fm, cp_fm_to_fm_submat, &
50 cp_fm_type, cp_fm_write_unformatted
51 USE cp_log_handling, ONLY: cp_get_default_logger,&
52 cp_logger_get_default_io_unit,&
53 cp_logger_type,&
54 cp_to_string
55 USE cp_output_handling, ONLY: cp_p_file,&
56 cp_print_key_finished_output,&
57 cp_print_key_generate_filename,&
58 cp_print_key_should_output,&
59 cp_print_key_unit_nr,&
60 debug_print_level
61 USE input_constants, ONLY: &
62 do_admm_purify_cauchy_subspace, do_admm_purify_mo_diag, do_admm_purify_none, do_loc_none, &
63 do_potential_coulomb, do_potential_id, do_potential_short, do_potential_truncated, &
64 op_loc_berry, state_loc_list, tddfpt_singlet, tddfpt_spin_cons, tddfpt_spin_flip, &
65 tddfpt_triplet, xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, xas_dip_vel, &
66 xas_not_excited, xas_tdp_by_index, xas_tdp_by_kind
67 USE input_cp2k_loc, ONLY: create_localize_section
68 USE input_section_types, ONLY: section_release,&
69 section_type,&
70 section_vals_create,&
71 section_vals_get_subs_vals,&
72 section_vals_type,&
73 section_vals_val_get,&
74 section_vals_val_set
75 USE kinds, ONLY: default_path_length,&
76 default_string_length,&
77 dp
78 USE libint_wrapper, ONLY: cp_libint_static_init
79 USE machine, ONLY: m_flush
80 USE mathlib, ONLY: get_diag
81 USE memory_utilities, ONLY: reallocate
82 USE message_passing, ONLY: mp_comm_type,&
83 mp_para_env_type
84 USE parallel_gemm_api, ONLY: parallel_gemm
85 USE parallel_rng_types, ONLY: uniform,&
86 rng_stream_type
87 USE particle_methods, ONLY: get_particle_set
88 USE particle_types, ONLY: particle_type
89 USE periodic_table, ONLY: ptable
90 USE physcon, ONLY: a_fine,&
91 angstrom,&
92 evolt
93 USE qs_environment_types, ONLY: get_qs_env,&
94 qs_environment_type
95 USE qs_interactions, ONLY: init_interaction_radii_orb_basis
96 USE qs_kind_types, ONLY: get_qs_kind,&
97 qs_kind_type
98 USE qs_loc_main, ONLY: qs_loc_driver
99 USE qs_loc_methods, ONLY: centers_spreads_berry,&
100 qs_print_cubes
101 USE qs_loc_types, ONLY: get_qs_loc_env,&
102 localized_wfn_control_create,&
103 localized_wfn_control_type,&
104 qs_loc_env_create,&
105 qs_loc_env_release,&
106 qs_loc_env_type
107 USE qs_loc_utils, ONLY: qs_loc_control_init,&
108 qs_loc_env_init,&
109 set_loc_centers
110 USE qs_mo_io, ONLY: write_mo_set_low
111 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
112 USE qs_mo_types, ONLY: allocate_mo_set,&
113 deallocate_mo_set,&
114 duplicate_mo_set,&
115 get_mo_set,&
116 init_mo_set,&
117 mo_set_type
118 USE qs_operators_ao, ONLY: p_xyz_ao,&
119 rrc_xyz_ao
120 USE qs_pdos, ONLY: calculate_projected_dos
121 USE qs_rho_types, ONLY: qs_rho_get,&
122 qs_rho_type
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, rho_ao
598 TYPE(dbcsr_type) :: matrix_tmp
599 TYPE(dft_control_type), POINTER :: dft_control
600 TYPE(gto_basis_set_type), POINTER :: tmp_basis
601 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
602 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
603 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
604 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
605 TYPE(qs_rho_type), POINTER :: rho
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, rho, rho_ao, mo_evals, cell)
611 NULLIFY (mo_coeff, matrix_ks, admm_env, dummy_section)
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, rho=rho, matrix_s=matrix_s)
839 CALL qs_rho_get(rho, rho_ao=rho_ao)
840
841 ALLOCATE (xas_tdp_env%q_projector(nspins))
842 ALLOCATE (xas_tdp_env%q_projector(1)%matrix)
843 CALL dbcsr_create(xas_tdp_env%q_projector(1)%matrix, name="Q PROJECTOR ALPHA", &
844 template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
845 IF (do_os) THEN
846 ALLOCATE (xas_tdp_env%q_projector(2)%matrix)
847 CALL dbcsr_create(xas_tdp_env%q_projector(2)%matrix, name="Q PROJECTOR BETA", &
848 template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
849 END IF
850
851! In the case of spin-restricted calculations, rho_ao includes double occupency => 0.5 prefactor
852 fact = -0.5_dp; IF (do_os) fact = -1.0_dp
853 CALL dbcsr_multiply('N', 'N', fact, matrix_s(1)%matrix, rho_ao(1)%matrix, 0.0_dp, &
854 xas_tdp_env%q_projector(1)%matrix, filter_eps=xas_tdp_control%eps_filter)
855 CALL dbcsr_add_on_diag(xas_tdp_env%q_projector(1)%matrix, 1.0_dp)
856 CALL dbcsr_finalize(xas_tdp_env%q_projector(1)%matrix)
857
858 IF (do_os) THEN
859 CALL dbcsr_multiply('N', 'N', fact, matrix_s(1)%matrix, rho_ao(2)%matrix, 0.0_dp, &
860 xas_tdp_env%q_projector(2)%matrix, filter_eps=xas_tdp_control%eps_filter)
861 CALL dbcsr_add_on_diag(xas_tdp_env%q_projector(2)%matrix, 1.0_dp)
862 CALL dbcsr_finalize(xas_tdp_env%q_projector(2)%matrix)
863 END IF
864
865! Create the structure for the dipole in the AO basis
866 ALLOCATE (xas_tdp_env%dipmat(3))
867 DO i = 1, 3
868 ALLOCATE (xas_tdp_env%dipmat(i)%matrix)
869 CALL dbcsr_copy(matrix_tmp, matrix_s(1)%matrix, name="XAS TDP dipole matrix")
870 IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
871 CALL dbcsr_create(xas_tdp_env%dipmat(i)%matrix, template=matrix_s(1)%matrix, &
872 matrix_type=dbcsr_type_antisymmetric)
873 CALL dbcsr_complete_redistribute(matrix_tmp, xas_tdp_env%dipmat(i)%matrix)
874 ELSE
875 CALL dbcsr_create(xas_tdp_env%dipmat(i)%matrix, template=matrix_s(1)%matrix, &
876 matrix_type=dbcsr_type_symmetric)
877 CALL dbcsr_copy(xas_tdp_env%dipmat(i)%matrix, matrix_tmp)
878 END IF
879 CALL dbcsr_set(xas_tdp_env%dipmat(i)%matrix, 0.0_dp)
880 CALL dbcsr_release(matrix_tmp)
881 END DO
882
883! Create the structure for the electric quadrupole in the AO basis, if required
884 IF (xas_tdp_control%do_quad) THEN
885 ALLOCATE (xas_tdp_env%quadmat(6))
886 DO i = 1, 6
887 ALLOCATE (xas_tdp_env%quadmat(i)%matrix)
888 CALL dbcsr_copy(xas_tdp_env%quadmat(i)%matrix, matrix_s(1)%matrix, name="XAS TDP quadrupole matrix")
889 CALL dbcsr_set(xas_tdp_env%quadmat(i)%matrix, 0.0_dp)
890 END DO
891 END IF
892
893! Precompute it in the velocity representation, if so chosen
894 IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
895 !enforce minimum image to avoid any PBCs related issues. Ok because very localized densities
896 CALL p_xyz_ao(xas_tdp_env%dipmat, qs_env, minimum_image=.true.)
897 END IF
898
899! Allocate SOC in AO basis matrices
900 IF (xas_tdp_control%do_soc .OR. xas_tdp_control%do_gw2x) THEN
901 ALLOCATE (xas_tdp_env%orb_soc(3))
902 DO i = 1, 3
903 ALLOCATE (xas_tdp_env%orb_soc(i)%matrix)
904 END DO
905 END IF
906
907! Check that everything is allowed
908 CALL safety_check(xas_tdp_control, qs_env)
909
910! Initialize libint for the 3-center integrals
911 CALL cp_libint_static_init()
912
913! Compute LUMOs as guess for OT solver and/or for GW2X correction
914 IF (xas_tdp_control%do_ot .OR. xas_tdp_control%do_gw2x) THEN
915 CALL make_lumo_guess(xas_tdp_env, xas_tdp_control, qs_env)
916 END IF
917
918 END SUBROUTINE xas_tdp_init
919
920! **************************************************************************************************
921!> \brief splits the excited atoms of a kind into batches for RI 3c integrals load balance
922!> \param ex_atoms_of_kind the excited atoms for the current kind, randomly shuffled
923!> \param nbatch number of batches to loop over
924!> \param batch_size standard size of a batch
925!> \param atoms_of_kind number of atoms for the current kind (excited or not)
926!> \param xas_tdp_env ...
927! **************************************************************************************************
928 SUBROUTINE get_ri_3c_batches(ex_atoms_of_kind, nbatch, batch_size, atoms_of_kind, xas_tdp_env)
929
930 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: ex_atoms_of_kind
931 INTEGER, INTENT(OUT) :: nbatch
932 INTEGER, INTENT(IN) :: batch_size
933 INTEGER, DIMENSION(:), INTENT(IN) :: atoms_of_kind
934 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
935
936 INTEGER :: iat, iatom, nex_atom
937 TYPE(rng_stream_type), ALLOCATABLE :: rng_stream
938
939 !Get the atoms from atoms_of_kind that are excited
940 nex_atom = 0
941 DO iat = 1, SIZE(atoms_of_kind)
942 iatom = atoms_of_kind(iat)
943 IF (.NOT. any(xas_tdp_env%ex_atom_indices == iatom)) cycle
944 nex_atom = nex_atom + 1
945 END DO
946
947 ALLOCATE (ex_atoms_of_kind(nex_atom))
948 nex_atom = 0
949 DO iat = 1, SIZE(atoms_of_kind)
950 iatom = atoms_of_kind(iat)
951 IF (.NOT. any(xas_tdp_env%ex_atom_indices == iatom)) cycle
952 nex_atom = nex_atom + 1
953 ex_atoms_of_kind(nex_atom) = iatom
954 END DO
955
956 !We shuffle those atoms to spread them
957 rng_stream = rng_stream_type(name="uniform_rng", distribution_type=uniform)
958 CALL rng_stream%shuffle(ex_atoms_of_kind(1:nex_atom))
959
960 nbatch = nex_atom/batch_size
961 IF (nbatch*batch_size /= nex_atom) nbatch = nbatch + 1
962
963 END SUBROUTINE get_ri_3c_batches
964
965! **************************************************************************************************
966!> \brief Checks for forbidden keywords combinations
967!> \param xas_tdp_control ...
968!> \param qs_env ...
969! **************************************************************************************************
970 SUBROUTINE safety_check(xas_tdp_control, qs_env)
971
972 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
973 TYPE(qs_environment_type), POINTER :: qs_env
974
975 TYPE(dft_control_type), POINTER :: dft_control
976
977 !PB only available without exact exchange
978 IF (xas_tdp_control%is_periodic .AND. xas_tdp_control%do_hfx &
979 .AND. xas_tdp_control%x_potential%potential_type == do_potential_coulomb) THEN
980 cpabort("XAS TDP with Coulomb operator for exact exchange only supports non-periodic BCs")
981 END IF
982
983 !open-shell/closed-shell tests
984 IF (xas_tdp_control%do_roks .OR. xas_tdp_control%do_uks) THEN
985
986 IF (.NOT. (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_spin_flip)) THEN
987 cpabort("Need spin-conserving and/or spin-flip excitations for open-shell systems")
988 END IF
989
990 IF (xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet) THEN
991 cpabort("Singlet/triplet excitations only for restricted closed-shell systems")
992 END IF
993
994 IF (xas_tdp_control%do_soc .AND. .NOT. &
995 (xas_tdp_control%do_spin_flip .AND. xas_tdp_control%do_spin_cons)) THEN
996
997 cpabort("Both spin-conserving and spin-flip excitations are required for SOC")
998 END IF
999 ELSE
1000
1001 IF (.NOT. (xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet)) THEN
1002 cpabort("Need singlet and/or triplet excitations for closed-shell systems")
1003 END IF
1004
1005 IF (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_spin_flip) THEN
1006 cpabort("Spin-conserving/spin-flip excitations only for open-shell systems")
1007 END IF
1008
1009 IF (xas_tdp_control%do_soc .AND. .NOT. &
1010 (xas_tdp_control%do_singlet .AND. xas_tdp_control%do_triplet)) THEN
1011
1012 cpabort("Both singlet and triplet excitations are needed for SOC")
1013 END IF
1014 END IF
1015
1016 !Warn against using E_RANGE with SOC
1017 IF (xas_tdp_control%do_soc .AND. xas_tdp_control%e_range > 0.0_dp) THEN
1018 cpwarn("Using E_RANGE and SOC together may lead to crashes, use N_EXCITED for safety.")
1019 END IF
1020
1021 !TDA, full-TDDFT and diagonalization
1022 IF (.NOT. xas_tdp_control%tamm_dancoff) THEN
1023
1024 IF (xas_tdp_control%do_spin_flip) THEN
1025 cpabort("Spin-flip kernel only implemented for Tamm-Dancoff approximation")
1026 END IF
1027
1028 IF (xas_tdp_control%do_ot) THEN
1029 cpabort("OT diagonalization only available within the Tamm-Dancoff approximation")
1030 END IF
1031 END IF
1032
1033 !GW2X, need hfx kernel and LOCALIZE
1034 IF (xas_tdp_control%do_gw2x) THEN
1035 IF (.NOT. xas_tdp_control%do_hfx) THEN
1036 cpabort("GW2x requires the definition of the EXACT_EXCHANGE kernel")
1037 END IF
1038 IF (.NOT. xas_tdp_control%do_loc) THEN
1039 cpabort("GW2X requires the LOCALIZE keyword in DONOR_STATES")
1040 END IF
1041 END IF
1042
1043 !Only allow ADMM schemes that correct for eigenvalues
1044 CALL get_qs_env(qs_env, dft_control=dft_control)
1045 IF (dft_control%do_admm) THEN
1046 IF ((.NOT. qs_env%admm_env%purification_method == do_admm_purify_none) .AND. &
1047 (.NOT. qs_env%admm_env%purification_method == do_admm_purify_cauchy_subspace) .AND. &
1048 (.NOT. qs_env%admm_env%purification_method == do_admm_purify_mo_diag)) THEN
1049
1050 cpabort("XAS_TDP only compatible with ADMM purification NONE, CAUCHY_SUBSPACE and MO_DIAG")
1051
1052 END IF
1053 END IF
1054
1055 END SUBROUTINE safety_check
1056
1057! **************************************************************************************************
1058!> \brief Prints some basic info about the chosen parameters
1059!> \param ou the output unis
1060!> \param xas_tdp_control ...
1061!> \param qs_env ...
1062! **************************************************************************************************
1063 SUBROUTINE print_info(ou, xas_tdp_control, qs_env)
1064
1065 INTEGER, INTENT(IN) :: ou
1066 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1067 TYPE(qs_environment_type), POINTER :: qs_env
1068
1069 INTEGER :: i
1070 REAL(dp) :: occ
1071 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1072 TYPE(dft_control_type), POINTER :: dft_control
1073 TYPE(section_vals_type), POINTER :: input, kernel_section
1074
1075 NULLIFY (input, kernel_section, dft_control, matrix_s)
1076
1077 CALL get_qs_env(qs_env, input=input, dft_control=dft_control, matrix_s=matrix_s)
1078
1079 !Overlap matrix sparsity
1080 occ = dbcsr_get_occupation(matrix_s(1)%matrix)
1081
1082 IF (ou <= 0) RETURN
1083
1084 !Reference calculation
1085 IF (xas_tdp_control%do_uks) THEN
1086 WRITE (unit=ou, fmt="(/,T3,A)") &
1087 "XAS_TDP| Reference calculation: Unrestricted Kohn-Sham"
1088 ELSE IF (xas_tdp_control%do_roks) THEN
1089 WRITE (unit=ou, fmt="(/,T3,A)") &
1090 "XAS_TDP| Reference calculation: Restricted Open-Shell Kohn-Sham"
1091 ELSE
1092 WRITE (unit=ou, fmt="(/,T3,A)") &
1093 "XAS_TDP| Reference calculation: Restricted Closed-Shell Kohn-Sham"
1094 END IF
1095
1096 !TDA
1097 IF (xas_tdp_control%tamm_dancoff) THEN
1098 WRITE (unit=ou, fmt="(T3,A)") &
1099 "XAS_TDP| Tamm-Dancoff Approximation (TDA): On"
1100 ELSE
1101 WRITE (unit=ou, fmt="(T3,A)") &
1102 "XAS_TDP| Tamm-Dancoff Approximation (TDA): Off"
1103 END IF
1104
1105 !Dipole form
1106 IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
1107 WRITE (unit=ou, fmt="(T3,A)") &
1108 "XAS_TDP| Transition Dipole Representation: VELOCITY"
1109 ELSE
1110 WRITE (unit=ou, fmt="(T3,A)") &
1111 "XAS_TDP| Transition Dipole Representation: LENGTH"
1112 END IF
1113
1114 !Quadrupole
1115 IF (xas_tdp_control%do_quad) THEN
1116 WRITE (unit=ou, fmt="(T3,A)") &
1117 "XAS_TDP| Transition Quadrupole: On"
1118 END IF
1119
1120 !EPS_PGF
1121 IF (xas_tdp_control%eps_pgf > 0.0_dp) THEN
1122 WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1123 "XAS_TDP| EPS_PGF_XAS: ", xas_tdp_control%eps_pgf
1124 ELSE
1125 WRITE (unit=ou, fmt="(T3,A,ES7.1,A)") &
1126 "XAS_TDP| EPS_PGF_XAS: ", dft_control%qs_control%eps_pgf_orb, " (= EPS_PGF_ORB)"
1127 END IF
1128
1129 !EPS_FILTER
1130 WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1131 "XAS_TDP| EPS_FILTER: ", xas_tdp_control%eps_filter
1132
1133 !Grid info
1134 IF (xas_tdp_control%do_xc) THEN
1135 WRITE (unit=ou, fmt="(T3,A)") &
1136 "XAS_TDP| Radial Grid(s) Info: Kind, na, nr"
1137 DO i = 1, SIZE(xas_tdp_control%grid_info, 1)
1138 WRITE (unit=ou, fmt="(T3,A,A6,A,A,A,A)") &
1139 " ", trim(xas_tdp_control%grid_info(i, 1)), ", ", &
1140 trim(xas_tdp_control%grid_info(i, 2)), ", ", trim(xas_tdp_control%grid_info(i, 3))
1141 END DO
1142 END IF
1143
1144 !No kernel
1145 IF (.NOT. xas_tdp_control%do_coulomb) THEN
1146 WRITE (unit=ou, fmt="(/,T3,A)") &
1147 "XAS_TDP| No kernel (standard DFT)"
1148 END IF
1149
1150 !XC kernel
1151 IF (xas_tdp_control%do_xc) THEN
1152
1153 WRITE (unit=ou, fmt="(/,T3,A,F5.2,A)") &
1154 "XAS_TDP| RI Region's Radius: ", xas_tdp_control%ri_radius*angstrom, " Ang"
1155
1156 WRITE (unit=ou, fmt="(T3,A,/)") &
1157 "XAS_TDP| XC Kernel Functional(s) used for the kernel:"
1158
1159 IF (qs_env%do_rixs) THEN
1160 kernel_section => section_vals_get_subs_vals(input, "PROPERTIES%RIXS%XAS_TDP%KERNEL")
1161 ELSE
1162 kernel_section => section_vals_get_subs_vals(input, "DFT%XAS_TDP%KERNEL")
1163 END IF
1164 CALL xc_write(ou, kernel_section, lsd=.true.)
1165 END IF
1166
1167 !HFX kernel
1168 IF (xas_tdp_control%do_hfx) THEN
1169 WRITE (unit=ou, fmt="(/,T3,A,/,/,T3,A,F5.3)") &
1170 "XAS_TDP| Exact Exchange Kernel: Yes ", &
1171 "EXACT_EXCHANGE| Scale: ", xas_tdp_control%sx
1172 IF (xas_tdp_control%x_potential%potential_type == do_potential_coulomb) THEN
1173 WRITE (unit=ou, fmt="(T3,A)") &
1174 "EXACT_EXCHANGE| Potential : Coulomb"
1175 ELSE IF (xas_tdp_control%x_potential%potential_type == do_potential_truncated) THEN
1176 WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,A)") &
1177 "EXACT_EXCHANGE| Potential: Truncated Coulomb", &
1178 "EXACT_EXCHANGE| Range: ", xas_tdp_control%x_potential%cutoff_radius*angstrom, ", (Ang)", &
1179 "EXACT_EXCHANGE| T_C_G_DATA: ", trim(xas_tdp_control%x_potential%filename)
1180 ELSE IF (xas_tdp_control%x_potential%potential_type == do_potential_short) THEN
1181 WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,F5.2,A,/,T3,A,ES7.1)") &
1182 "EXACT_EXCHANGE| Potential: Short Range", &
1183 "EXACT_EXCHANGE| Omega: ", xas_tdp_control%x_potential%omega, ", (1/a0)", &
1184 "EXACT_EXCHANGE| Effective Range: ", xas_tdp_control%x_potential%cutoff_radius*angstrom, ", (Ang)", &
1185 "EXACT_EXCHANGE| EPS_RANGE: ", xas_tdp_control%eps_range
1186 END IF
1187 IF (xas_tdp_control%eps_screen > 1.0e-16) THEN
1188 WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1189 "EXACT_EXCHANGE| EPS_SCREENING: ", xas_tdp_control%eps_screen
1190 END IF
1191
1192 !RI metric
1193 IF (xas_tdp_control%do_ri_metric) THEN
1194
1195 WRITE (unit=ou, fmt="(/,T3,A)") &
1196 "EXACT_EXCHANGE| Using a RI metric"
1197 IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_id) THEN
1198 WRITE (unit=ou, fmt="(T3,A)") &
1199 "EXACT_EXCHANGE RI_METRIC| Potential : Overlap"
1200 ELSE IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_truncated) THEN
1201 WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,A)") &
1202 "EXACT_EXCHANGE RI_METRIC| Potential: Truncated Coulomb", &
1203 "EXACT_EXCHANGE RI_METRIC| Range: ", xas_tdp_control%ri_m_potential%cutoff_radius &
1204 *angstrom, ", (Ang)", &
1205 "EXACT_EXCHANGE RI_METRIC| T_C_G_DATA: ", trim(xas_tdp_control%ri_m_potential%filename)
1206 ELSE IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_short) THEN
1207 WRITE (unit=ou, fmt="(T3,A,/,T3,A,F5.2,A,/,T3,A,F5.2,A,/,T3,A,ES7.1)") &
1208 "EXACT_EXCHANGE RI_METRIC| Potential: Short Range", &
1209 "EXACT_EXCHANGE RI_METRIC| Omega: ", xas_tdp_control%ri_m_potential%omega, ", (1/a0)", &
1210 "EXACT_EXCHANGE RI_METRIC| Effective Range: ", &
1211 xas_tdp_control%ri_m_potential%cutoff_radius*angstrom, ", (Ang)", &
1212 "EXACT_EXCHANGE RI_METRIC| EPS_RANGE: ", xas_tdp_control%eps_range
1213 END IF
1214 END IF
1215 ELSE
1216 WRITE (unit=ou, fmt="(/,T3,A,/)") &
1217 "XAS_TDP| Exact Exchange Kernel: No "
1218 END IF
1219
1220 !overlap mtrix occupation
1221 WRITE (unit=ou, fmt="(/,T3,A,F5.2)") &
1222 "XAS_TDP| Overlap matrix occupation: ", occ
1223
1224 !GW2X parameter
1225 IF (xas_tdp_control%do_gw2x) THEN
1226 WRITE (unit=ou, fmt="(T3,A,/)") &
1227 "XAS_TDP| GW2X correction enabled"
1228
1229 IF (xas_tdp_control%xps_only) THEN
1230 WRITE (unit=ou, fmt="(T3,A)") &
1231 "GW2X| Only computing ionizations potentials for XPS"
1232 END IF
1233
1234 IF (xas_tdp_control%pseudo_canonical) THEN
1235 WRITE (unit=ou, fmt="(T3,A)") &
1236 "GW2X| Using the pseudo-canonical scheme"
1237 ELSE
1238 WRITE (unit=ou, fmt="(T3,A)") &
1239 "GW2X| Using the GW2X* scheme"
1240 END IF
1241
1242 WRITE (unit=ou, fmt="(T3,A,ES7.1)") &
1243 "GW2X| EPS_GW2X: ", xas_tdp_control%gw2x_eps
1244
1245 WRITE (unit=ou, fmt="(T3,A,I5)") &
1246 "GW2X| contraction batch size: ", xas_tdp_control%batch_size
1247
1248 IF ((int(xas_tdp_control%c_os) /= 1) .OR. (int(xas_tdp_control%c_ss) /= 1)) THEN
1249 WRITE (unit=ou, fmt="(T3,A,F7.4,/,T3,A,F7.4)") &
1250 "GW2X| Same-spin scaling factor: ", xas_tdp_control%c_ss, &
1251 "GW2X| Opposite-spin scaling factor: ", xas_tdp_control%c_os
1252 END IF
1253
1254 END IF
1255
1256 END SUBROUTINE print_info
1257
1258! **************************************************************************************************
1259!> \brief Assosciate (possibly localized) lowest energy MOs to each excited atoms. The procedure
1260!> looks for MOs "centered" on the excited atoms by comparing distances. It
1261!> then fills the mos_of_ex_atoms arrays of the xas_tdp_env. Only the xas_tdp_control%n_search
1262!> lowest energy MOs are considered. Largely inspired by MI's implementation of XAS
1263!> It is assumed that the Berry phase is used to compute centers.
1264!> \param xas_tdp_env ...
1265!> \param xas_tdp_control ...
1266!> \param qs_env ...
1267!> \note Whether localization took place or not, the procedure is the same as centers are stored in
1268!> xas_tdp_env%qs_loc_env%localized_wfn_control%centers_set
1269!> Assumes that find_mo_centers has been run previously
1270! **************************************************************************************************
1271 SUBROUTINE assign_mos_to_ex_atoms(xas_tdp_env, xas_tdp_control, qs_env)
1272
1273 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1274 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1275 TYPE(qs_environment_type), POINTER :: qs_env
1276
1277 INTEGER :: at_index, iat, iat_memo, imo, ispin, &
1278 n_atoms, n_search, nex_atoms, nspins
1279 INTEGER, DIMENSION(3) :: perd_init
1280 INTEGER, DIMENSION(:, :, :), POINTER :: mos_of_ex_atoms
1281 REAL(dp) :: dist, dist_min
1282 REAL(dp), DIMENSION(3) :: at_pos, r_ac, wfn_center
1283 TYPE(cell_type), POINTER :: cell
1284 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1285 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1286
1287 NULLIFY (localized_wfn_control, mos_of_ex_atoms, cell, particle_set)
1288
1289! Initialization. mos_of_ex_atoms filled with -1, meaning no assigned state
1290 mos_of_ex_atoms => xas_tdp_env%mos_of_ex_atoms
1291 mos_of_ex_atoms(:, :, :) = -1
1292 n_search = xas_tdp_control%n_search
1293 nex_atoms = xas_tdp_env%nex_atoms
1294 localized_wfn_control => xas_tdp_env%qs_loc_env%localized_wfn_control
1295 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
1296 n_atoms = SIZE(particle_set)
1297 nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1298
1299! Temporarly impose periodic BCs because of Berry's phase operator used for localization
1300 perd_init = cell%perd
1301 cell%perd = 1
1302
1303! Loop over n_search lowest energy MOs and all atoms, for each spin
1304 DO ispin = 1, nspins
1305 DO imo = 1, n_search
1306! retrieve MO wave function center coordinates.
1307 wfn_center(1:3) = localized_wfn_control%centers_set(ispin)%array(1:3, imo)
1308 iat_memo = 0
1309
1310! a large enough value to avoid bad surprises
1311 dist_min = 10000.0_dp
1312 DO iat = 1, n_atoms
1313 at_pos = particle_set(iat)%r
1314 r_ac = pbc(at_pos, wfn_center, cell)
1315 dist = norm2(r_ac)
1316
1317! keep memory of which atom is the closest to the wave function center
1318 IF (dist < dist_min) THEN
1319 iat_memo = iat
1320 dist_min = dist
1321 END IF
1322 END DO
1323
1324! Verify that the closest atom is actually excited and assign the MO if so
1325 IF (any(xas_tdp_env%ex_atom_indices == iat_memo)) THEN
1326 at_index = locate(xas_tdp_env%ex_atom_indices, iat_memo)
1327 mos_of_ex_atoms(imo, at_index, ispin) = 1
1328 END IF
1329 END DO !imo
1330 END DO !ispin
1331
1332! Go back to initial BCs
1333 cell%perd = perd_init
1334
1335 END SUBROUTINE assign_mos_to_ex_atoms
1336
1337! **************************************************************************************************
1338!> \brief Re-initialize the qs_loc_env to the current MOs.
1339!> \param qs_loc_env the env to re-initialize
1340!> \param n_loc_states the number of states to include
1341!> \param do_uks in cas of spin unrestricted calculation, initialize for both spins
1342!> \param qs_env ...
1343!> \note Useful when one needs to make use of qs_loc features and it is either with canonical MOs
1344!> or the localized MOs have been modified. do_localize is overwritten.
1345!> Same loc range for both spins
1346! **************************************************************************************************
1347 SUBROUTINE reinit_qs_loc_env(qs_loc_env, n_loc_states, do_uks, qs_env)
1348
1349 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1350 INTEGER, INTENT(IN) :: n_loc_states
1351 LOGICAL, INTENT(IN) :: do_uks
1352 TYPE(qs_environment_type), POINTER :: qs_env
1353
1354 INTEGER :: i, nspins
1355 TYPE(localized_wfn_control_type), POINTER :: loc_wfn_control
1356
1357! First, release the old env
1358 CALL qs_loc_env_release(qs_loc_env)
1359
1360! Re-create it
1361 CALL qs_loc_env_create(qs_loc_env)
1362 CALL localized_wfn_control_create(qs_loc_env%localized_wfn_control)
1363 loc_wfn_control => qs_loc_env%localized_wfn_control
1364
1365! Initialize it
1366 loc_wfn_control%localization_method = do_loc_none
1367 loc_wfn_control%operator_type = op_loc_berry
1368 loc_wfn_control%nloc_states(:) = n_loc_states
1369 loc_wfn_control%eps_occ = 0.0_dp
1370 loc_wfn_control%lu_bound_states(1, :) = 1
1371 loc_wfn_control%lu_bound_states(2, :) = n_loc_states
1372 loc_wfn_control%set_of_states = state_loc_list
1373 loc_wfn_control%do_homo = .true.
1374 ALLOCATE (loc_wfn_control%loc_states(n_loc_states, 2))
1375 DO i = 1, n_loc_states
1376 loc_wfn_control%loc_states(i, :) = i
1377 END DO
1378
1379 nspins = 1; IF (do_uks) nspins = 2
1380 CALL set_loc_centers(loc_wfn_control, loc_wfn_control%nloc_states, nspins=nspins)
1381 ! need to set do_localize=.TRUE. because otherwise no routine works
1382 IF (do_uks) THEN
1383 CALL qs_loc_env_init(qs_loc_env, loc_wfn_control, qs_env, do_localize=.true.)
1384 ELSE
1385 CALL qs_loc_env_init(qs_loc_env, loc_wfn_control, qs_env, myspin=1, do_localize=.true.)
1386 END IF
1387
1388 END SUBROUTINE reinit_qs_loc_env
1389
1390! *************************************************************************************************
1391!> \brief Diagonalize the subset of previously localized MOs that are associated to each excited
1392!> atoms. Updates the MO coeffs accordingly.
1393!> \param xas_tdp_env ...
1394!> \param xas_tdp_control ...
1395!> \param qs_env ...
1396!> \note Needed because after localization, the MOs loose their identity (1s, 2s , 2p, etc)
1397! **************************************************************************************************
1398 SUBROUTINE diagonalize_assigned_mo_subset(xas_tdp_env, xas_tdp_control, qs_env)
1399
1400 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1401 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1402 TYPE(qs_environment_type), POINTER :: qs_env
1403
1404 INTEGER :: i, iat, ilmo, ispin, nao, nlmo, nspins
1405 REAL(dp), ALLOCATABLE, DIMENSION(:) :: evals
1406 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1407 TYPE(cp_fm_struct_type), POINTER :: ks_struct, lmo_struct
1408 TYPE(cp_fm_type) :: evecs, ks_fm, lmo_fm, work
1409 TYPE(cp_fm_type), POINTER :: mo_coeff
1410 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1411 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1412 TYPE(mp_para_env_type), POINTER :: para_env
1413
1414 NULLIFY (mos, mo_coeff, matrix_ks, para_env, blacs_env, lmo_struct, ks_struct)
1415
1416 ! Get what we need from qs_env
1417 CALL get_qs_env(qs_env, mos=mos, matrix_ks=matrix_ks, para_env=para_env, blacs_env=blacs_env)
1418
1419 nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1420
1421 ! Loop over the excited atoms and spin
1422 DO ispin = 1, nspins
1423 DO iat = 1, xas_tdp_env%nex_atoms
1424
1425 ! get the MOs
1426 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao)
1427
1428 ! count how many MOs are associated to this atom and create a fm/struct
1429 nlmo = count(xas_tdp_env%mos_of_ex_atoms(:, iat, ispin) == 1)
1430 CALL cp_fm_struct_create(lmo_struct, nrow_global=nao, ncol_global=nlmo, &
1431 para_env=para_env, context=blacs_env)
1432 CALL cp_fm_create(lmo_fm, lmo_struct)
1433 CALL cp_fm_create(work, lmo_struct)
1434
1435 CALL cp_fm_struct_create(ks_struct, nrow_global=nlmo, ncol_global=nlmo, &
1436 para_env=para_env, context=blacs_env)
1437 CALL cp_fm_create(ks_fm, ks_struct)
1438 CALL cp_fm_create(evecs, ks_struct)
1439
1440 ! Loop over the localized MOs associated to this atom
1441 i = 0
1442 DO ilmo = 1, xas_tdp_control%n_search
1443 IF (xas_tdp_env%mos_of_ex_atoms(ilmo, iat, ispin) == -1) cycle
1444
1445 i = i + 1
1446 ! put the coeff in our atom-restricted lmo_fm
1447 CALL cp_fm_to_fm_submat(mo_coeff, lmo_fm, nrow=nao, ncol=1, s_firstrow=1, &
1448 s_firstcol=ilmo, t_firstrow=1, t_firstcol=i)
1449
1450 END DO !ilmo
1451
1452 ! Computing the KS matrix in the subset of MOs
1453 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, lmo_fm, work, ncol=nlmo)
1454 CALL parallel_gemm('T', 'N', nlmo, nlmo, nao, 1.0_dp, lmo_fm, work, 0.0_dp, ks_fm)
1455
1456 ! Diagonalizing the KS matrix in the subset of MOs
1457 ALLOCATE (evals(nlmo))
1458 CALL cp_fm_syevd(ks_fm, evecs, evals)
1459 DEALLOCATE (evals)
1460
1461 ! Express the MOs in the basis that diagonalizes KS
1462 CALL parallel_gemm('N', 'N', nao, nlmo, nlmo, 1.0_dp, lmo_fm, evecs, 0.0_dp, work)
1463
1464 ! Replacing the new MOs back in the MO coeffs
1465 i = 0
1466 DO ilmo = 1, xas_tdp_control%n_search
1467 IF (xas_tdp_env%mos_of_ex_atoms(ilmo, iat, ispin) == -1) cycle
1468
1469 i = i + 1
1470 CALL cp_fm_to_fm_submat(work, mo_coeff, nrow=nao, ncol=1, s_firstrow=1, &
1471 s_firstcol=i, t_firstrow=1, t_firstcol=ilmo)
1472
1473 END DO
1474
1475 ! Excited atom clean-up
1476 CALL cp_fm_release(lmo_fm)
1477 CALL cp_fm_release(work)
1478 CALL cp_fm_struct_release(lmo_struct)
1479 CALL cp_fm_release(ks_fm)
1480 CALL cp_fm_release(evecs)
1481 CALL cp_fm_struct_release(ks_struct)
1482 END DO !iat
1483 END DO !ispin
1484
1485 END SUBROUTINE diagonalize_assigned_mo_subset
1486
1487! **************************************************************************************************
1488!> \brief Assign core MO(s) to a given donor_state, taking the type (1S, 2S, etc) into account.
1489!> The projection on a representative Slater-type orbital basis is used as a indicator.
1490!> It is assumed that MOs are already assigned to excited atoms based on their center
1491!> \param donor_state the donor_state to which a MO must be assigned
1492!> \param xas_tdp_env ...
1493!> \param xas_tdp_control ...
1494!> \param qs_env ...
1495! **************************************************************************************************
1496 SUBROUTINE assign_mos_to_donor_state(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
1497
1498 TYPE(donor_state_type), POINTER :: donor_state
1499 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1500 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1501 TYPE(qs_environment_type), POINTER :: qs_env
1502
1503 INTEGER :: at_index, i, iat, imo, ispin, l, my_mo, &
1504 n_search, n_states, nao, ndo_so, nj, &
1505 nsgf_kind, nsgf_sto, nspins, &
1506 output_unit, zval
1507 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: my_mos
1508 INTEGER, DIMENSION(2) :: next_best_overlap_ind
1509 INTEGER, DIMENSION(4, 7) :: ne
1510 INTEGER, DIMENSION(:), POINTER :: first_sgf, lq, nq
1511 INTEGER, DIMENSION(:, :, :), POINTER :: mos_of_ex_atoms
1512 LOGICAL :: unique
1513 REAL(dp) :: zeff
1514 REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, overlap, sto_overlap
1515 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: max_overlap
1516 REAL(dp), DIMENSION(2) :: next_best_overlap
1517 REAL(dp), DIMENSION(:), POINTER :: mo_evals, zeta
1518 REAL(dp), DIMENSION(:, :), POINTER :: overlap_matrix, tmp_coeff
1519 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1520 TYPE(cp_fm_struct_type), POINTER :: eval_mat_struct, gs_struct, matrix_struct
1521 TYPE(cp_fm_type) :: eval_mat, work_mat
1522 TYPE(cp_fm_type), POINTER :: gs_coeffs, mo_coeff
1523 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1524 TYPE(gto_basis_set_type), POINTER :: kind_basis_set, sto_to_gto_basis_set
1525 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1526 TYPE(mp_para_env_type), POINTER :: para_env
1527 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1528 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1529 TYPE(sto_basis_set_type), POINTER :: sto_basis_set
1530
1531 NULLIFY (sto_basis_set, sto_to_gto_basis_set, qs_kind_set, kind_basis_set, lq, nq, zeta)
1532 NULLIFY (overlap_matrix, mos, mo_coeff, mos_of_ex_atoms, tmp_coeff, first_sgf, particle_set)
1533 NULLIFY (mo_evals, matrix_ks, para_env, blacs_env)
1534 NULLIFY (eval_mat_struct, gs_struct, gs_coeffs)
1535
1536 output_unit = cp_logger_get_default_io_unit()
1537
1538 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, mos=mos, particle_set=particle_set, &
1539 matrix_ks=matrix_ks, para_env=para_env, blacs_env=blacs_env)
1540
1541 nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1542
1543! Construction of a STO that fits the type of orbital we look for
1544 ALLOCATE (zeta(1))
1545 ALLOCATE (lq(1))
1546 ALLOCATE (nq(1))
1547! Retrieving quantum numbers
1548 IF (donor_state%state_type == xas_1s_type) THEN
1549 nq(1) = 1
1550 lq(1) = 0
1551 n_states = 1
1552 ELSE IF (donor_state%state_type == xas_2s_type) THEN
1553 nq(1) = 2
1554 lq(1) = 0
1555 n_states = 1
1556 ELSE IF (donor_state%state_type == xas_2p_type) THEN
1557 nq(1) = 2
1558 lq(1) = 1
1559 n_states = 3
1560 ELSE
1561 cpabort("Procedure for required type not implemented")
1562 END IF
1563 ALLOCATE (my_mos(n_states, nspins))
1564 ALLOCATE (max_overlap(n_states, nspins))
1565
1566! Getting the atomic number
1567 CALL get_qs_kind(qs_kind_set(donor_state%kind_index), zeff=zeff)
1568 zval = int(zeff)
1569
1570! Electronic configuration (copied from MI's XAS)
1571 ne = 0
1572 DO l = 1, 4
1573 nj = 2*(l - 1) + 1
1574 DO i = l, 7
1575 ne(l, i) = ptable(zval)%e_conv(l - 1) - 2*nj*(i - l)
1576 ne(l, i) = max(ne(l, i), 0)
1577 ne(l, i) = min(ne(l, i), 2*nj)
1578 END DO
1579 END DO
1580
1581! computing zeta with the Slater sum rules
1582 zeta(1) = srules(zval, ne, nq(1), lq(1))
1583
1584! Allocating memory and initiate STO
1585 CALL allocate_sto_basis_set(sto_basis_set)
1586 CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, lq=lq, zet=zeta)
1587
1588! Some clean-up
1589 DEALLOCATE (nq, lq, zeta)
1590
1591! Expanding the STO into (normalized) GTOs for later calculations, use standard 3 gaussians
1592 CALL create_gto_from_sto_basis(sto_basis_set=sto_basis_set, &
1593 gto_basis_set=sto_to_gto_basis_set, &
1594 ngauss=3)
1595 sto_to_gto_basis_set%norm_type = 2
1596 CALL init_orb_basis_set(sto_to_gto_basis_set)
1597
1598! Retrieving the atomic kind related GTO in which MOs are expanded
1599 CALL get_qs_kind(qs_kind_set(donor_state%kind_index), basis_set=kind_basis_set)
1600
1601! Allocating and computing the overlap between the two basis (they share the same center)
1602 CALL get_gto_basis_set(gto_basis_set=kind_basis_set, nsgf=nsgf_kind)
1603 CALL get_gto_basis_set(gto_basis_set=sto_to_gto_basis_set, nsgf=nsgf_sto)
1604 ALLOCATE (overlap_matrix(nsgf_sto, nsgf_kind))
1605
1606! Making use of MI's subroutine
1607 CALL calc_stogto_overlap(sto_to_gto_basis_set, kind_basis_set, overlap_matrix)
1608
1609! Some clean-up
1610 CALL deallocate_sto_basis_set(sto_basis_set)
1611 CALL deallocate_gto_basis_set(sto_to_gto_basis_set)
1612
1613! Looping over the potential donor states to compute overlap with STO basis
1614 mos_of_ex_atoms => xas_tdp_env%mos_of_ex_atoms
1615 n_search = xas_tdp_control%n_search
1616 at_index = donor_state%at_index
1617 iat = locate(xas_tdp_env%ex_atom_indices, at_index)
1618 ALLOCATE (first_sgf(SIZE(particle_set))) !probably do not need that
1619 CALL get_particle_set(particle_set=particle_set, qs_kind_set=qs_kind_set, first_sgf=first_sgf)
1620 ALLOCATE (tmp_coeff(nsgf_kind, 1))
1621 ALLOCATE (sto_overlap(nsgf_kind))
1622 ALLOCATE (overlap(n_search))
1623
1624 next_best_overlap = 0.0_dp
1625 max_overlap = 0.0_dp
1626
1627 DO ispin = 1, nspins
1628
1629 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao)
1630 overlap = 0.0_dp
1631
1632 my_mo = 0
1633 DO imo = 1, n_search
1634 IF (mos_of_ex_atoms(imo, iat, ispin) > 0) THEN
1635
1636 sto_overlap = 0.0_dp
1637 tmp_coeff = 0.0_dp
1638
1639! Getting the relevant coefficients for the candidate state
1640 CALL cp_fm_get_submatrix(fm=mo_coeff, target_m=tmp_coeff, start_row=first_sgf(at_index), &
1641 start_col=imo, n_rows=nsgf_kind, n_cols=1, transpose=.false.)
1642
1643! Computing the product overlap_matrix*coeffs
1644 CALL dgemm('N', 'N', nsgf_sto, 1, nsgf_kind, 1.0_dp, overlap_matrix, nsgf_sto, &
1645 tmp_coeff, nsgf_kind, 0.0_dp, sto_overlap, nsgf_sto)
1646
1647! Each element of column vector sto_overlap is the overlap of a basis element of the
1648! generated STO basis with the kind specific orbital basis. Take the sum of the absolute
1649! values so that rotation (of the px, py, pz for example) does not hinder our search
1650 overlap(imo) = sum(abs(sto_overlap))
1651
1652 END IF
1653 END DO
1654
1655! Finding the best overlap(s)
1656 DO i = 1, n_states
1657 my_mo = maxloc(overlap, 1)
1658 my_mos(i, ispin) = my_mo
1659 max_overlap(i, ispin) = maxval(overlap, 1)
1660 overlap(my_mo) = 0.0_dp
1661 END DO
1662! Getting the next best overlap (for validation purposes)
1663 next_best_overlap(ispin) = maxval(overlap, 1)
1664 next_best_overlap_ind(ispin) = maxloc(overlap, 1)
1665
1666! Sort MO indices
1667 CALL sort_unique(my_mos(:, ispin), unique)
1668
1669 END DO !ispin
1670
1671! Some clean-up
1672 DEALLOCATE (overlap_matrix, tmp_coeff)
1673
1674! Dealing with the result
1675 IF (all(my_mos > 0) .AND. all(my_mos <= n_search)) THEN
1676! Assigning the MO indices to the donor_state
1677 ALLOCATE (donor_state%mo_indices(n_states, nspins))
1678 donor_state%mo_indices = my_mos
1679 donor_state%ndo_mo = n_states
1680
1681! Storing the MOs in the donor_state, as vectors column: first columns alpha spin, then beta
1682 CALL cp_fm_struct_create(gs_struct, nrow_global=nao, ncol_global=n_states*nspins, &
1683 para_env=para_env, context=blacs_env)
1684 ALLOCATE (donor_state%gs_coeffs)
1685 CALL cp_fm_create(donor_state%gs_coeffs, gs_struct)
1686
1687 IF (.NOT. ASSOCIATED(xas_tdp_env%mo_coeff)) THEN
1688 ALLOCATE (xas_tdp_env%mo_coeff(nspins))
1689 END IF
1690
1691 DO ispin = 1, nspins
1692 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1693 ! check if mo_coeff is copied before for another donor_state
1694 IF (.NOT. ASSOCIATED(xas_tdp_env%mo_coeff(ispin)%local_data)) THEN
1695 ! copy mo_coeff
1696 CALL cp_fm_get_info(matrix=mo_coeff, &
1697 matrix_struct=matrix_struct)
1698 CALL cp_fm_create(xas_tdp_env%mo_coeff(ispin), matrix_struct)
1699 CALL cp_fm_to_fm(mo_coeff, xas_tdp_env%mo_coeff(ispin))
1700 END IF
1701
1702 DO i = 1, n_states
1703 CALL cp_fm_to_fm_submat(msource=mo_coeff, mtarget=donor_state%gs_coeffs, nrow=nao, &
1704 ncol=1, s_firstrow=1, s_firstcol=my_mos(i, ispin), &
1705 t_firstrow=1, t_firstcol=(ispin - 1)*n_states + i)
1706 END DO
1707 END DO
1708 gs_coeffs => donor_state%gs_coeffs
1709
1710 !Keep the subset of the coeffs centered on the excited atom as global array (used a lot)
1711 ALLOCATE (donor_state%contract_coeffs(nsgf_kind, n_states*nspins))
1712 CALL cp_fm_get_submatrix(gs_coeffs, donor_state%contract_coeffs, start_row=first_sgf(at_index), &
1713 start_col=1, n_rows=nsgf_kind, n_cols=n_states*nspins)
1714
1715! Assigning corresponding energy eigenvalues and writing some info in standard input file
1716
1717 !standard eigenvalues as gotten from the KS diagonalization in the ground state
1718 IF (.NOT. xas_tdp_control%do_loc .AND. .NOT. xas_tdp_control%do_roks) THEN
1719 IF (output_unit > 0) THEN
1720 WRITE (unit=output_unit, fmt="(T5,A,/,T5,A,/,T5,A)") &
1721 "The following canonical MO(s) have been associated with the donor state(s)", &
1722 "based on the overlap with the components of a minimal STO basis: ", &
1723 " Spin MO index overlap(sum)"
1724 END IF
1725
1726 ALLOCATE (donor_state%energy_evals(n_states, nspins))
1727 donor_state%energy_evals = 0.0_dp
1728
1729! Canonical MO, no change in eigenvalues, only diagonal elements
1730 DO ispin = 1, nspins
1731 CALL get_mo_set(mos(ispin), eigenvalues=mo_evals)
1732 DO i = 1, n_states
1733 donor_state%energy_evals(i, ispin) = mo_evals(my_mos(i, ispin))
1734
1735 IF (output_unit > 0) THEN
1736 WRITE (unit=output_unit, fmt="(T46,I4,I11,F17.5)") &
1737 ispin, my_mos(i, ispin), max_overlap(i, ispin)
1738 END IF
1739 END DO
1740 END DO
1741
1742 !either localization of MOs or ROKS, in both cases the MO eigenvalues from the KS
1743 !digonalization mat have changed
1744 ELSE
1745 IF (output_unit > 0) THEN
1746 WRITE (unit=output_unit, fmt="(T5,A,/,T5,A,/,T5,A)") &
1747 "The following localized MO(s) have been associated with the donor state(s)", &
1748 "based on the overlap with the components of a minimal STO basis: ", &
1749 " Spin MO index overlap(sum)"
1750 END IF
1751
1752! Loop over the donor states and print
1753 DO ispin = 1, nspins
1754 DO i = 1, n_states
1755
1756! Print info
1757 IF (output_unit > 0) THEN
1758 WRITE (unit=output_unit, fmt="(T46,I4,I11,F17.5)") &
1759 ispin, my_mos(i, ispin), max_overlap(i, ispin)
1760 END IF
1761 END DO
1762 END DO
1763
1764! MO have been rotated or non-physical ROKS MO eigrenvalues:
1765! => need epsilon_ij = <psi_i|F|psi_j> = sum_{pq} c_{qi}c_{pj} F_{pq}
1766! Note: only have digonal elements by construction
1767 ndo_so = nspins*n_states
1768 CALL cp_fm_create(work_mat, gs_struct)
1769 CALL cp_fm_struct_create(eval_mat_struct, nrow_global=ndo_so, ncol_global=ndo_so, &
1770 para_env=para_env, context=blacs_env)
1771 CALL cp_fm_create(eval_mat, eval_mat_struct)
1772 ALLOCATE (diag(ndo_so))
1773
1774 IF (.NOT. xas_tdp_control%do_roks) THEN
1775
1776 ALLOCATE (donor_state%energy_evals(n_states, nspins))
1777 donor_state%energy_evals = 0.0_dp
1778
1779! Compute gs_coeff^T * matrix_ks * gs_coeff to get the epsilon_ij matrix
1780 DO ispin = 1, nspins
1781 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, gs_coeffs, work_mat, ncol=ndo_so)
1782 CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, work_mat, 0.0_dp, eval_mat)
1783
1784! Put the epsilon_ii into the donor_state. No off-diagonal element because of subset diag
1785 CALL cp_fm_get_diag(eval_mat, diag)
1786 donor_state%energy_evals(:, ispin) = diag((ispin - 1)*n_states + 1:ispin*n_states)
1787
1788 END DO
1789
1790 ELSE
1791 ! If ROKS, slightly different procedure => 2 KS matrices but one type of MOs
1792 ALLOCATE (donor_state%energy_evals(n_states, 2))
1793 donor_state%energy_evals = 0.0_dp
1794
1795! Compute gs_coeff^T * matrix_ks * gs_coeff to get the epsilon_ij matrix
1796 DO ispin = 1, 2
1797 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, gs_coeffs, work_mat, ncol=ndo_so)
1798 CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, work_mat, 0.0_dp, eval_mat)
1799
1800 CALL cp_fm_get_diag(eval_mat, diag)
1801 donor_state%energy_evals(:, ispin) = diag(:)
1802
1803 END DO
1804
1805 DEALLOCATE (diag)
1806 END IF
1807
1808! Clean-up
1809 CALL cp_fm_release(work_mat)
1810 CALL cp_fm_release(eval_mat)
1811 CALL cp_fm_struct_release(eval_mat_struct)
1812
1813 END IF ! do_localize and/or ROKS
1814
1815! Allocate and initialize GW2X corrected IPs as energy_evals
1816 ALLOCATE (donor_state%gw2x_evals(SIZE(donor_state%energy_evals, 1), SIZE(donor_state%energy_evals, 2)))
1817 donor_state%gw2x_evals(:, :) = donor_state%energy_evals(:, :)
1818
1819! Clean-up
1820 CALL cp_fm_struct_release(gs_struct)
1821 DEALLOCATE (first_sgf)
1822
1823 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T5,A)") " "
1824
1825 DO ispin = 1, nspins
1826 IF (output_unit > 0) THEN
1827 WRITE (unit=output_unit, fmt="(T5,A,I1,A,F7.5,A,I4)") &
1828 "The next best overlap for spin ", ispin, " is ", next_best_overlap(ispin), &
1829 " for MO with index ", next_best_overlap_ind(ispin)
1830 END IF
1831 END DO
1832 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T5,A)") " "
1833
1834 ELSE
1835 cpabort("A core donor state could not be assigned MO(s). Increasing NSEARCH might help.")
1836 END IF
1837
1838 END SUBROUTINE assign_mos_to_donor_state
1839
1840! **************************************************************************************************
1841!> \brief Compute the centers and spreads of (core) MOs using the Berry phase operator
1842!> \param xas_tdp_env ...
1843!> \param xas_tdp_control ...
1844!> \param qs_env ...
1845!> \note xas_tdp_env%qs_loc_env is used and modified. OK since no localization done after this
1846!> subroutine is used.
1847! **************************************************************************************************
1848 SUBROUTINE find_mo_centers(xas_tdp_env, xas_tdp_control, qs_env)
1849
1850 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1851 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1852 TYPE(qs_environment_type), POINTER :: qs_env
1853
1854 INTEGER :: dim_op, i, ispin, j, n_centers, nao, &
1855 nspins
1856 REAL(dp), DIMENSION(6) :: weights
1857 TYPE(cell_type), POINTER :: cell
1858 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1859 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1860 TYPE(cp_fm_type) :: opvec
1861 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij_fm_set
1862 TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
1863 TYPE(cp_fm_type), POINTER :: vectors
1864 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
1865 TYPE(mp_para_env_type), POINTER :: para_env
1866 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1867 TYPE(section_vals_type), POINTER :: print_loc_section, prog_run_info
1868
1869 NULLIFY (qs_loc_env, cell, print_loc_section, op_sm_set, moloc_coeff, vectors)
1870 NULLIFY (tmp_fm_struct, para_env, blacs_env, prog_run_info)
1871
1872! Initialization
1873 print_loc_section => xas_tdp_control%print_loc_subsection
1874 n_centers = xas_tdp_control%n_search
1875 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env, cell=cell)
1876
1877! Set print option to debug to keep clean output file
1878 prog_run_info => section_vals_get_subs_vals(print_loc_section, "PROGRAM_RUN_INFO")
1879 CALL section_vals_val_set(prog_run_info, keyword_name="_SECTION_PARAMETERS_", &
1880 i_val=debug_print_level)
1881
1882! Re-initialize the qs_loc_env to get the current MOs. Use force_loc because needed for centers
1883 CALL reinit_qs_loc_env(xas_tdp_env%qs_loc_env, n_centers, xas_tdp_control%do_uks, qs_env)
1884 qs_loc_env => xas_tdp_env%qs_loc_env
1885
1886! Get what we need from the qs_lovc_env
1887 CALL get_qs_loc_env(qs_loc_env=qs_loc_env, weights=weights, op_sm_set=op_sm_set, &
1888 moloc_coeff=moloc_coeff)
1889
1890! Prepare for zij
1891 vectors => moloc_coeff(1)
1892 CALL cp_fm_get_info(vectors, nrow_global=nao)
1893 CALL cp_fm_create(opvec, vectors%matrix_struct)
1894
1895 CALL cp_fm_struct_create(tmp_fm_struct, para_env=para_env, context=blacs_env, &
1896 ncol_global=n_centers, nrow_global=n_centers)
1897
1898 IF (cell%orthorhombic) THEN
1899 dim_op = 3
1900 ELSE
1901 dim_op = 6
1902 END IF
1903 ALLOCATE (zij_fm_set(2, dim_op))
1904 DO i = 1, dim_op
1905 DO j = 1, 2
1906 CALL cp_fm_create(zij_fm_set(j, i), tmp_fm_struct)
1907 END DO
1908 END DO
1909
1910 ! If spin-unrestricted, need to go spin by spin
1911 nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
1912
1913 DO ispin = 1, nspins
1914! zij computation, copied from qs_loc_methods:optimize_loc_berry
1915 vectors => moloc_coeff(ispin)
1916 DO i = 1, dim_op
1917 DO j = 1, 2
1918 CALL cp_fm_set_all(zij_fm_set(j, i), 0.0_dp)
1919 CALL cp_dbcsr_sm_fm_multiply(op_sm_set(j, i)%matrix, vectors, opvec, ncol=n_centers)
1920 CALL parallel_gemm("T", "N", n_centers, n_centers, nao, 1.0_dp, vectors, opvec, 0.0_dp, &
1921 zij_fm_set(j, i))
1922 END DO
1923 END DO
1924
1925! Compute centers (and spread)
1926 CALL centers_spreads_berry(qs_loc_env=qs_loc_env, zij=zij_fm_set, nmoloc=n_centers, &
1927 cell=cell, weights=weights, ispin=ispin, &
1928 print_loc_section=print_loc_section, only_initial_out=.true.)
1929 END DO !ispins
1930
1931! Clean-up
1932 CALL cp_fm_release(opvec)
1933 CALL cp_fm_struct_release(tmp_fm_struct)
1934 CALL cp_fm_release(zij_fm_set)
1935
1936! Make sure we leave with the correct do_loc value
1937 qs_loc_env%do_localize = xas_tdp_control%do_loc
1938
1939 END SUBROUTINE find_mo_centers
1940
1941! **************************************************************************************************
1942!> \brief Prints the MO to donor_state assocaition with overlap and Mulliken population analysis
1943!> \param xas_tdp_env ...
1944!> \param xas_tdp_control ...
1945!> \param qs_env ...
1946!> \note Called only in case of CHECK_ONLY run
1947! **************************************************************************************************
1948 SUBROUTINE print_checks(xas_tdp_env, xas_tdp_control, qs_env)
1949
1950 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1951 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1952 TYPE(qs_environment_type), POINTER :: qs_env
1953
1954 CHARACTER(LEN=default_string_length) :: kind_name
1955 INTEGER :: current_state_index, iat, iatom, ikind, &
1956 istate, output_unit, tmp_index
1957 INTEGER, DIMENSION(:), POINTER :: atoms_of_kind
1958 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1959 TYPE(donor_state_type), POINTER :: current_state
1960
1961 NULLIFY (atomic_kind_set, atoms_of_kind, current_state)
1962
1963 output_unit = cp_logger_get_default_io_unit()
1964
1965 IF (output_unit > 0) THEN
1966 WRITE (output_unit, "(/,T3,A,/,T3,A,/,T3,A)") &
1967 "# Check the donor states for their quality. They need to have a well defined type ", &
1968 " (1s, 2s, etc) which is indicated by the overlap. They also need to be localized, ", &
1969 " for which the Mulliken population analysis is one indicator (must be close to 1.0)"
1970 END IF
1971
1972! Loop over the donor states (as in the main xas_tdp loop)
1973 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
1974 current_state_index = 1
1975
1976 !loop over atomic kinds
1977 DO ikind = 1, SIZE(atomic_kind_set)
1978
1979 CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), name=kind_name, &
1980 atom_list=atoms_of_kind)
1981
1982 IF (.NOT. any(xas_tdp_env%ex_kind_indices == ikind)) cycle
1983
1984 !loop over atoms of kind
1985 DO iat = 1, SIZE(atoms_of_kind)
1986 iatom = atoms_of_kind(iat)
1987
1988 IF (.NOT. any(xas_tdp_env%ex_atom_indices == iatom)) cycle
1989 tmp_index = locate(xas_tdp_env%ex_atom_indices, iatom)
1990
1991 !loop over states of excited atom
1992 DO istate = 1, SIZE(xas_tdp_env%state_types, 1)
1993
1994 IF (xas_tdp_env%state_types(istate, tmp_index) == xas_not_excited) cycle
1995
1996 current_state => xas_tdp_env%donor_states(current_state_index)
1997 CALL set_donor_state(current_state, at_index=iatom, &
1998 at_symbol=kind_name, kind_index=ikind, &
1999 state_type=xas_tdp_env%state_types(istate, tmp_index))
2000
2001 IF (output_unit > 0) THEN
2002 WRITE (output_unit, "(/,T4,A,A2,A,I4,A,A,A)") &
2003 "-Donor state of type ", xas_tdp_env%state_type_char(current_state%state_type), &
2004 " for atom", current_state%at_index, " of kind ", trim(current_state%at_symbol), ":"
2005 END IF
2006
2007 !Assign the MOs and perform Mulliken
2008 CALL assign_mos_to_donor_state(current_state, xas_tdp_env, xas_tdp_control, qs_env)
2009 CALL perform_mulliken_on_donor_state(current_state, qs_env)
2010
2011 current_state_index = current_state_index + 1
2012 NULLIFY (current_state)
2013
2014 END DO !istate
2015 END DO !iat
2016 END DO !ikind
2017
2018 IF (output_unit > 0) THEN
2019 WRITE (output_unit, "(/,T5,A)") &
2020 "Use LOCALIZE and/or increase N_SEARCH for better results, if so required."
2021 END IF
2022
2023 END SUBROUTINE print_checks
2024
2025! **************************************************************************************************
2026!> \brief Computes the required multipole moment in the length representation for a given atom
2027!> \param iatom index of the given atom
2028!> \param xas_tdp_env ...
2029!> \param xas_tdp_control ...
2030!> \param qs_env ...
2031!> \note Assumes that wither dipole or quadrupole in length rep is required
2032! **************************************************************************************************
2033 SUBROUTINE compute_lenrep_multipole(iatom, xas_tdp_env, xas_tdp_control, qs_env)
2034
2035 INTEGER, INTENT(IN) :: iatom
2036 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2037 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2038 TYPE(qs_environment_type), POINTER :: qs_env
2039
2040 INTEGER :: i, order
2041 REAL(dp), DIMENSION(3) :: rc
2042 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: work
2043 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2044
2045 NULLIFY (work, particle_set)
2046
2047 CALL get_qs_env(qs_env, particle_set=particle_set)
2048 rc = particle_set(iatom)%r
2049
2050 ALLOCATE (work(9))
2051 IF (xas_tdp_control%dipole_form == xas_dip_len) THEN
2052 DO i = 1, 3
2053 CALL dbcsr_set(xas_tdp_env%dipmat(i)%matrix, 0.0_dp)
2054 work(i)%matrix => xas_tdp_env%dipmat(i)%matrix
2055 END DO
2056 order = 1
2057 END IF
2058 IF (xas_tdp_control%do_quad) THEN
2059 DO i = 1, 6
2060 CALL dbcsr_set(xas_tdp_env%quadmat(i)%matrix, 0.0_dp)
2061 work(3 + i)%matrix => xas_tdp_env%quadmat(i)%matrix
2062 END DO
2063 order = 2
2064 IF (xas_tdp_control%dipole_form == xas_dip_vel) order = -2
2065 END IF
2066
2067 !enforce minimum image to avoid PBCs related issues, ok because localized densities
2068 CALL rrc_xyz_ao(work, qs_env, rc, order=order, minimum_image=.true.)
2069 DEALLOCATE (work)
2070
2071 END SUBROUTINE compute_lenrep_multipole
2072
2073! **************************************************************************************************
2074!> \brief Computes the oscillator strength based on the dipole moment (velocity or length rep) for
2075!> all available excitation energies and store the results in the donor_state. There is no
2076!> triplet dipole in the spin-restricted ground state.
2077!> \param donor_state the donor state which is excited
2078!> \param xas_tdp_control ...
2079!> \param xas_tdp_env ...
2080!> \note The oscillator strength is a scalar: osc_str = 2/(3*omega)*(dipole_v)^2 in the velocity rep
2081!> or : osc_str = 2/3*omega*(dipole_r)^2 in the length representation
2082!> The formulae for the dipoles come from the trace of the dipole operator with the transition
2083!> densities, i.e. what we get from solving the xas_tdp problem. Same procedure with or wo TDA
2084! **************************************************************************************************
2085 SUBROUTINE compute_dipole_fosc(donor_state, xas_tdp_control, xas_tdp_env)
2086
2087 TYPE(donor_state_type), POINTER :: donor_state
2088 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2089 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2090
2091 CHARACTER(len=*), PARAMETER :: routinen = 'compute_dipole_fosc'
2092
2093 INTEGER :: handle, iosc, j, nao, ndo_mo, ndo_so, &
2094 ngs, nosc, nspins
2095 LOGICAL :: do_sc, do_sg
2096 REAL(dp) :: alpha_xyz, beta_xyz, osc_xyz, pref
2097 REAL(dp), ALLOCATABLE, DIMENSION(:) :: alpha_contr, beta_contr, tot_contr
2098 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dip_block
2099 REAL(dp), DIMENSION(:), POINTER :: lr_evals
2100 REAL(dp), DIMENSION(:, :), POINTER :: alpha_osc, beta_osc, osc_str
2101 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2102 TYPE(cp_fm_struct_type), POINTER :: col_struct, mat_struct
2103 TYPE(cp_fm_type) :: col_work, mat_work
2104 TYPE(cp_fm_type), POINTER :: lr_coeffs
2105 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat
2106 TYPE(mp_para_env_type), POINTER :: para_env
2107
2108 NULLIFY (dipmat, col_struct, mat_struct, para_env, blacs_env, lr_coeffs)
2109 NULLIFY (lr_evals, osc_str, alpha_osc, beta_osc)
2110
2111 CALL timeset(routinen, handle)
2112
2113! Initialization
2114 do_sc = xas_tdp_control%do_spin_cons
2115 do_sg = xas_tdp_control%do_singlet
2116 IF (do_sc) THEN
2117 nspins = 2
2118 lr_evals => donor_state%sc_evals
2119 lr_coeffs => donor_state%sc_coeffs
2120 ELSE IF (do_sg) THEN
2121 nspins = 1
2122 lr_evals => donor_state%sg_evals
2123 lr_coeffs => donor_state%sg_coeffs
2124 ELSE
2125 cpabort("Dipole oscilaltor strength only for singlets and spin-conserving excitations.")
2126 END IF
2127 ndo_mo = donor_state%ndo_mo
2128 ndo_so = ndo_mo*nspins
2129 ngs = ndo_so; IF (xas_tdp_control%do_roks) ngs = ndo_mo !in ROKS, same gs coeffs
2130 nosc = SIZE(lr_evals)
2131 ALLOCATE (donor_state%osc_str(nosc, 4), donor_state%alpha_osc(nosc, 4), donor_state%beta_osc(nosc, 4))
2132 osc_str => donor_state%osc_str
2133 alpha_osc => donor_state%alpha_osc
2134 beta_osc => donor_state%beta_osc
2135 osc_str = 0.0_dp
2136 alpha_osc = 0.0_dp
2137 beta_osc = 0.0_dp
2138 dipmat => xas_tdp_env%dipmat
2139
2140 ! do some work matrix initialization
2141 CALL cp_fm_get_info(donor_state%gs_coeffs, matrix_struct=col_struct, para_env=para_env, &
2142 context=blacs_env, nrow_global=nao)
2143 CALL cp_fm_struct_create(mat_struct, para_env=para_env, context=blacs_env, &
2144 nrow_global=ndo_so*nosc, ncol_global=ngs)
2145 CALL cp_fm_create(mat_work, mat_struct)
2146 CALL cp_fm_create(col_work, col_struct)
2147
2148 ALLOCATE (tot_contr(ndo_mo), dip_block(ndo_so, ngs), alpha_contr(ndo_mo), beta_contr(ndo_mo))
2149 pref = 2.0_dp; IF (do_sc) pref = 1.0_dp !because of singlet definition u = 1/sqrt(2)(c_a+c_b)
2150
2151! Looping over cartesian coord
2152 DO j = 1, 3
2153
2154 !Compute dip*gs_coeffs
2155 CALL cp_dbcsr_sm_fm_multiply(dipmat(j)%matrix, donor_state%gs_coeffs, col_work, ncol=ngs)
2156 !compute lr_coeffs*dip*gs_coeffs
2157 CALL parallel_gemm('T', 'N', ndo_so*nosc, ngs, nao, 1.0_dp, lr_coeffs, col_work, 0.0_dp, mat_work)
2158
2159 !Loop over the excited states
2160 DO iosc = 1, nosc
2161
2162 tot_contr = 0.0_dp
2163 CALL cp_fm_get_submatrix(fm=mat_work, target_m=dip_block, start_row=(iosc - 1)*ndo_so + 1, &
2164 start_col=1, n_rows=ndo_so, n_cols=ngs)
2165 IF (do_sg) THEN
2166 tot_contr(:) = get_diag(dip_block)
2167 ELSE IF (do_sc .AND. xas_tdp_control%do_uks) THEN
2168 alpha_contr(:) = get_diag(dip_block(1:ndo_mo, 1:ndo_mo))
2169 beta_contr(:) = get_diag(dip_block(ndo_mo + 1:ndo_so, ndo_mo + 1:ndo_so))
2170 tot_contr(:) = alpha_contr(:) + beta_contr(:)
2171 ELSE
2172 !roks
2173 alpha_contr(:) = get_diag(dip_block(1:ndo_mo, :))
2174 beta_contr(:) = get_diag(dip_block(ndo_mo + 1:ndo_so, :))
2175 tot_contr(:) = alpha_contr(:) + beta_contr(:)
2176 END IF
2177
2178 osc_xyz = sum(tot_contr)**2
2179 alpha_xyz = sum(alpha_contr)**2
2180 beta_xyz = sum(beta_contr)**2
2181
2182 alpha_osc(iosc, 4) = alpha_osc(iosc, 4) + alpha_xyz
2183 alpha_osc(iosc, j) = alpha_xyz
2184
2185 beta_osc(iosc, 4) = beta_osc(iosc, 4) + beta_xyz
2186 beta_osc(iosc, j) = beta_xyz
2187
2188 osc_str(iosc, 4) = osc_str(iosc, 4) + osc_xyz
2189 osc_str(iosc, j) = osc_xyz
2190
2191 END DO !iosc
2192 END DO !j
2193
2194 !compute the prefactor
2195 DO j = 1, 4
2196 IF (xas_tdp_control%dipole_form == xas_dip_len) THEN
2197 osc_str(:, j) = pref*2.0_dp/3.0_dp*lr_evals(:)*osc_str(:, j)
2198 alpha_osc(:, j) = pref*2.0_dp/3.0_dp*lr_evals(:)*alpha_osc(:, j)
2199 beta_osc(:, j) = pref*2.0_dp/3.0_dp*lr_evals(:)*beta_osc(:, j)
2200 ELSE
2201 osc_str(:, j) = pref*2.0_dp/3.0_dp/lr_evals(:)*osc_str(:, j)
2202 alpha_osc(:, j) = pref*2.0_dp/3.0_dp/lr_evals(:)*alpha_osc(:, j)
2203 beta_osc(:, j) = pref*2.0_dp/3.0_dp/lr_evals(:)*beta_osc(:, j)
2204 END IF
2205 END DO
2206
2207 !clean-up
2208 CALL cp_fm_release(mat_work)
2209 CALL cp_fm_release(col_work)
2210 CALL cp_fm_struct_release(mat_struct)
2211
2212 CALL timestop(handle)
2213
2214 END SUBROUTINE compute_dipole_fosc
2215
2216! **************************************************************************************************
2217!> \brief Computes the oscillator strength due to the electric quadrupole moment and store it in
2218!> the donor_state (for singlet or spin-conserving)
2219!> \param donor_state the donor state which is excited
2220!> \param xas_tdp_control ...
2221!> \param xas_tdp_env ...
2222!> \note Formula: 1/20*a_fine^2*omega^3 * sum_ab (sum_i r_ia*r_ib - 1/3*ri^2*delta_ab)
2223! **************************************************************************************************
2224 SUBROUTINE compute_quadrupole_fosc(donor_state, xas_tdp_control, xas_tdp_env)
2225
2226 TYPE(donor_state_type), POINTER :: donor_state
2227 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2228 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2229
2230 CHARACTER(len=*), PARAMETER :: routinen = 'compute_quadrupole_fosc'
2231
2232 INTEGER :: handle, iosc, j, nao, ndo_mo, ndo_so, &
2233 ngs, nosc, nspins
2234 LOGICAL :: do_sc, do_sg
2235 REAL(dp) :: pref
2236 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tot_contr, trace
2237 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: quad_block
2238 REAL(dp), DIMENSION(:), POINTER :: lr_evals, osc_str
2239 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2240 TYPE(cp_fm_struct_type), POINTER :: col_struct, mat_struct
2241 TYPE(cp_fm_type) :: col_work, mat_work
2242 TYPE(cp_fm_type), POINTER :: lr_coeffs
2243 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: quadmat
2244 TYPE(mp_para_env_type), POINTER :: para_env
2245
2246 NULLIFY (lr_evals, osc_str, lr_coeffs, col_struct, mat_struct, para_env)
2247 NULLIFY (blacs_env)
2248
2249 CALL timeset(routinen, handle)
2250
2251 ! Initialization
2252 do_sc = xas_tdp_control%do_spin_cons
2253 do_sg = xas_tdp_control%do_singlet
2254 IF (do_sc) THEN
2255 nspins = 2
2256 lr_evals => donor_state%sc_evals
2257 lr_coeffs => donor_state%sc_coeffs
2258 ELSE IF (do_sg) THEN
2259 nspins = 1
2260 lr_evals => donor_state%sg_evals
2261 lr_coeffs => donor_state%sg_coeffs
2262 ELSE
2263 cpabort("Quadrupole oscillator strengths only for singlet and spin-conserving excitations")
2264 END IF
2265 ndo_mo = donor_state%ndo_mo
2266 ndo_so = ndo_mo*nspins
2267 ngs = ndo_so; IF (xas_tdp_control%do_roks) ngs = ndo_mo !only alpha do_mo in ROKS
2268 nosc = SIZE(lr_evals)
2269 ALLOCATE (donor_state%quad_osc_str(nosc))
2270 osc_str => donor_state%quad_osc_str
2271 osc_str = 0.0_dp
2272 quadmat => xas_tdp_env%quadmat
2273
2274 !work matrices init
2275 CALL cp_fm_get_info(donor_state%gs_coeffs, matrix_struct=col_struct, para_env=para_env, &
2276 context=blacs_env, nrow_global=nao)
2277 CALL cp_fm_struct_create(mat_struct, para_env=para_env, context=blacs_env, &
2278 nrow_global=ndo_so*nosc, ncol_global=ngs)
2279 CALL cp_fm_create(mat_work, mat_struct)
2280 CALL cp_fm_create(col_work, col_struct)
2281
2282 ALLOCATE (quad_block(ndo_so, ngs), tot_contr(ndo_mo))
2283 pref = 2.0_dp; IF (do_sc) pref = 1.0_dp !because of singlet definition u = 1/sqrt(2)*...
2284 ALLOCATE (trace(nosc))
2285 trace = 0.0_dp
2286
2287 !Loop over the cartesioan coord :x2, xy, xz, y2, yz, z2
2288 DO j = 1, 6
2289
2290 !Compute quad*gs_coeffs
2291 CALL cp_dbcsr_sm_fm_multiply(quadmat(j)%matrix, donor_state%gs_coeffs, col_work, ncol=ngs)
2292 !compute lr_coeffs*quadmat*gs_coeffs
2293 CALL parallel_gemm('T', 'N', ndo_so*nosc, ngs, nao, 1.0_dp, lr_coeffs, col_work, 0.0_dp, mat_work)
2294
2295 !Loop over the excited states
2296 DO iosc = 1, nosc
2297
2298 tot_contr = 0.0_dp
2299 CALL cp_fm_get_submatrix(fm=mat_work, target_m=quad_block, start_row=(iosc - 1)*ndo_so + 1, &
2300 start_col=1, n_rows=ndo_so, n_cols=ngs)
2301
2302 IF (do_sg) THEN
2303 tot_contr(:) = get_diag(quad_block)
2304 ELSE IF (do_sc .AND. xas_tdp_control%do_uks) THEN
2305 tot_contr(:) = get_diag(quad_block(1:ndo_mo, 1:ndo_mo)) !alpha
2306 tot_contr(:) = tot_contr(:) + get_diag(quad_block(ndo_mo + 1:ndo_so, ndo_mo + 1:ndo_so)) !beta
2307 ELSE
2308 !roks
2309 tot_contr(:) = get_diag(quad_block(1:ndo_mo, :)) !alpha
2310 tot_contr(:) = tot_contr(:) + get_diag(quad_block(ndo_mo + 1:ndo_so, :)) !beta
2311 END IF
2312
2313 !if x2, y2, or z2 direction, need to update the trace (for later)
2314 IF (j == 1 .OR. j == 4 .OR. j == 6) THEN
2315 osc_str(iosc) = osc_str(iosc) + sum(tot_contr)**2
2316 trace(iosc) = trace(iosc) + sum(tot_contr)
2317
2318 !if xy, xz or yz, need to count twice the contribution (for yx, zx and zy)
2319 ELSE
2320 osc_str(iosc) = osc_str(iosc) + 2.0_dp*sum(tot_contr)**2
2321 END IF
2322
2323 END DO !iosc
2324 END DO !j
2325
2326 !compute the prefactor, and remove 1/3*trace^2
2327 osc_str(:) = pref*1._dp/20._dp*a_fine**2*lr_evals(:)**3*(osc_str(:) - 1._dp/3._dp*trace(:)**2)
2328
2329 !clean-up
2330 CALL cp_fm_release(mat_work)
2331 CALL cp_fm_release(col_work)
2332 CALL cp_fm_struct_release(mat_struct)
2333
2334 CALL timestop(handle)
2335
2336 END SUBROUTINE compute_quadrupole_fosc
2337
2338! **************************************************************************************************
2339!> \brief Writes the core MOs to excited atoms associations in the main output file
2340!> \param xas_tdp_env ...
2341!> \param xas_tdp_control ...
2342!> \param qs_env ...
2343!> \note Look at alpha spin MOs, as we are dealing with core states and alpha/beta MOs are the same
2344! **************************************************************************************************
2345 SUBROUTINE write_mos_to_ex_atoms_association(xas_tdp_env, xas_tdp_control, qs_env)
2346
2347 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2348 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2349 TYPE(qs_environment_type), POINTER :: qs_env
2350
2351 CHARACTER(LEN=default_string_length) :: kind_name
2352 INTEGER :: at_index, imo, ispin, nmo, nspins, &
2353 output_unit, tmp_index
2354 INTEGER, DIMENSION(3) :: perd_init
2355 INTEGER, DIMENSION(:), POINTER :: ex_atom_indices
2356 INTEGER, DIMENSION(:, :, :), POINTER :: mos_of_ex_atoms
2357 REAL(dp) :: dist, mo_spread
2358 REAL(dp), DIMENSION(3) :: at_pos, r_ac, wfn_center
2359 TYPE(cell_type), POINTER :: cell
2360 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2361
2362 NULLIFY (cell, particle_set, mos_of_ex_atoms, ex_atom_indices)
2363
2364 output_unit = cp_logger_get_default_io_unit()
2365
2366 IF (output_unit > 0) THEN
2367 WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,/,T3,A)") &
2368 " Associated Associated Distance to MO spread (Ang^2)", &
2369 "Spin MO index atom index atom kind MO center (Ang) -w_i ln(|z_ij|^2)", &
2370 "---------------------------------------------------------------------------------"
2371 END IF
2372
2373! Initialization
2374 nspins = 1; IF (xas_tdp_control%do_uks) nspins = 2
2375 mos_of_ex_atoms => xas_tdp_env%mos_of_ex_atoms
2376 ex_atom_indices => xas_tdp_env%ex_atom_indices
2377 nmo = xas_tdp_control%n_search
2378 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell)
2379
2380! because the use of Berry's phase operator implies PBCs
2381 perd_init = cell%perd
2382 cell%perd = 1
2383
2384! Retrieving all the info for each MO and spin
2385 DO imo = 1, nmo
2386 DO ispin = 1, nspins
2387
2388! each Mo is associated to at most one atom (only 1 in array of -1)
2389 IF (any(mos_of_ex_atoms(imo, :, ispin) == 1)) THEN
2390 tmp_index = maxloc(mos_of_ex_atoms(imo, :, ispin), 1)
2391 at_index = ex_atom_indices(tmp_index)
2392 kind_name = particle_set(at_index)%atomic_kind%name
2393
2394 at_pos = particle_set(at_index)%r
2395 wfn_center = xas_tdp_env%qs_loc_env%localized_wfn_control%centers_set(ispin)%array(1:3, imo)
2396 r_ac = pbc(at_pos, wfn_center, cell)
2397 dist = norm2(r_ac)
2398! convert distance from a.u. to Angstrom
2399 dist = dist*angstrom
2400
2401 mo_spread = xas_tdp_env%qs_loc_env%localized_wfn_control%centers_set(ispin)%array(4, imo)
2402 mo_spread = mo_spread*angstrom*angstrom
2403
2404 IF (output_unit > 0) THEN
2405 WRITE (unit=output_unit, fmt="(T3,I4,I10,I14,A14,ES19.3,ES20.3)") &
2406 ispin, imo, at_index, trim(kind_name), dist, mo_spread
2407 END IF
2408
2409 END IF
2410 END DO !ispin
2411 END DO !imo
2412
2413 IF (output_unit > 0) THEN
2414 WRITE (unit=output_unit, fmt="(T3,A,/)") &
2415 "---------------------------------------------------------------------------------"
2416 END IF
2417
2418! Go back to initial BCs
2419 cell%perd = perd_init
2420
2421 END SUBROUTINE write_mos_to_ex_atoms_association
2422
2423! **************************************************************************************************
2424!> \brief Performs Mulliken population analysis for the MO(s) of a donor_state_type so that user
2425!> can verify it is indeed a core state
2426!> \param donor_state ...
2427!> \param qs_env ...
2428!> \note This is a specific case of Mulliken analysis. In general one computes sum_i (SP)_ii, where
2429!> i labels the basis function centered on the atom of interest. For a specific MO with index
2430!> j, one need to compute sum_{ik} c_{ij} S_{ik} c_{kj}, k = 1,nao
2431! **************************************************************************************************
2432 SUBROUTINE perform_mulliken_on_donor_state(donor_state, qs_env)
2433 TYPE(donor_state_type), POINTER :: donor_state
2434 TYPE(qs_environment_type), POINTER :: qs_env
2435
2436 INTEGER :: at_index, i, ispin, nao, natom, ndo_mo, &
2437 ndo_so, nsgf, nspins, output_unit
2438 INTEGER, DIMENSION(:), POINTER :: first_sgf, last_sgf
2439 INTEGER, DIMENSION(:, :), POINTER :: mo_indices
2440 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: mul_pop, pop_mat
2441 REAL(dp), DIMENSION(:, :), POINTER :: work_array
2442 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2443 TYPE(cp_fm_struct_type), POINTER :: col_vect_struct
2444 TYPE(cp_fm_type) :: work_vect
2445 TYPE(cp_fm_type), POINTER :: gs_coeffs
2446 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2447 TYPE(mp_para_env_type), POINTER :: para_env
2448 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2449 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2450
2451 NULLIFY (mo_indices, qs_kind_set, particle_set, first_sgf, work_array)
2452 NULLIFY (matrix_s, para_env, blacs_env, col_vect_struct, last_sgf)
2453
2454! Initialization
2455 at_index = donor_state%at_index
2456 mo_indices => donor_state%mo_indices
2457 ndo_mo = donor_state%ndo_mo
2458 gs_coeffs => donor_state%gs_coeffs
2459 output_unit = cp_logger_get_default_io_unit()
2460 nspins = 1; IF (SIZE(mo_indices, 2) == 2) nspins = 2
2461 ndo_so = ndo_mo*nspins
2462 ALLOCATE (mul_pop(ndo_mo, nspins))
2463 mul_pop = 0.0_dp
2464
2465 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set, &
2466 para_env=para_env, blacs_env=blacs_env, matrix_s=matrix_s)
2467 CALL cp_fm_get_info(gs_coeffs, nrow_global=nao, matrix_struct=col_vect_struct)
2468
2469 natom = SIZE(particle_set, 1)
2470 ALLOCATE (first_sgf(natom))
2471 ALLOCATE (last_sgf(natom))
2472
2473 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
2474 nsgf = last_sgf(at_index) - first_sgf(at_index) + 1
2475
2476 CALL cp_fm_create(work_vect, col_vect_struct)
2477
2478! Take the product of S*coeffs
2479 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, gs_coeffs, work_vect, ncol=ndo_so)
2480
2481! Only consider the product coeffs^T * S * coeffs on the atom of interest
2482 ALLOCATE (work_array(nsgf, ndo_so))
2483 ALLOCATE (pop_mat(ndo_so, ndo_so))
2484
2485 CALL cp_fm_get_submatrix(fm=work_vect, target_m=work_array, start_row=first_sgf(at_index), &
2486 start_col=1, n_rows=nsgf, n_cols=ndo_so, transpose=.false.)
2487
2488 CALL dgemm('T', 'N', ndo_so, ndo_so, nsgf, 1.0_dp, donor_state%contract_coeffs, nsgf, &
2489 work_array, nsgf, 0.0_dp, pop_mat, ndo_so)
2490
2491! The Mullikan population for the MOs in on the diagonal.
2492 DO ispin = 1, nspins
2493 DO i = 1, ndo_mo
2494 mul_pop(i, ispin) = pop_mat((ispin - 1)*ndo_mo + i, (ispin - 1)*ndo_mo + i)
2495 END DO
2496 END DO
2497
2498! Printing in main output file
2499 IF (output_unit > 0) THEN
2500 WRITE (unit=output_unit, fmt="(T5,A,/,T5,A)") &
2501 "Mulliken population analysis retricted to the associated MO(s) yields: ", &
2502 " Spin MO index charge"
2503 DO ispin = 1, nspins
2504 DO i = 1, ndo_mo
2505 WRITE (unit=output_unit, fmt="(T51,I4,I10,F11.3)") &
2506 ispin, mo_indices(i, ispin), mul_pop(i, ispin)
2507 END DO
2508 END DO
2509 END IF
2510
2511! Clean-up
2512 DEALLOCATE (first_sgf, last_sgf, work_array)
2513 CALL cp_fm_release(work_vect)
2514
2515 END SUBROUTINE perform_mulliken_on_donor_state
2516
2517! **************************************************************************************************
2518!> \brief write the PDOS wrt the LR-orbitals for the current donor_state and/or the CUBES files
2519!> \param ex_type the excitation type: singlet, triplet, spin-conserving, etc
2520!> \param donor_state ...
2521!> \param xas_tdp_env ...
2522!> \param xas_tdp_section ...
2523!> \param qs_env ...
2524! **************************************************************************************************
2525 SUBROUTINE xas_tdp_post(ex_type, donor_state, xas_tdp_env, xas_tdp_section, qs_env)
2526
2527 INTEGER, INTENT(IN) :: ex_type
2528 TYPE(donor_state_type), POINTER :: donor_state
2529 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2530 TYPE(section_vals_type), POINTER :: xas_tdp_section
2531 TYPE(qs_environment_type), POINTER :: qs_env
2532
2533 CHARACTER(len=*), PARAMETER :: routinen = 'xas_tdp_post'
2534
2535 CHARACTER(len=default_string_length) :: domo, domon, excite, pos, xas_mittle
2536 INTEGER :: ex_state_idx, handle, ic, ido_mo, imo, irep, ispin, n_dependent, n_rep, nao, &
2537 ncubes, ndo_mo, ndo_so, nlumo, nmo, nspins, output_unit
2538 INTEGER, DIMENSION(:), POINTER :: bounds, list, state_list
2539 LOGICAL :: append_cube, do_cubes, do_pdos, &
2540 do_wfn_restart
2541 REAL(dp), DIMENSION(:), POINTER :: lr_evals
2542 REAL(dp), DIMENSION(:, :), POINTER :: centers
2543 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2544 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2545 TYPE(cp_fm_struct_type), POINTER :: fm_struct, mo_struct
2546 TYPE(cp_fm_type) :: mo_coeff, work_fm
2547 TYPE(cp_fm_type), POINTER :: lr_coeffs
2548 TYPE(cp_logger_type), POINTER :: logger
2549 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2550 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2551 TYPE(mo_set_type), POINTER :: mo_set
2552 TYPE(mp_para_env_type), POINTER :: para_env
2553 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2554 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2555 TYPE(section_vals_type), POINTER :: print_key
2556
2557 NULLIFY (atomic_kind_set, particle_set, qs_kind_set, mo_set, lr_evals, lr_coeffs)
2558 NULLIFY (mo_struct, para_env, blacs_env, fm_struct, matrix_s, print_key, logger)
2559 NULLIFY (bounds, state_list, list, mos)
2560
2561 !Tests on what to do
2562 logger => cp_get_default_logger()
2563 do_pdos = .false.; do_cubes = .false.; do_wfn_restart = .false.
2564
2565 IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
2566 "PRINT%PDOS"), cp_p_file)) do_pdos = .true.
2567
2568 IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
2569 "PRINT%CUBES"), cp_p_file)) do_cubes = .true.
2570
2571 IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
2572 "PRINT%RESTART_WFN"), cp_p_file)) do_wfn_restart = .true.
2573
2574 IF (.NOT. (do_pdos .OR. do_cubes .OR. do_wfn_restart)) RETURN
2575
2576 CALL timeset(routinen, handle)
2577
2578 !Initialization
2579 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
2580 qs_kind_set=qs_kind_set, para_env=para_env, blacs_env=blacs_env, &
2581 matrix_s=matrix_s, mos=mos)
2582
2583 SELECT CASE (ex_type)
2584 CASE (tddfpt_spin_cons)
2585 lr_evals => donor_state%sc_evals
2586 lr_coeffs => donor_state%sc_coeffs
2587 nspins = 2
2588 excite = "spincons"
2589 CASE (tddfpt_spin_flip)
2590 lr_evals => donor_state%sf_evals
2591 lr_coeffs => donor_state%sf_coeffs
2592 nspins = 2
2593 excite = "spinflip"
2594 CASE (tddfpt_singlet)
2595 lr_evals => donor_state%sg_evals
2596 lr_coeffs => donor_state%sg_coeffs
2597 nspins = 1
2598 excite = "singlet"
2599 CASE (tddfpt_triplet)
2600 lr_evals => donor_state%tp_evals
2601 lr_coeffs => donor_state%tp_coeffs
2602 nspins = 1
2603 excite = "triplet"
2604 END SELECT
2605
2606 SELECT CASE (donor_state%state_type)
2607 CASE (xas_1s_type)
2608 domo = "1s"
2609 CASE (xas_2s_type)
2610 domo = "2s"
2611 CASE (xas_2p_type)
2612 domo = "2p"
2613 END SELECT
2614
2615 ndo_mo = donor_state%ndo_mo
2616 ndo_so = ndo_mo*nspins
2617 nmo = SIZE(lr_evals)
2618 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
2619
2620 CALL cp_fm_struct_create(mo_struct, context=blacs_env, para_env=para_env, &
2621 nrow_global=nao, ncol_global=nmo)
2622 CALL cp_fm_create(mo_coeff, mo_struct)
2623
2624 !Dump the TDDFT excited state AMEW wavefunction into a file for restart in RTP
2625 IF (do_wfn_restart) THEN
2626 block
2627 TYPE(mo_set_type), DIMENSION(2) :: restart_mos
2628 IF (.NOT. (nspins == 1 .AND. donor_state%state_type == xas_1s_type)) THEN
2629 cpabort("RESTART.wfn file only available for RKS K-edge XAS spectroscopy")
2630 END IF
2631
2632 CALL section_vals_val_get(xas_tdp_section, "PRINT%RESTART_WFN%EXCITED_STATE_INDEX", n_rep_val=n_rep)
2633
2634 DO irep = 1, n_rep
2635 CALL section_vals_val_get(xas_tdp_section, "PRINT%RESTART_WFN%EXCITED_STATE_INDEX", &
2636 i_rep_val=irep, i_val=ex_state_idx)
2637 cpassert(ex_state_idx <= SIZE(lr_evals))
2638
2639 DO ispin = 1, 2
2640 CALL duplicate_mo_set(restart_mos(ispin), mos(1))
2641 ! Set the new occupation number in the case of spin-independent based calculaltion since the restart is spin-depedent
2642 IF (SIZE(mos) == 1) &
2643 restart_mos(ispin)%occupation_numbers = mos(1)%occupation_numbers/2
2644 END DO
2645
2646 CALL cp_fm_to_fm_submat(msource=lr_coeffs, mtarget=restart_mos(1)%mo_coeff, nrow=nao, &
2647 ncol=1, s_firstrow=1, s_firstcol=ex_state_idx, t_firstrow=1, &
2648 t_firstcol=donor_state%mo_indices(1, 1))
2649
2650 xas_mittle = 'xasat'//trim(adjustl(cp_to_string(donor_state%at_index)))//'_'//trim(domo)// &
2651 '_'//trim(excite)//'_idx'//trim(adjustl(cp_to_string(ex_state_idx)))
2652 output_unit = cp_print_key_unit_nr(logger, xas_tdp_section, "PRINT%RESTART_WFN", &
2653 extension=".wfn", file_status="REPLACE", &
2654 file_action="WRITE", file_form="UNFORMATTED", &
2655 middle_name=xas_mittle)
2656
2657 CALL write_mo_set_low(restart_mos, particle_set=particle_set, &
2658 qs_kind_set=qs_kind_set, ires=output_unit)
2659
2660 CALL cp_print_key_finished_output(output_unit, logger, xas_tdp_section, "PRINT%RESTART_WFN")
2661
2662 DO ispin = 1, 2
2663 CALL deallocate_mo_set(restart_mos(ispin))
2664 END DO
2665 END DO
2666 END block
2667 END IF
2668
2669 !PDOS related stuff
2670 IF (do_pdos) THEN
2671
2672 !If S^0.5 not yet stored, compute it once and for all
2673 IF (.NOT. ASSOCIATED(xas_tdp_env%matrix_shalf) .AND. do_pdos) THEN
2674 CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
2675 nrow_global=nao, ncol_global=nao)
2676 ALLOCATE (xas_tdp_env%matrix_shalf)
2677 CALL cp_fm_create(xas_tdp_env%matrix_shalf, fm_struct)
2678 CALL cp_fm_create(work_fm, fm_struct)
2679
2680 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, xas_tdp_env%matrix_shalf)
2681 CALL cp_fm_power(xas_tdp_env%matrix_shalf, work_fm, 0.5_dp, epsilon(0.0_dp), n_dependent)
2682
2683 CALL cp_fm_release(work_fm)
2684 CALL cp_fm_struct_release(fm_struct)
2685 END IF
2686
2687 !Giving some PDOS info
2688 output_unit = cp_logger_get_default_io_unit()
2689 IF (output_unit > 0) THEN
2690 WRITE (unit=output_unit, fmt="(/,T5,A,/,T5,A,/,T5,A)") &
2691 "Computing the PDOS of linear-response orbitals for spectral features analysis", &
2692 "Note: using standard PDOS routines => ignore mentions of KS states and MO ", &
2693 " occupation numbers. Eigenvalues in *.pdos files are excitations energies."
2694 END IF
2695
2696 !Check on NLUMO
2697 CALL section_vals_val_get(xas_tdp_section, "PRINT%PDOS%NLUMO", i_val=nlumo)
2698 IF (nlumo /= 0) THEN
2699 cpwarn("NLUMO is irrelevant for XAS_TDP PDOS. It was overwritten to 0.")
2700 END IF
2701 CALL section_vals_val_set(xas_tdp_section, "PRINT%PDOS%NLUMO", i_val=0)
2702 END IF
2703
2704 !CUBES related stuff
2705 IF (do_cubes) THEN
2706
2707 print_key => section_vals_get_subs_vals(xas_tdp_section, "PRINT%CUBES")
2708
2709 CALL section_vals_val_get(print_key, "CUBES_LU_BOUNDS", i_vals=bounds)
2710 ncubes = bounds(2) - bounds(1) + 1
2711 IF (ncubes > 0) THEN
2712 ALLOCATE (state_list(ncubes))
2713 DO ic = 1, ncubes
2714 state_list(ic) = bounds(1) + ic - 1
2715 END DO
2716 END IF
2717
2718 IF (.NOT. ASSOCIATED(state_list)) THEN
2719 CALL section_vals_val_get(print_key, "CUBES_LIST", n_rep_val=n_rep)
2720
2721 ncubes = 0
2722 DO irep = 1, n_rep
2723 NULLIFY (list)
2724 CALL section_vals_val_get(print_key, "CUBES_LIST", i_rep_val=irep, i_vals=list)
2725 IF (ASSOCIATED(list)) THEN
2726 CALL reallocate(state_list, 1, ncubes + SIZE(list))
2727 DO ic = 1, SIZE(list)
2728 state_list(ncubes + ic) = list(ic)
2729 END DO
2730 ncubes = ncubes + SIZE(list)
2731 END IF
2732 END DO
2733 END IF
2734
2735 IF (.NOT. ASSOCIATED(state_list)) THEN
2736 ncubes = 1
2737 ALLOCATE (state_list(1))
2738 state_list(1) = 1
2739 END IF
2740
2741 CALL section_vals_val_get(print_key, "APPEND", l_val=append_cube)
2742 pos = "REWIND"
2743 IF (append_cube) pos = "APPEND"
2744
2745 ALLOCATE (centers(6, ncubes))
2746 centers = 0.0_dp
2747
2748 END IF
2749
2750 !Loop over MOs and spin, one PDOS/CUBE for each
2751 DO ido_mo = 1, ndo_mo
2752 DO ispin = 1, nspins
2753
2754 !need to create a mo set for the LR-orbitals
2755 ALLOCATE (mo_set)
2756 CALL allocate_mo_set(mo_set, nao=nao, nmo=nmo, nelectron=nmo, n_el_f=real(nmo, dp), &
2757 maxocc=1.0_dp, flexible_electron_count=0.0_dp)
2758 CALL init_mo_set(mo_set, fm_ref=mo_coeff, name="PDOS XAS_TDP MOs")
2759 mo_set%eigenvalues(:) = lr_evals(:)
2760
2761 !get the actual coeff => most common case: closed-shell K-edge, can directly take lr_coeffs
2762 IF (nspins == 1 .AND. ndo_mo == 1) THEN
2763 CALL cp_fm_to_fm(lr_coeffs, mo_set%mo_coeff)
2764 ELSE
2765 DO imo = 1, nmo
2766 CALL cp_fm_to_fm_submat(msource=lr_coeffs, mtarget=mo_set%mo_coeff, &
2767 nrow=nao, ncol=1, s_firstrow=1, &
2768 s_firstcol=(imo - 1)*ndo_so + (ispin - 1)*ndo_mo + ido_mo, &
2769 t_firstrow=1, t_firstcol=imo)
2770 END DO
2771 END IF
2772
2773 !naming the output
2774 domon = domo
2775 IF (donor_state%state_type == xas_2p_type) domon = trim(domo)//trim(adjustl(cp_to_string(ido_mo)))
2776 xas_mittle = 'xasat'//trim(adjustl(cp_to_string(donor_state%at_index)))//'_'// &
2777 trim(domon)//'_'//trim(excite)
2778
2779 IF (do_pdos) THEN
2780 CALL calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
2781 qs_env, xas_tdp_section, ispin, xas_mittle, &
2782 external_matrix_shalf=xas_tdp_env%matrix_shalf)
2783 END IF
2784
2785 IF (do_cubes) THEN
2786 CALL qs_print_cubes(qs_env, mo_set%mo_coeff, ncubes, state_list, centers, &
2787 print_key=print_key, root=xas_mittle, ispin=ispin, &
2788 file_position=pos)
2789 END IF
2790
2791 !clean-up
2792 CALL deallocate_mo_set(mo_set)
2793 DEALLOCATE (mo_set)
2794
2795 END DO
2796 END DO
2797
2798 !clean-up
2799 CALL cp_fm_release(mo_coeff)
2800 CALL cp_fm_struct_release(mo_struct)
2801 IF (do_cubes) DEALLOCATE (centers, state_list)
2802
2803 CALL timestop(handle)
2804
2805 END SUBROUTINE xas_tdp_post
2806
2807! **************************************************************************************************
2808!> \brief Computed the LUMOs for the OT eigensolver guesses
2809!> \param xas_tdp_env ...
2810!> \param xas_tdp_control ...
2811!> \param qs_env ...
2812!> \note Uses stendard diagonalization. Do not use the stendard make_lumo subroutine as it uses
2813!> the OT eigensolver and there is no guarantee that it will converge fast
2814! **************************************************************************************************
2815 SUBROUTINE make_lumo_guess(xas_tdp_env, xas_tdp_control, qs_env)
2816
2817 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2818 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2819 TYPE(qs_environment_type), POINTER :: qs_env
2820
2821 CHARACTER(len=*), PARAMETER :: routinen = 'make_lumo_guess'
2822
2823 INTEGER :: handle, ispin, nao, nelec_spin(2), &
2824 nlumo(2), nocc(2), nspins
2825 LOGICAL :: do_os
2826 REAL(dp), ALLOCATABLE, DIMENSION(:) :: evals
2827 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2828 TYPE(cp_fm_struct_type), POINTER :: fm_struct, lumo_struct
2829 TYPE(cp_fm_type) :: amatrix, bmatrix, evecs, work_fm
2830 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
2831 TYPE(mp_para_env_type), POINTER :: para_env
2832
2833 NULLIFY (matrix_ks, matrix_s, para_env, blacs_env)
2834 NULLIFY (lumo_struct, fm_struct)
2835
2836 CALL timeset(routinen, handle)
2837
2838 do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
2839 nspins = 1; IF (do_os) nspins = 2
2840 ALLOCATE (xas_tdp_env%lumo_evecs(nspins))
2841 ALLOCATE (xas_tdp_env%lumo_evals(nspins))
2842 CALL get_qs_env(qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s, nelectron_spin=nelec_spin, &
2843 para_env=para_env, blacs_env=blacs_env)
2844 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
2845
2846 IF (do_os) THEN
2847 nlumo = nao - nelec_spin
2848 nocc = nelec_spin
2849 ELSE
2850 nlumo = nao - nelec_spin(1)/2
2851 nocc = nelec_spin(1)/2
2852 END IF
2853
2854 ALLOCATE (xas_tdp_env%ot_prec(nspins))
2855
2856 DO ispin = 1, nspins
2857
2858 !Going through fm to diagonalize
2859 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2860 nrow_global=nao, ncol_global=nao)
2861 CALL cp_fm_create(amatrix, fm_struct)
2862 CALL cp_fm_create(bmatrix, fm_struct)
2863 CALL cp_fm_create(evecs, fm_struct)
2864 CALL cp_fm_create(work_fm, fm_struct)
2865 ALLOCATE (evals(nao))
2866 ALLOCATE (xas_tdp_env%lumo_evals(ispin)%array(nlumo(ispin)))
2867
2868 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, amatrix)
2869 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, bmatrix)
2870
2871 !The actual diagonalization through Cholesky decomposition
2872 CALL cp_fm_geeig(amatrix, bmatrix, evecs, evals, work_fm)
2873
2874 !Storing results
2875 CALL cp_fm_struct_create(lumo_struct, para_env=para_env, context=blacs_env, &
2876 nrow_global=nao, ncol_global=nlumo(ispin))
2877 CALL cp_fm_create(xas_tdp_env%lumo_evecs(ispin), lumo_struct)
2878
2879 CALL cp_fm_to_fm_submat(evecs, xas_tdp_env%lumo_evecs(ispin), nrow=nao, &
2880 ncol=nlumo(ispin), s_firstrow=1, s_firstcol=nocc(ispin) + 1, &
2881 t_firstrow=1, t_firstcol=1)
2882
2883 xas_tdp_env%lumo_evals(ispin)%array(1:nlumo(ispin)) = evals(nocc(ispin) + 1:nao)
2884
2885 CALL build_ot_spin_prec(evecs, evals, ispin, xas_tdp_env, xas_tdp_control, qs_env)
2886
2887 !clean-up
2888 CALL cp_fm_release(amatrix)
2889 CALL cp_fm_release(bmatrix)
2890 CALL cp_fm_release(evecs)
2891 CALL cp_fm_release(work_fm)
2892 CALL cp_fm_struct_release(fm_struct)
2893 CALL cp_fm_struct_release(lumo_struct)
2894 DEALLOCATE (evals)
2895 END DO
2896
2897 CALL timestop(handle)
2898
2899 END SUBROUTINE make_lumo_guess
2900
2901! **************************************************************************************************
2902!> \brief Builds a preconditioner for the OT eigensolver, based on some heurstics that prioritize
2903!> LUMOs with lower eigenvalues
2904!> \param evecs all the ground state eigenvectors
2905!> \param evals all the ground state eigenvalues
2906!> \param ispin ...
2907!> \param xas_tdp_env ...
2908!> \param xas_tdp_control ...
2909!> \param qs_env ...
2910!> \note assumes that the preconditioner matrix array is allocated
2911! **************************************************************************************************
2912 SUBROUTINE build_ot_spin_prec(evecs, evals, ispin, xas_tdp_env, xas_tdp_control, qs_env)
2913
2914 TYPE(cp_fm_type), INTENT(IN) :: evecs
2915 REAL(dp), DIMENSION(:) :: evals
2916 INTEGER :: ispin
2917 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2918 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2919 TYPE(qs_environment_type), POINTER :: qs_env
2920
2921 CHARACTER(len=*), PARAMETER :: routinen = 'build_ot_spin_prec'
2922
2923 INTEGER :: handle, nao, nelec_spin(2), nguess, &
2924 nocc, nspins
2925 LOGICAL :: do_os
2926 REAL(dp) :: shift
2927 REAL(dp), ALLOCATABLE, DIMENSION(:) :: scaling
2928 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2929 TYPE(cp_fm_type) :: fm_prec, work_fm
2930 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2931 TYPE(mp_para_env_type), POINTER :: para_env
2932
2933 NULLIFY (fm_struct, para_env, matrix_s)
2934
2935 CALL timeset(routinen, handle)
2936
2937 do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
2938 CALL get_qs_env(qs_env, para_env=para_env, nelectron_spin=nelec_spin, matrix_s=matrix_s)
2939 CALL cp_fm_get_info(evecs, nrow_global=nao, matrix_struct=fm_struct)
2940 CALL cp_fm_create(fm_prec, fm_struct)
2941 ALLOCATE (scaling(nao))
2942 nocc = nelec_spin(1)/2
2943 nspins = 1
2944 IF (do_os) THEN
2945 nocc = nelec_spin(ispin)
2946 nspins = 2
2947 END IF
2948
2949 !rough estimate of the number of required evals
2950 nguess = nao - nocc
2951 IF (xas_tdp_control%n_excited > 0 .AND. xas_tdp_control%n_excited < nguess) THEN
2952 nguess = xas_tdp_control%n_excited/nspins
2953 ELSE IF (xas_tdp_control%e_range > 0.0_dp) THEN
2954 nguess = count(evals(nocc + 1:nao) - evals(nocc + 1) <= xas_tdp_control%e_range)
2955 END IF
2956
2957 !Give max weight to the first LUMOs
2958 scaling(nocc + 1:nocc + nguess) = 100.0_dp
2959 !Then gradually decrease weight
2960 shift = evals(nocc + 1) - 0.01_dp
2961 scaling(nocc + nguess:nao) = 1.0_dp/(evals(nocc + nguess:nao) - shift)
2962 !HOMOs do not matter, but need well behaved matrix
2963 scaling(1:nocc) = 1.0_dp
2964
2965 !Building the precond as an fm
2966 CALL cp_fm_create(work_fm, fm_struct)
2967
2968 CALL cp_fm_copy_general(evecs, work_fm, para_env)
2969 CALL cp_fm_column_scale(work_fm, scaling)
2970
2971 CALL parallel_gemm('N', 'T', nao, nao, nao, 1.0_dp, work_fm, evecs, 0.0_dp, fm_prec)
2972
2973 !Copy into dbcsr format
2974 ALLOCATE (xas_tdp_env%ot_prec(ispin)%matrix)
2975 CALL dbcsr_create(xas_tdp_env%ot_prec(ispin)%matrix, template=matrix_s(1)%matrix, name="OT_PREC")
2976 CALL copy_fm_to_dbcsr(fm_prec, xas_tdp_env%ot_prec(ispin)%matrix)
2977 CALL dbcsr_filter(xas_tdp_env%ot_prec(ispin)%matrix, xas_tdp_control%eps_filter)
2978
2979 CALL cp_fm_release(work_fm)
2980 CALL cp_fm_release(fm_prec)
2981
2982 CALL timestop(handle)
2983
2984 END SUBROUTINE build_ot_spin_prec
2985
2986! **************************************************************************************************
2987!> \brief Prints GW2X corrected ionization potentials to main output file, including SOC splitting
2988!> \param donor_state ...
2989!> \param xas_tdp_env ...
2990!> \param xas_tdp_control ...
2991!> \param qs_env ...
2992! **************************************************************************************************
2993 SUBROUTINE print_xps(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
2994
2995 TYPE(donor_state_type), POINTER :: donor_state
2996 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2997 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2998 TYPE(qs_environment_type), POINTER :: qs_env
2999
3000 INTEGER :: ido_mo, ispin, nspins, output_unit
3001 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ips, soc_shifts
3002
3003 output_unit = cp_logger_get_default_io_unit()
3004
3005 nspins = 1; IF (xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks) nspins = 2
3006
3007 ALLOCATE (ips(SIZE(donor_state%gw2x_evals, 1), SIZE(donor_state%gw2x_evals, 2)))
3008 ips(:, :) = donor_state%gw2x_evals(:, :)
3009
3010 !IPs in PBCs cannot be trusted because of a lack of a potential reference
3011 IF (.NOT. xas_tdp_control%is_periodic) THEN
3012
3013 !Apply SOC splitting
3014 IF (donor_state%ndo_mo > 1) THEN
3015 CALL get_soc_splitting(soc_shifts, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
3016 ips(:, :) = ips(:, :) + soc_shifts
3017
3018 IF (output_unit > 0) THEN
3019 WRITE (output_unit, fmt="(/,T5,A,F23.6)") &
3020 "Ionization potentials for XPS (GW2X + SOC): ", -ips(1, 1)*evolt
3021
3022 DO ispin = 1, nspins
3023 DO ido_mo = 1, donor_state%ndo_mo
3024
3025 IF (ispin == 1 .AND. ido_mo == 1) cycle
3026
3027 WRITE (output_unit, fmt="(T5,A,F23.6)") &
3028 " ", -ips(ido_mo, ispin)*evolt
3029
3030 END DO
3031 END DO
3032 END IF
3033
3034 ELSE
3035
3036 ! No SOC, only 1 donor MO per spin
3037 IF (output_unit > 0) THEN
3038 WRITE (output_unit, fmt="(/,T5,A,F29.6)") &
3039 "Ionization potentials for XPS (GW2X): ", -ips(1, 1)*evolt
3040
3041 IF (nspins == 2) THEN
3042 WRITE (output_unit, fmt="(T5,A,F29.6)") &
3043 " ", -ips(1, 2)*evolt
3044 END IF
3045 END IF
3046
3047 END IF
3048 END IF
3049
3050 END SUBROUTINE print_xps
3051
3052! **************************************************************************************************
3053!> \brief Prints the excitation energies and the oscillator strengths for a given donor_state in a file
3054!> \param donor_state the donor_state to print
3055!> \param xas_tdp_env ...
3056!> \param xas_tdp_control ...
3057!> \param xas_tdp_section ...
3058! **************************************************************************************************
3059 SUBROUTINE print_xas_tdp_to_file(donor_state, xas_tdp_env, xas_tdp_control, xas_tdp_section)
3060
3061 TYPE(donor_state_type), POINTER :: donor_state
3062 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
3063 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
3064 TYPE(section_vals_type), POINTER :: xas_tdp_section
3065
3066 INTEGER :: i, output_unit, xas_tdp_unit
3067 TYPE(cp_logger_type), POINTER :: logger
3068
3069 NULLIFY (logger)
3070 logger => cp_get_default_logger()
3071
3072 xas_tdp_unit = cp_print_key_unit_nr(logger, xas_tdp_section, "PRINT%SPECTRUM", &
3073 extension=".spectrum", file_position="APPEND", &
3074 file_action="WRITE", file_form="FORMATTED")
3075
3076 output_unit = cp_logger_get_default_io_unit()
3077
3078 IF (output_unit > 0) THEN
3079 WRITE (output_unit, fmt="(/,T5,A,/)") &
3080 "Calculations done: "
3081 END IF
3082
3083 IF (xas_tdp_control%do_spin_cons) THEN
3084 IF (xas_tdp_unit > 0) THEN
3085
3086! Printing the general donor state information
3087 WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3088 "==================================================================================", &
3089 "XAS TDP open-shell spin-conserving (no SOC) excitations for DONOR STATE: ", &
3090 xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3091 "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3092 donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3093 "=================================================================================="
3094
3095! Simply dump the excitation energies/ oscillator strength as they come
3096
3097 IF (xas_tdp_control%do_quad) THEN
3098 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3099 " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3100 DO i = 1, SIZE(donor_state%sc_evals)
3101 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3102 i, donor_state%sc_evals(i)*evolt, donor_state%osc_str(i, 4), &
3103 donor_state%quad_osc_str(i)
3104 END DO
3105 ELSE IF (xas_tdp_control%xyz_dip) THEN
3106 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3107 " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3108 DO i = 1, SIZE(donor_state%sc_evals)
3109 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3110 i, donor_state%sc_evals(i)*evolt, donor_state%osc_str(i, 4), &
3111 donor_state%osc_str(i, 1), donor_state%osc_str(i, 2), donor_state%osc_str(i, 3)
3112 END DO
3113 ELSE IF (xas_tdp_control%spin_dip) THEN
3114 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3115 " Index Excitation energy (eV) fosc dipole (a.u.) alpha-comp beta-comp"
3116 DO i = 1, SIZE(donor_state%sc_evals)
3117 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6)") &
3118 i, donor_state%sc_evals(i)*evolt, donor_state%osc_str(i, 4), &
3119 donor_state%alpha_osc(i, 4), donor_state%beta_osc(i, 4)
3120 END DO
3121 ELSE
3122 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3123 " Index Excitation energy (eV) fosc dipole (a.u.)"
3124 DO i = 1, SIZE(donor_state%sc_evals)
3125 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3126 i, donor_state%sc_evals(i)*evolt, donor_state%osc_str(i, 4)
3127 END DO
3128 END IF
3129
3130 WRITE (xas_tdp_unit, fmt="(A,/)") " "
3131 END IF !xas_tdp_unit > 0
3132
3133 IF (output_unit > 0) THEN
3134 WRITE (output_unit, fmt="(T5,A,F17.6)") &
3135 "First spin-conserving XAS excitation energy (eV): ", donor_state%sc_evals(1)*evolt
3136 END IF
3137
3138 END IF ! do_spin_cons
3139
3140 IF (xas_tdp_control%do_spin_flip) THEN
3141 IF (xas_tdp_unit > 0) THEN
3142
3143! Printing the general donor state information
3144 WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3145 "==================================================================================", &
3146 "XAS TDP open-shell spin-flip (no SOC) excitations for DONOR STATE: ", &
3147 xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3148 "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3149 donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3150 "=================================================================================="
3151
3152! Simply dump the excitation energies/ oscillator strength as they come
3153
3154 IF (xas_tdp_control%do_quad) THEN
3155 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3156 " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3157 DO i = 1, SIZE(donor_state%sf_evals)
3158 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3159 i, donor_state%sf_evals(i)*evolt, 0.0_dp, 0.0_dp !spin-forbidden !
3160 END DO
3161 ELSE IF (xas_tdp_control%xyz_dip) THEN
3162 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3163 " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3164 DO i = 1, SIZE(donor_state%sf_evals)
3165 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3166 i, donor_state%sf_evals(i)*evolt, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp
3167 END DO
3168 ELSE IF (xas_tdp_control%spin_dip) THEN
3169 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3170 " Index Excitation energy (eV) fosc dipole (a.u.) alpha-comp beta-comp"
3171 DO i = 1, SIZE(donor_state%sf_evals)
3172 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6)") &
3173 i, donor_state%sf_evals(i)*evolt, 0.0_dp, 0.0_dp, 0.0_dp
3174 END DO
3175 ELSE
3176 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3177 " Index Excitation energy (eV) fosc dipole (a.u.)"
3178 DO i = 1, SIZE(donor_state%sf_evals)
3179 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3180 i, donor_state%sf_evals(i)*evolt, 0.0_dp
3181 END DO
3182 END IF
3183
3184 WRITE (xas_tdp_unit, fmt="(A,/)") " "
3185 END IF !xas_tdp_unit
3186
3187 IF (output_unit > 0) THEN
3188 WRITE (output_unit, fmt="(T5,A,F23.6)") &
3189 "First spin-flip XAS excitation energy (eV): ", donor_state%sf_evals(1)*evolt
3190 END IF
3191 END IF ! do_spin_flip
3192
3193 IF (xas_tdp_control%do_singlet) THEN
3194 IF (xas_tdp_unit > 0) THEN
3195
3196! Printing the general donor state information
3197 WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3198 "==================================================================================", &
3199 "XAS TDP singlet excitations (no SOC) for DONOR STATE: ", &
3200 xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3201 "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3202 donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3203 "=================================================================================="
3204
3205! Simply dump the excitation energies/ oscillator strength as they come
3206
3207 IF (xas_tdp_control%do_quad) THEN
3208 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3209 " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3210 DO i = 1, SIZE(donor_state%sg_evals)
3211 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3212 i, donor_state%sg_evals(i)*evolt, donor_state%osc_str(i, 4), &
3213 donor_state%quad_osc_str(i)
3214 END DO
3215 ELSE IF (xas_tdp_control%xyz_dip) THEN
3216 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3217 " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3218 DO i = 1, SIZE(donor_state%sg_evals)
3219 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3220 i, donor_state%sg_evals(i)*evolt, donor_state%osc_str(i, 4), &
3221 donor_state%osc_str(i, 1), donor_state%osc_str(i, 2), donor_state%osc_str(i, 3)
3222 END DO
3223 ELSE
3224 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3225 " Index Excitation energy (eV) fosc dipole (a.u.)"
3226 DO i = 1, SIZE(donor_state%sg_evals)
3227 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3228 i, donor_state%sg_evals(i)*evolt, donor_state%osc_str(i, 4)
3229 END DO
3230 END IF
3231
3232 WRITE (xas_tdp_unit, fmt="(A,/)") " "
3233 END IF !xas_tdp_unit
3234
3235 IF (output_unit > 0) THEN
3236 WRITE (output_unit, fmt="(T5,A,F25.6)") &
3237 "First singlet XAS excitation energy (eV): ", donor_state%sg_evals(1)*evolt
3238 END IF
3239 END IF ! do_singlet
3240
3241 IF (xas_tdp_control%do_triplet) THEN
3242 IF (xas_tdp_unit > 0) THEN
3243
3244! Printing the general donor state information
3245 WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3246 "==================================================================================", &
3247 "XAS TDP triplet excitations (no SOC) for DONOR STATE: ", &
3248 xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3249 "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3250 donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3251 "=================================================================================="
3252
3253! Simply dump the excitation energies/ oscillator strength as they come
3254
3255 IF (xas_tdp_control%do_quad) THEN
3256 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3257 " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3258 DO i = 1, SIZE(donor_state%tp_evals)
3259 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3260 i, donor_state%tp_evals(i)*evolt, 0.0_dp, 0.0_dp !spin-forbidden !
3261 END DO
3262 ELSE IF (xas_tdp_control%xyz_dip) THEN
3263 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3264 " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3265 DO i = 1, SIZE(donor_state%tp_evals)
3266 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3267 i, donor_state%tp_evals(i)*evolt, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp
3268 END DO
3269 ELSE IF (xas_tdp_control%spin_dip) THEN
3270 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3271 " Index Excitation energy (eV) fosc dipole (a.u.) alpha-comp beta-comp"
3272 DO i = 1, SIZE(donor_state%tp_evals)
3273 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6)") &
3274 i, donor_state%tp_evals(i)*evolt, 0.0_dp, 0.0_dp, 0.0_dp
3275 END DO
3276 ELSE
3277 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3278 " Index Excitation energy (eV) fosc dipole (a.u.)"
3279 DO i = 1, SIZE(donor_state%tp_evals)
3280 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3281 i, donor_state%tp_evals(i)*evolt, 0.0_dp
3282 END DO
3283 END IF
3284
3285 WRITE (xas_tdp_unit, fmt="(A,/)") " "
3286 END IF !xas_tdp_unit
3287
3288 IF (output_unit > 0) THEN
3289 WRITE (output_unit, fmt="(T5,A,F25.6)") &
3290 "First triplet XAS excitation energy (eV): ", donor_state%tp_evals(1)*evolt
3291 END IF
3292 END IF ! do_triplet
3293
3294 IF (xas_tdp_control%do_soc .AND. donor_state%state_type == xas_2p_type) THEN
3295 IF (xas_tdp_unit > 0) THEN
3296
3297! Printing the general donor state information
3298 WRITE (xas_tdp_unit, fmt="(A,/,A,A,A/,A,I5,A,I5,A,A,/,A)") &
3299 "==================================================================================", &
3300 "XAS TDP excitations after spin-orbit coupling for DONOR STATE: ", &
3301 xas_tdp_env%state_type_char(donor_state%state_type), ",", &
3302 "from EXCITED ATOM: ", donor_state%at_index, ", of KIND (index/symbol): ", &
3303 donor_state%kind_index, "/", trim(donor_state%at_symbol), &
3304 "=================================================================================="
3305
3306! Simply dump the excitation energies/ oscillator strength as they come
3307 IF (xas_tdp_control%do_quad) THEN
3308 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3309 " Index Excitation energy (eV) fosc dipole (a.u.) fosc quadrupole (a.u.)"
3310 DO i = 1, SIZE(donor_state%soc_evals)
3311 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F25.6)") &
3312 i, donor_state%soc_evals(i)*evolt, donor_state%soc_osc_str(i, 4), &
3313 donor_state%soc_quad_osc_str(i)
3314 END DO
3315 ELSE IF (xas_tdp_control%xyz_dip) THEN
3316 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3317 " Index Excitation energy (eV) fosc dipole (a.u.) x-component y-component z-component"
3318 DO i = 1, SIZE(donor_state%soc_evals)
3319 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6,F14.6,F14.6,F14.6)") &
3320 i, donor_state%soc_evals(i)*evolt, donor_state%soc_osc_str(i, 4), &
3321 donor_state%soc_osc_str(i, 1), donor_state%soc_osc_str(i, 2), donor_state%soc_osc_str(i, 3)
3322 END DO
3323 ELSE
3324 WRITE (xas_tdp_unit, fmt="(T3,A)") &
3325 " Index Excitation energy (eV) fosc dipole (a.u.)"
3326 DO i = 1, SIZE(donor_state%soc_evals)
3327 WRITE (xas_tdp_unit, fmt="(T3,I6,F27.6,F22.6)") &
3328 i, donor_state%soc_evals(i)*evolt, donor_state%soc_osc_str(i, 4)
3329 END DO
3330 END IF
3331
3332 WRITE (xas_tdp_unit, fmt="(A,/)") " "
3333 END IF !xas_tdp_unit
3334
3335 IF (output_unit > 0) THEN
3336 WRITE (output_unit, fmt="(T5,A,F29.6)") &
3337 "First SOC XAS excitation energy (eV): ", donor_state%soc_evals(1)*evolt
3338 END IF
3339 END IF !do_soc
3340
3341 CALL cp_print_key_finished_output(xas_tdp_unit, logger, xas_tdp_section, "PRINT%SPECTRUM")
3342
3343 END SUBROUTINE print_xas_tdp_to_file
3344
3345! **************************************************************************************************
3346!> \brief Prints the donor_state and excitation_type info into a RESTART file for cheap PDOS and/or
3347!> CUBE printing without expensive computation
3348!> \param ex_type singlet, triplet, etc.
3349!> \param donor_state ...
3350!> \param xas_tdp_section ...
3351!> \param qs_env ...
3352! **************************************************************************************************
3353 SUBROUTINE write_donor_state_restart(ex_type, donor_state, xas_tdp_section, qs_env)
3354
3355 INTEGER, INTENT(IN) :: ex_type
3356 TYPE(donor_state_type), POINTER :: donor_state
3357 TYPE(section_vals_type), POINTER :: xas_tdp_section
3358 TYPE(qs_environment_type), POINTER :: qs_env
3359
3360 CHARACTER(len=*), PARAMETER :: routinen = 'write_donor_state_restart'
3361
3362 CHARACTER(len=default_path_length) :: filename
3363 CHARACTER(len=default_string_length) :: domo, excite, my_middle
3364 INTEGER :: ex_atom, handle, ispin, nao, ndo_mo, &
3365 nex, nspins, output_unit, rst_unit, &
3366 state_type
3367 INTEGER, DIMENSION(:, :), POINTER :: mo_indices
3368 LOGICAL :: do_print
3369 REAL(dp), DIMENSION(:), POINTER :: lr_evals
3370 TYPE(cp_fm_type), POINTER :: lr_coeffs
3371 TYPE(cp_logger_type), POINTER :: logger
3372 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3373 TYPE(section_vals_type), POINTER :: print_key
3374
3375 NULLIFY (logger, lr_coeffs, lr_evals, print_key, mos)
3376
3377 !Initialization
3378 logger => cp_get_default_logger()
3379 do_print = .false.
3380 IF (btest(cp_print_key_should_output(logger%iter_info, xas_tdp_section, &
3381 "PRINT%RESTART", used_print_key=print_key), cp_p_file)) do_print = .true.
3382
3383 IF (.NOT. do_print) RETURN
3384
3385 CALL timeset(routinen, handle)
3386
3387 output_unit = cp_logger_get_default_io_unit()
3388
3389 !Get general info
3390 SELECT CASE (ex_type)
3391 CASE (tddfpt_spin_cons)
3392 lr_evals => donor_state%sc_evals
3393 lr_coeffs => donor_state%sc_coeffs
3394 excite = "spincons"
3395 nspins = 2
3396 CASE (tddfpt_spin_flip)
3397 lr_evals => donor_state%sf_evals
3398 lr_coeffs => donor_state%sf_coeffs
3399 excite = "spinflip"
3400 nspins = 2
3401 CASE (tddfpt_singlet)
3402 lr_evals => donor_state%sg_evals
3403 lr_coeffs => donor_state%sg_coeffs
3404 excite = "singlet"
3405 nspins = 1
3406 CASE (tddfpt_triplet)
3407 lr_evals => donor_state%tp_evals
3408 lr_coeffs => donor_state%tp_coeffs
3409 excite = "triplet"
3410 nspins = 1
3411 END SELECT
3412
3413 SELECT CASE (donor_state%state_type)
3414 CASE (xas_1s_type)
3415 domo = "1s"
3416 CASE (xas_2s_type)
3417 domo = "2s"
3418 CASE (xas_2p_type)
3419 domo = "2p"
3420 END SELECT
3421
3422 ndo_mo = donor_state%ndo_mo
3423 nex = SIZE(lr_evals)
3424 CALL cp_fm_get_info(lr_coeffs, nrow_global=nao)
3425 state_type = donor_state%state_type
3426 ex_atom = donor_state%at_index
3427 mo_indices => donor_state%mo_indices
3428
3429 !Opening restart file
3430 rst_unit = -1
3431 my_middle = 'xasat'//trim(adjustl(cp_to_string(ex_atom)))//'_'//trim(domo)//'_'//trim(excite)
3432 rst_unit = cp_print_key_unit_nr(logger, xas_tdp_section, "PRINT%RESTART", extension=".rst", &
3433 file_status="REPLACE", file_action="WRITE", &
3434 file_form="UNFORMATTED", middle_name=trim(my_middle))
3435
3436 filename = cp_print_key_generate_filename(logger, print_key, middle_name=trim(my_middle), &
3437 extension=".rst", my_local=.false.)
3438
3439 IF (output_unit > 0) THEN
3440 WRITE (unit=output_unit, fmt="(/,T5,A,/T5,A,A,A)") &
3441 "Linear-response orbitals and excitation energies are written in: ", &
3442 '"', trim(filename), '"'
3443 END IF
3444
3445 !Writing
3446 IF (rst_unit > 0) THEN
3447 WRITE (rst_unit) ex_atom, state_type, ndo_mo, ex_type
3448 WRITE (rst_unit) nao, nex, nspins
3449 WRITE (rst_unit) mo_indices(:, :)
3450 WRITE (rst_unit) lr_evals(:)
3451 END IF
3452 CALL cp_fm_write_unformatted(lr_coeffs, rst_unit)
3453
3454 !The MOs as well (because the may have been localized)
3455 CALL get_qs_env(qs_env, mos=mos)
3456 DO ispin = 1, nspins
3457 CALL cp_fm_write_unformatted(mos(ispin)%mo_coeff, rst_unit)
3458 END DO
3459
3460 !closing
3461 CALL cp_print_key_finished_output(rst_unit, logger, xas_tdp_section, "PRINT%RESTART")
3462
3463 CALL timestop(handle)
3464
3465 END SUBROUTINE write_donor_state_restart
3466
3467! **************************************************************************************************
3468!> \brief Reads donor_state info from a restart file
3469!> \param donor_state the pre-allocated donor_state
3470!> \param ex_type the excitations stored in this specific file
3471!> \param filename the restart file to read from
3472!> \param qs_env ...
3473! **************************************************************************************************
3474 SUBROUTINE read_donor_state_restart(donor_state, ex_type, filename, qs_env)
3475
3476 TYPE(donor_state_type), POINTER :: donor_state
3477 INTEGER, INTENT(OUT) :: ex_type
3478 CHARACTER(len=*), INTENT(IN) :: filename
3479 TYPE(qs_environment_type), POINTER :: qs_env
3480
3481 CHARACTER(len=*), PARAMETER :: routinen = 'read_donor_state_restart'
3482
3483 INTEGER :: handle, ispin, nao, nex, nspins, &
3484 output_unit, read_params(7), rst_unit
3485 INTEGER, DIMENSION(:, :), POINTER :: mo_indices
3486 LOGICAL :: file_exists
3487 REAL(dp), DIMENSION(:), POINTER :: lr_evals
3488 TYPE(cp_blacs_env_type), POINTER :: blacs_env
3489 TYPE(cp_fm_struct_type), POINTER :: fm_struct
3490 TYPE(cp_fm_type), POINTER :: lr_coeffs
3491 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3492 TYPE(mp_comm_type) :: group
3493 TYPE(mp_para_env_type), POINTER :: para_env
3494
3495 NULLIFY (lr_evals, lr_coeffs, para_env, fm_struct, blacs_env, mos)
3496
3497 CALL timeset(routinen, handle)
3498
3499 output_unit = cp_logger_get_default_io_unit()
3500 cpassert(ASSOCIATED(donor_state))
3501 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
3502 group = para_env
3503
3504 file_exists = .false.
3505 rst_unit = -1
3506
3507 IF (para_env%is_source()) THEN
3508
3509 INQUIRE (file=filename, exist=file_exists)
3510 IF (.NOT. file_exists) cpabort("Trying to read non-existing XAS_TDP restart file")
3511
3512 CALL open_file(file_name=trim(filename), file_action="READ", file_form="UNFORMATTED", &
3513 file_position="REWIND", file_status="OLD", unit_number=rst_unit)
3514 END IF
3515
3516 IF (output_unit > 0) THEN
3517 WRITE (unit=output_unit, fmt="(/,T5,A,/,T5,A,A,A)") &
3518 "Reading linear-response orbitals and excitation energies from file: ", &
3519 '"', filename, '"'
3520 END IF
3521
3522 !read general params
3523 IF (rst_unit > 0) THEN
3524 READ (rst_unit) read_params(1:4)
3525 READ (rst_unit) read_params(5:7)
3526 END IF
3527 CALL group%bcast(read_params)
3528 donor_state%at_index = read_params(1)
3529 donor_state%state_type = read_params(2)
3530 donor_state%ndo_mo = read_params(3)
3531 ex_type = read_params(4)
3532 nao = read_params(5)
3533 nex = read_params(6)
3534 nspins = read_params(7)
3535
3536 ALLOCATE (mo_indices(donor_state%ndo_mo, nspins))
3537 IF (rst_unit > 0) THEN
3538 READ (rst_unit) mo_indices(1:donor_state%ndo_mo, 1:nspins)
3539 END IF
3540 CALL group%bcast(mo_indices)
3541 donor_state%mo_indices => mo_indices
3542
3543 !read evals
3544 ALLOCATE (lr_evals(nex))
3545 IF (rst_unit > 0) READ (rst_unit) lr_evals(1:nex)
3546 CALL group%bcast(lr_evals)
3547
3548 !read evecs
3549 CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
3550 nrow_global=nao, ncol_global=nex*donor_state%ndo_mo*nspins)
3551 ALLOCATE (lr_coeffs)
3552 CALL cp_fm_create(lr_coeffs, fm_struct)
3553 CALL cp_fm_read_unformatted(lr_coeffs, rst_unit)
3554 CALL cp_fm_struct_release(fm_struct)
3555
3556 !read MO coeffs and replace in qs_env
3557 CALL get_qs_env(qs_env, mos=mos)
3558 DO ispin = 1, nspins
3559 CALL cp_fm_read_unformatted(mos(ispin)%mo_coeff, rst_unit)
3560 END DO
3561
3562 !closing file
3563 IF (para_env%is_source()) THEN
3564 CALL close_file(unit_number=rst_unit)
3565 END IF
3566
3567 !case study on excitation type
3568 SELECT CASE (ex_type)
3569 CASE (tddfpt_spin_cons)
3570 donor_state%sc_evals => lr_evals
3571 donor_state%sc_coeffs => lr_coeffs
3572 CASE (tddfpt_spin_flip)
3573 donor_state%sf_evals => lr_evals
3574 donor_state%sf_coeffs => lr_coeffs
3575 CASE (tddfpt_singlet)
3576 donor_state%sg_evals => lr_evals
3577 donor_state%sg_coeffs => lr_coeffs
3578 CASE (tddfpt_triplet)
3579 donor_state%tp_evals => lr_evals
3580 donor_state%tp_coeffs => lr_coeffs
3581 END SELECT
3582
3583 CALL timestop(handle)
3584
3585 END SUBROUTINE read_donor_state_restart
3586
3587! **************************************************************************************************
3588!> \brief Checks whether this is a restart calculation and runs it if so
3589!> \param rst_filename the file to read for restart
3590!> \param xas_tdp_section ...
3591!> \param qs_env ...
3592! **************************************************************************************************
3593 SUBROUTINE restart_calculation(rst_filename, xas_tdp_section, qs_env)
3594
3595 CHARACTER(len=*), INTENT(IN) :: rst_filename
3596 TYPE(section_vals_type), POINTER :: xas_tdp_section
3597 TYPE(qs_environment_type), POINTER :: qs_env
3598
3599 INTEGER :: ex_type
3600 TYPE(donor_state_type), POINTER :: donor_state
3601 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
3602
3603 NULLIFY (xas_tdp_env, donor_state)
3604
3605 !create a donor_state that we fill with the information we read
3606 ALLOCATE (donor_state)
3607 CALL donor_state_create(donor_state)
3608 CALL read_donor_state_restart(donor_state, ex_type, rst_filename, qs_env)
3609
3610 !create a dummy xas_tdp_env and compute the post XAS_TDP stuff
3611 CALL xas_tdp_env_create(xas_tdp_env)
3612 CALL xas_tdp_post(ex_type, donor_state, xas_tdp_env, xas_tdp_section, qs_env)
3613
3614 !clean-up
3615 CALL xas_tdp_env_release(xas_tdp_env)
3616 CALL free_ds_memory(donor_state)
3617 DEALLOCATE (donor_state%mo_indices)
3618 DEALLOCATE (donor_state)
3619
3620 END SUBROUTINE restart_calculation
3621
3622END 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:214
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:67
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)
...
Definition basis_set_types.F:631
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:2646
basis_set_types::deallocate_sto_basis_set
subroutine, public deallocate_sto_basis_set(sto_basis_set)
...
Definition basis_set_types.F:2077
basis_set_types::allocate_sto_basis_set
subroutine, public allocate_sto_basis_set(sto_basis_set)
...
Definition basis_set_types.F:2063
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:2367
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:2154
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:1089
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:1086
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:1194
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1132
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_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:552
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:214
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:1879
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:997
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:443
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_fm_diag.F:1355
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:1537
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:1009
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2146
cp_fm_types::cp_fm_read_unformatted
subroutine, public cp_fm_read_unformatted(fm, unit)
...
Definition cp_fm_types.F:2377
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:1472
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:528
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:894
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
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:1063
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::xas_1s_type
integer, parameter, public xas_1s_type
Definition input_constants.F:681
input_constants::do_admm_purify_none
integer, parameter, public do_admm_purify_none
Definition input_constants.F:858
input_constants::do_loc_none
integer, parameter, public do_loc_none
Definition input_constants.F:453
input_constants::op_loc_berry
integer, parameter, public op_loc_berry
Definition input_constants.F:468
input_constants::xas_not_excited
integer, parameter, public xas_not_excited
Definition input_constants.F:681
input_constants::tddfpt_singlet
integer, parameter, public tddfpt_singlet
Definition input_constants.F:577
input_constants::xas_2p_type
integer, parameter, public xas_2p_type
Definition input_constants.F:681
input_constants::xas_dip_len
integer, parameter, public xas_dip_len
Definition input_constants.F:692
input_constants::xas_dip_vel
integer, parameter, public xas_dip_vel
Definition input_constants.F:692
input_constants::tddfpt_triplet
integer, parameter, public tddfpt_triplet
Definition input_constants.F:577
input_constants::xas_2s_type
integer, parameter, public xas_2s_type
Definition input_constants.F:681
input_constants::tddfpt_spin_flip
integer, parameter, public tddfpt_spin_flip
Definition input_constants.F:577
input_constants::xas_tdp_by_kind
integer, parameter, public xas_tdp_by_kind
Definition input_constants.F:698
input_constants::do_admm_purify_cauchy_subspace
integer, parameter, public do_admm_purify_cauchy_subspace
Definition input_constants.F:858
input_constants::state_loc_list
integer, parameter, public state_loc_list
Definition input_constants.F:468
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:811
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:811
input_constants::xas_tdp_by_index
integer, parameter, public xas_tdp_by_index
Definition input_constants.F:698
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:811
input_constants::do_admm_purify_mo_diag
integer, parameter, public do_admm_purify_mo_diag
Definition input_constants.F:858
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:811
input_constants::tddfpt_spin_cons
integer, parameter, public tddfpt_spin_cons
Definition input_constants.F:577
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:1236
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:136
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::get_diag
pure real(kind=dp) function, dimension(min(size(a, 1), size(a, 2))), public get_diag(a)
Return the diagonal elements of matrix a as a vector.
Definition mathlib.F:493
mathlib::diag
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition mathlib.F:1496
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition parallel_rng_types.F:51
parallel_rng_types::uniform
integer, parameter, public uniform
Definition parallel_rng_types.F:73
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::a_fine
real(kind=dp), parameter, public a_fine
Definition physcon.F:119
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_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, 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:526
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::centers_spreads_berry
subroutine, public centers_spreads_berry(qs_loc_env, zij, nmoloc, cell, weights, ispin, print_loc_section, only_initial_out)
...
Definition qs_loc_methods.F:578
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:786
qs_loc_types
New version of the module for the localization of the molecular orbitals This should be able to use d...
Definition qs_loc_types.F:25
qs_loc_types::localized_wfn_control_create
subroutine, public localized_wfn_control_create(localized_wfn_control)
create the localized_wfn_control_type
Definition qs_loc_types.F:234
qs_loc_types::qs_loc_env_release
subroutine, public qs_loc_env_release(qs_loc_env)
...
Definition qs_loc_types.F:192
qs_loc_types::get_qs_loc_env
subroutine, public get_qs_loc_env(qs_loc_env, cell, local_molecules, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set, para_env, particle_set, weights, dim_op)
...
Definition qs_loc_types.F:317
qs_loc_types::qs_loc_env_create
subroutine, public qs_loc_env_create(qs_loc_env)
...
Definition qs_loc_types.F:166
qs_loc_utils
Some utilities for the construction of the localization environment.
Definition qs_loc_utils.F:13
qs_loc_utils::qs_loc_env_init
subroutine, public qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, loc_coeff, mo_loc_history)
allocates the data, and initializes the operators
Definition qs_loc_utils.F:232
qs_loc_utils::set_loc_centers
subroutine, public set_loc_centers(localized_wfn_control, nmoloc, nspins)
create the center and spread array and the file names for the output
Definition qs_loc_utils.F:1601
qs_loc_utils::qs_loc_control_init
subroutine, public qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, do_xas, nloc_xas, spin_xas)
initializes everything needed for localization of the HOMOs
Definition qs_loc_utils.F:1006
qs_mo_io
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
qs_mo_io::write_mo_set_low
subroutine, public write_mo_set_low(mo_array, qs_kind_set, particle_set, ires, rt_mos, matrix_ks)
...
Definition qs_mo_io.F:287
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::duplicate_mo_set
subroutine, public duplicate_mo_set(mo_set_new, mo_set_old)
allocate a new mo_set, and copy the old data
Definition qs_mo_types.F:149
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:206
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition qs_mo_types.F:352
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:245
qs_operators_ao
Definition qs_operators_ao.F:13
qs_operators_ao::p_xyz_ao
subroutine, public p_xyz_ao(op, qs_env, minimum_image)
Calculation of the components of the dipole operator in the velocity form The elements of the sparse ...
Definition qs_operators_ao.F:818
qs_operators_ao::rrc_xyz_ao
subroutine, public rrc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
Calculation of the components of the dipole operator in the length form by taking the relative positi...
Definition qs_operators_ao.F:1084
qs_pdos
Calculation and writing of projected density of states The DOS is computed per angular momentum and p...
Definition qs_pdos.F:15
qs_pdos::calculate_projected_dos
subroutine, public calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, dft_section, ispin, xas_mittle, external_matrix_shalf)
Compute and write projected density of states.
Definition qs_pdos.F:139
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
qs_scf_types::ot_method_nr
integer, parameter, public ot_method_nr
Definition qs_scf_types.F:51
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:160
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:3376
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:1878
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:1386
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:1132
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:1227
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:1045
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:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:622
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:113
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:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F: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:121
qs_loc_types::qs_loc_env_type
contains all the info needed by quickstep to calculate the spread of a selected set of orbitals and i...
Definition qs_loc_types.F:80
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
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