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