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