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