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