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