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 Resonant Inelastic XRAY Scattering (RIXS) calculations
10 : !> \author BSG (02.2025)
11 : ! **************************************************************************************************
12 : MODULE rixs_methods
13 : USE bibliography, ONLY: VazdaCruz2021,&
14 : cite_reference
15 : USE cp_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_control_types, ONLY: dft_control_type,&
17 : rixs_control_create,&
18 : rixs_control_release,&
19 : rixs_control_type
20 : USE cp_control_utils, ONLY: read_rixs_control
21 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
22 : dbcsr_type
23 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
24 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 : cp_fm_struct_release,&
26 : cp_fm_struct_type
27 : USE cp_fm_types, ONLY: cp_fm_create,&
28 : cp_fm_get_info,&
29 : cp_fm_get_submatrix,&
30 : cp_fm_release,&
31 : cp_fm_to_fm,&
32 : cp_fm_to_fm_submat,&
33 : cp_fm_type
34 : USE cp_log_handling, ONLY: cp_get_default_logger,&
35 : cp_logger_get_default_io_unit,&
36 : cp_logger_type
37 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
38 : cp_print_key_unit_nr
39 : USE header, ONLY: rixs_header
40 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
41 : section_vals_type
42 : USE kinds, ONLY: dp
43 : USE message_passing, ONLY: mp_para_env_type
44 : USE parallel_gemm_api, ONLY: parallel_gemm
45 : USE physcon, ONLY: evolt
46 : USE qs_environment_types, ONLY: get_qs_env,&
47 : qs_environment_type
48 : USE qs_tddfpt2_methods, ONLY: tddfpt
49 : USE rixs_types, ONLY: rixs_env_create,&
50 : rixs_env_release,&
51 : rixs_env_type,&
52 : tddfpt2_valence_type
53 : USE xas_tdp_methods, ONLY: xas_tdp
54 : USE xas_tdp_types, ONLY: donor_state_type,&
55 : xas_tdp_env_type
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 : PRIVATE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rixs_methods'
62 :
63 : PUBLIC :: rixs, rixs_core
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief Driver for RIXS calculations.
69 : !> \param qs_env the inherited qs_environment
70 : !> \author BSG
71 : ! **************************************************************************************************
72 :
73 14 : SUBROUTINE rixs(qs_env)
74 :
75 : TYPE(qs_environment_type), POINTER :: qs_env
76 :
77 : CHARACTER(len=*), PARAMETER :: routineN = 'rixs'
78 :
79 : INTEGER :: handle, output_unit
80 : TYPE(dft_control_type), POINTER :: dft_control
81 : TYPE(section_vals_type), POINTER :: rixs_section, tddfp2_section, &
82 : xas_tdp_section
83 :
84 14 : CALL timeset(routineN, handle)
85 :
86 14 : NULLIFY (rixs_section)
87 14 : rixs_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS")
88 14 : output_unit = cp_logger_get_default_io_unit()
89 :
90 14 : qs_env%do_rixs = .TRUE.
91 :
92 14 : CALL cite_reference(VazdaCruz2021)
93 :
94 14 : CALL get_qs_env(qs_env, dft_control=dft_control)
95 :
96 14 : xas_tdp_section => section_vals_get_subs_vals(rixs_section, "XAS_TDP")
97 14 : tddfp2_section => section_vals_get_subs_vals(rixs_section, "TDDFPT")
98 :
99 14 : CALL rixs_core(rixs_section, qs_env)
100 :
101 14 : IF (output_unit > 0) THEN
102 : WRITE (UNIT=output_unit, FMT="(/,(T2,A79))") &
103 7 : "*******************************************************************************", &
104 7 : "! Normal termination of Resonant Inelastic X-RAY Scattering calculation !", &
105 14 : "*******************************************************************************"
106 : END IF
107 :
108 14 : CALL timestop(handle)
109 :
110 14 : END SUBROUTINE rixs
111 :
112 : ! **************************************************************************************************
113 : !> \brief Perform RIXS calculation.
114 : !> \param rixs_section ...
115 : !> \param qs_env ...
116 : ! **************************************************************************************************
117 14 : SUBROUTINE rixs_core(rixs_section, qs_env)
118 :
119 : TYPE(section_vals_type), POINTER :: rixs_section
120 : TYPE(qs_environment_type), POINTER :: qs_env
121 :
122 : CHARACTER(len=*), PARAMETER :: routineN = 'rixs_core'
123 :
124 : INTEGER :: ax, c_nstates, current_state_index, fstate, handle, iatom, ispin, istate, &
125 : nactive_max, nao, ncol, nex_atoms, nocc, nspins, output_unit, td_state, v_nstates
126 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nactive
127 : LOGICAL :: do_sc, do_sg, roks, uks
128 14 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: w_i0, w_if
129 14 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dip_block, mu_i0
130 14 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: mu_if
131 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
132 : TYPE(cp_fm_struct_type), POINTER :: core_evect_struct, cRc_struct, dip_0_struct, &
133 : dip_f_struct, gs_coeff_struct, i_dip_0_struct, i_dip_f_struct, Rc_struct
134 : TYPE(cp_fm_type) :: dip_0
135 14 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: core_evects, cRc, dip_f, i_dip_0, &
136 14 : i_dip_f, Rc, state_gs_coeffs
137 14 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: local_gs_coeffs, mo_coeffs
138 14 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: valence_evects
139 : TYPE(cp_fm_type), POINTER :: target_ex_coeffs
140 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
141 : TYPE(dft_control_type), POINTER :: dft_control
142 : TYPE(donor_state_type), POINTER :: current_state
143 : TYPE(mp_para_env_type), POINTER :: para_env
144 : TYPE(rixs_control_type), POINTER :: rixs_control
145 : TYPE(rixs_env_type), POINTER :: rixs_env
146 : TYPE(tddfpt2_valence_type), POINTER :: valence_state
147 : TYPE(xas_tdp_env_type), POINTER :: core_state
148 :
149 14 : NULLIFY (rixs_control, dft_control, rixs_env)
150 14 : NULLIFY (valence_state, core_state)
151 14 : NULLIFY (para_env, blacs_env)
152 14 : NULLIFY (local_gs_coeffs, mo_coeffs, valence_evects)
153 14 : NULLIFY (dipmat, dip_0_struct, i_dip_0_struct, dip_f_struct, i_dip_f_struct, &
154 14 : core_evect_struct, gs_coeff_struct)
155 :
156 28 : output_unit = cp_logger_get_default_io_unit()
157 :
158 : CALL get_qs_env(qs_env, &
159 : dft_control=dft_control, &
160 : matrix_s=matrix_s, &
161 : para_env=para_env, &
162 14 : blacs_env=blacs_env)
163 14 : CALL rixs_control_create(rixs_control)
164 14 : CALL read_rixs_control(rixs_control, rixs_section, dft_control%qs_control)
165 :
166 : ! create rixs_env
167 14 : CALL rixs_env_create(rixs_env)
168 :
169 : ! first, xas_tdp calculation
170 14 : CALL xas_tdp(qs_env, rixs_env)
171 :
172 14 : do_sg = rixs_control%xas_tdp_control%do_singlet
173 14 : do_sc = rixs_control%xas_tdp_control%do_spin_cons
174 :
175 14 : IF (rixs_control%xas_tdp_control%check_only) THEN
176 0 : CPWARN("CHECK_ONLY run for XAS_TDP requested, RIXS and TDDFPT will not be performed.")
177 : ELSE
178 :
179 : ! then, tddfpt calculation
180 14 : CALL tddfpt(qs_env, calc_forces=.FALSE., rixs_env=rixs_env)
181 :
182 14 : IF (output_unit > 0) THEN
183 7 : CALL rixs_header(output_unit)
184 : END IF
185 :
186 : ! timings for rixs only, excluding xas_tdp and tddft calls
187 14 : CALL timeset(routineN, handle)
188 :
189 14 : IF (do_sg) THEN ! singlet
190 : nspins = 1
191 4 : ELSE IF (do_sc) THEN ! spin-conserving
192 : nspins = 2
193 : ELSE
194 0 : CPABORT("RIXS only implemented for singlet and spin-conserving excitations")
195 : END IF
196 :
197 14 : IF (output_unit > 0) THEN
198 7 : IF (dft_control%uks) THEN
199 1 : uks = .TRUE.
200 1 : WRITE (UNIT=output_unit, FMT="(T2,A)") "RIXS| Unrestricted Open-Shell Kohn-Sham"
201 6 : ELSE IF (dft_control%roks) THEN
202 1 : roks = .TRUE.
203 1 : WRITE (UNIT=output_unit, FMT="(T2,A)") "RIXS| Restricted Open-Shell Kohn-Sham"
204 : END IF
205 : END IF
206 :
207 14 : core_state => rixs_env%core_state
208 14 : valence_state => rixs_env%valence_state
209 :
210 : ! gs coefficients from tddfpt
211 14 : mo_coeffs => valence_state%mos_active
212 : ! localised gs coefficients from xas_tdp
213 14 : local_gs_coeffs => core_state%mo_coeff
214 14 : valence_evects => valence_state%evects
215 :
216 : ! res tddft
217 : IF (.NOT. roks) THEN
218 : CALL cp_fm_get_info(matrix=local_gs_coeffs(1), ncol_global=nocc)
219 : CALL cp_fm_get_info(matrix=mo_coeffs(1), ncol_global=ncol)
220 : IF (ncol /= nocc) THEN
221 : CPABORT("RIXS with restricted space excitations NYI")
222 : END IF
223 : END IF
224 :
225 14 : IF (rixs_control%xas_tdp_control%do_loc) THEN
226 2 : IF (output_unit > 0) THEN
227 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
228 1 : "RIXS| Found localised XAS_TDP orbitals"
229 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
230 1 : "RIXS| Rotating TDDFPT vectors..."
231 : END IF
232 2 : CALL rotate_vectors(valence_state%evects, local_gs_coeffs, mo_coeffs, matrix_s(1)%matrix, output_unit)
233 : END IF
234 :
235 : ! find max nactive for open-shell cases
236 28 : ALLOCATE (nactive(nspins))
237 32 : DO ispin = 1, nspins
238 32 : CALL cp_fm_get_info(matrix=mo_coeffs(ispin), nrow_global=nao, ncol_global=nactive(ispin))
239 : END DO
240 32 : nactive_max = MAXVAL(nactive)
241 :
242 14 : nex_atoms = core_state%nex_atoms
243 14 : v_nstates = valence_state%nstates
244 14 : c_nstates = core_state%nstates
245 :
246 14 : IF (rixs_control%core_states > 0) THEN
247 4 : rixs_control%core_states = MIN(rixs_control%core_states, c_nstates)
248 : ELSE
249 10 : rixs_control%core_states = c_nstates
250 : END IF
251 14 : c_nstates = rixs_control%core_states
252 :
253 14 : IF (rixs_control%valence_states > 0) THEN
254 2 : rixs_control%valence_states = MIN(rixs_control%valence_states, v_nstates)
255 : ELSE
256 12 : rixs_control%valence_states = v_nstates
257 : END IF
258 14 : v_nstates = rixs_control%valence_states
259 :
260 14 : IF (output_unit > 0) THEN
261 : WRITE (UNIT=output_unit, FMT="(T2,A,I5,A,I5)") &
262 7 : "RIXS| Using ", c_nstates, " core states out of ", core_state%nstates
263 : WRITE (UNIT=output_unit, FMT="(T2,A,I5,A,I5,/)") &
264 7 : "RIXS| Using ", v_nstates, " valence states out of ", valence_state%nstates
265 : END IF
266 :
267 14 : dipmat => core_state%dipmat
268 :
269 78 : ALLOCATE (core_evects(nspins), state_gs_coeffs(nspins))
270 154 : ALLOCATE (dip_block(1, nspins), mu_i0(4, c_nstates), mu_if(4, c_nstates, v_nstates), w_i0(c_nstates), w_if(v_nstates))
271 60 : w_if(:) = valence_state%evals(1:v_nstates)*evolt
272 46 : ALLOCATE (i_dip_0(nspins))
273 78 : ALLOCATE (dip_f(nspins), i_dip_f(nspins))
274 78 : ALLOCATE (Rc(nspins), cRc(nspins))
275 :
276 : CALL cp_fm_struct_create(core_evect_struct, para_env=para_env, context=blacs_env, &
277 14 : nrow_global=nao, ncol_global=c_nstates)
278 : CALL cp_fm_struct_create(gs_coeff_struct, para_env=para_env, context=blacs_env, &
279 14 : nrow_global=nao, ncol_global=1)
280 :
281 : ! looping over ex_atoms and ex_kinds is enough as excited atoms have to be unique
282 14 : current_state_index = 1
283 30 : DO iatom = 1, nex_atoms
284 16 : current_state => core_state%donor_states(current_state_index)
285 16 : IF (output_unit > 0) THEN
286 : WRITE (UNIT=output_unit, FMT="(T2,A,A,A,I3,A,A)") &
287 8 : "RIXS| Calculating dipole moment from core-excited state ", &
288 8 : core_state%state_type_char(current_state%state_type), " for atom ", &
289 16 : current_state%at_index, " of kind ", TRIM(current_state%at_symbol)
290 : END IF
291 :
292 846 : mu_i0 = 0.0_dp
293 2786 : mu_if = 0.0_dp
294 :
295 16 : IF (do_sg) THEN ! singlet
296 12 : target_ex_coeffs => current_state%sg_coeffs
297 130 : w_i0(:) = current_state%sg_evals(1:c_nstates)*evolt
298 4 : ELSE IF (do_sc) THEN ! spin-conserving
299 4 : target_ex_coeffs => current_state%sc_coeffs
300 52 : w_i0(:) = current_state%sc_evals(1:c_nstates)*evolt
301 : END IF
302 :
303 : ! reshape sc and sg coeffs (separate spins to columns)
304 36 : DO ispin = 1, nspins
305 20 : CALL cp_fm_create(core_evects(ispin), core_evect_struct)
306 : CALL cp_fm_to_fm_submat(msource=target_ex_coeffs, mtarget=core_evects(ispin), s_firstrow=1, &
307 36 : s_firstcol=(c_nstates*(ispin - 1) + 1), t_firstrow=1, t_firstcol=1, nrow=nao, ncol=c_nstates)
308 : END DO
309 36 : DO ispin = 1, nspins
310 20 : CALL cp_fm_create(state_gs_coeffs(ispin), gs_coeff_struct)
311 16 : IF (roks) THEN
312 : ! store same coeffs for both spins, easier later on
313 : CALL cp_fm_to_fm_submat(msource=current_state%gs_coeffs, mtarget=state_gs_coeffs(ispin), s_firstrow=1, &
314 20 : s_firstcol=1, t_firstrow=1, t_firstcol=1, nrow=nao, ncol=1)
315 : ELSE
316 : CALL cp_fm_to_fm_submat(msource=current_state%gs_coeffs, mtarget=state_gs_coeffs(ispin), s_firstrow=1, &
317 : s_firstcol=ispin, t_firstrow=1, t_firstcol=1, nrow=nao, ncol=1)
318 : END IF
319 : END DO
320 :
321 : ! initialise matrices for i->0
322 : CALL cp_fm_struct_create(dip_0_struct, para_env=para_env, context=blacs_env, &
323 16 : nrow_global=nao, ncol_global=1)
324 16 : CALL cp_fm_create(dip_0, dip_0_struct)
325 : CALL cp_fm_struct_create(i_dip_0_struct, para_env=para_env, context=blacs_env, &
326 16 : nrow_global=c_nstates, ncol_global=1)
327 36 : DO ispin = 1, nspins
328 36 : CALL cp_fm_create(i_dip_0(ispin), i_dip_0_struct)
329 : END DO
330 :
331 : ! initialise matrices for i->f
332 36 : DO ispin = 1, nspins
333 : CALL cp_fm_struct_create(dip_f_struct, para_env=para_env, context=blacs_env, &
334 20 : nrow_global=nao, ncol_global=nactive(ispin))
335 : CALL cp_fm_struct_create(i_dip_f_struct, para_env=para_env, context=blacs_env, &
336 20 : nrow_global=c_nstates, ncol_global=nactive(ispin))
337 20 : CALL cp_fm_create(dip_f(ispin), dip_f_struct)
338 20 : CALL cp_fm_create(i_dip_f(ispin), i_dip_f_struct)
339 20 : CALL cp_fm_struct_release(i_dip_f_struct)
340 36 : CALL cp_fm_struct_release(dip_f_struct)
341 : END DO
342 :
343 : ! initialise cRc
344 36 : DO ispin = 1, nspins
345 : CALL cp_fm_struct_create(Rc_struct, para_env=para_env, context=blacs_env, &
346 20 : nrow_global=nao, ncol_global=nactive(ispin))
347 : CALL cp_fm_struct_create(cRc_struct, para_env=para_env, context=blacs_env, &
348 20 : nrow_global=nactive(ispin), ncol_global=nactive(ispin))
349 20 : CALL cp_fm_create(Rc(ispin), Rc_struct)
350 20 : CALL cp_fm_create(cRc(ispin), cRc_struct)
351 20 : CALL cp_fm_struct_release(cRc_struct)
352 36 : CALL cp_fm_struct_release(Rc_struct)
353 : END DO
354 :
355 : ! 0 -> i
356 64 : DO ax = 1, 3
357 :
358 : ! i*R*0
359 108 : DO ispin = 1, nspins
360 60 : CALL cp_dbcsr_sm_fm_multiply(dipmat(ax)%matrix, state_gs_coeffs(ispin), dip_0, ncol=1)
361 108 : CALL parallel_gemm('T', 'N', c_nstates, 1, nao, 1.0_dp, core_evects(ispin), dip_0, 0.0_dp, i_dip_0(ispin))
362 : END DO
363 :
364 562 : DO istate = 1, c_nstates
365 1782 : dip_block = 0.0_dp
366 1140 : DO ispin = 1, nspins
367 : CALL cp_fm_get_submatrix(fm=i_dip_0(ispin), target_m=dip_block, start_row=istate, &
368 642 : start_col=1, n_rows=1, n_cols=1)
369 1140 : mu_i0(ax, istate) = mu_i0(ax, istate) + dip_block(1, 1)
370 : END DO ! ispin
371 546 : mu_i0(4, istate) = mu_i0(4, istate) + mu_i0(ax, istate)**2
372 : END DO ! istate
373 :
374 : END DO ! ax
375 :
376 : ! i -> f
377 66 : DO td_state = 1, v_nstates
378 :
379 50 : IF (output_unit > 0) THEN
380 : WRITE (UNIT=output_unit, FMT="(T9,A,I3,A,F10.4)") &
381 25 : "to valence-excited state ", td_state, " with energy ", w_if(td_state)
382 : END IF
383 :
384 216 : DO ax = 1, 3
385 :
386 360 : DO ispin = 1, nspins
387 :
388 : ! 1) build c^T*R*c
389 210 : CALL cp_dbcsr_sm_fm_multiply(dipmat(ax)%matrix, mo_coeffs(ispin), Rc(ispin), ncol=nactive(ispin))
390 : CALL parallel_gemm('T', 'N', nactive(ispin), nactive(ispin), nao, 1.0_dp, mo_coeffs(ispin), &
391 210 : Rc(ispin), 0.0_dp, cRc(ispin))
392 :
393 : ! 2) build dip_f (X*dip = X*c^T*R*c)
394 : CALL parallel_gemm('N', 'N', nao, nactive(ispin), nactive(ispin), 1.0_dp, valence_evects(ispin, td_state), &
395 210 : cRc(ispin), 0.0_dp, dip_f(ispin))
396 :
397 : ! 3) build i_dip_f (Y*dip_f = Y^T*X*c^T*R*c)
398 : CALL parallel_gemm('T', 'N', c_nstates, nactive(ispin), nao, 1.0_dp, core_evects(ispin), &
399 360 : dip_f(ispin), 0.0_dp, i_dip_f(ispin))
400 :
401 : END DO
402 :
403 1832 : DO istate = 1, c_nstates
404 9696 : DO fstate = 1, nactive_max
405 31392 : dip_block = 0.0_dp
406 21360 : DO ispin = 1, nspins
407 19728 : IF (fstate <= nactive(ispin)) THEN
408 : CALL cp_fm_get_submatrix(fm=i_dip_f(ispin), target_m=dip_block, start_row=istate, &
409 10944 : start_col=fstate, n_rows=1, n_cols=1)
410 10944 : mu_if(ax, istate, td_state) = mu_if(ax, istate, td_state) + dip_block(1, 1)
411 : END IF
412 : END DO ! ispin
413 : END DO ! fstate (tddft)
414 1782 : mu_if(4, istate, td_state) = mu_if(4, istate, td_state) + mu_if(ax, istate, td_state)**2
415 : END DO ! istate (core)
416 :
417 : END DO ! ax
418 :
419 : END DO ! td_state
420 :
421 16 : IF (output_unit > 0) THEN
422 8 : WRITE (UNIT=output_unit, FMT="(/,T2,A,/)") "RIXS| Printing spectrum to file"
423 : END IF
424 16 : CALL print_rixs_to_file(current_state, mu_i0, mu_if, w_i0, w_if, rixs_env, rixs_section, rixs_control)
425 :
426 16 : current_state_index = current_state_index + 1
427 :
428 : ! cleanup
429 36 : DO ispin = 1, nspins
430 20 : CALL cp_fm_release(core_evects(ispin))
431 20 : CALL cp_fm_release(state_gs_coeffs(ispin))
432 20 : CALL cp_fm_release(i_dip_0(ispin))
433 20 : CALL cp_fm_release(i_dip_f(ispin))
434 20 : CALL cp_fm_release(dip_f(ispin))
435 20 : CALL cp_fm_release(Rc(ispin))
436 36 : CALL cp_fm_release(cRc(ispin))
437 : END DO
438 16 : CALL cp_fm_struct_release(i_dip_0_struct)
439 16 : CALL cp_fm_struct_release(dip_0_struct)
440 46 : CALL cp_fm_release(dip_0)
441 :
442 : END DO ! iatom
443 :
444 : NULLIFY (current_state)
445 :
446 : ! cleanup
447 14 : CALL cp_fm_struct_release(core_evect_struct)
448 28 : CALL cp_fm_struct_release(gs_coeff_struct)
449 :
450 : END IF
451 :
452 : ! more cleanup
453 14 : CALL rixs_control_release(rixs_control)
454 14 : CALL rixs_env_release(rixs_env)
455 14 : NULLIFY (valence_state, core_state)
456 :
457 14 : CALL timestop(handle)
458 :
459 28 : END SUBROUTINE rixs_core
460 :
461 : ! **************************************************************************************************
462 : !> \brief Rotate vectors. Returns rotated mo_occ and evects.
463 : !> \param evects ...
464 : !> \param mo_ref ...
465 : !> \param mo_occ ...
466 : !> \param overlap_matrix ...
467 : !> \param unit_nr ...
468 : ! **************************************************************************************************
469 :
470 2 : SUBROUTINE rotate_vectors(evects, mo_ref, mo_occ, overlap_matrix, unit_nr)
471 : TYPE(cp_fm_type), DIMENSION(:, :) :: evects
472 : TYPE(cp_fm_type), DIMENSION(:) :: mo_ref, mo_occ
473 : TYPE(dbcsr_type), POINTER :: overlap_matrix
474 : INTEGER :: unit_nr
475 :
476 : INTEGER :: ispin, istate, ncol, nrow, nspins, &
477 : nstates
478 : REAL(kind=dp) :: diff
479 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
480 : TYPE(cp_fm_struct_type), POINTER :: emat_struct
481 : TYPE(cp_fm_type) :: emat, rotated_mo_coeffs, smo
482 : TYPE(cp_fm_type), POINTER :: current_evect
483 : TYPE(mp_para_env_type), POINTER :: para_env
484 :
485 2 : NULLIFY (emat_struct, para_env, blacs_env, current_evect)
486 :
487 2 : nspins = SIZE(evects, DIM=1)
488 4 : DO ispin = 1, nspins
489 :
490 : CALL cp_fm_get_info(matrix=mo_occ(ispin), nrow_global=nrow, ncol_global=ncol, &
491 2 : para_env=para_env, context=blacs_env)
492 2 : CALL cp_fm_create(smo, mo_occ(ispin)%matrix_struct)
493 :
494 : ! rotate mo_occ
495 : ! smo = matrix_s x mo_occ
496 2 : CALL cp_dbcsr_sm_fm_multiply(overlap_matrix, mo_occ(ispin), smo, ncol, alpha=1.0_dp, beta=0.0_dp)
497 : CALL cp_fm_struct_create(emat_struct, nrow_global=ncol, ncol_global=ncol, &
498 2 : para_env=para_env, context=blacs_env)
499 2 : CALL cp_fm_create(emat, emat_struct)
500 : ! emat = mo_ref^T x smo
501 2 : CALL parallel_gemm('T', 'N', ncol, ncol, nrow, 1.0_dp, mo_ref(ispin), smo, 0.0_dp, emat)
502 2 : CALL cp_fm_create(rotated_mo_coeffs, mo_occ(ispin)%matrix_struct)
503 : ! rotated_mo_coeffs = cpmos x emat
504 2 : CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, mo_occ(ispin), emat, 0.0_dp, rotated_mo_coeffs)
505 :
506 77 : diff = MAXVAL(ABS(rotated_mo_coeffs%local_data - mo_occ(ispin)%local_data))
507 2 : IF (unit_nr > 0) THEN
508 1 : WRITE (unit_nr, FMT="(T9,A,I2,A,F10.6,/)") "For spin ", ispin, ": Max difference between orbitals = ", diff
509 : END IF
510 :
511 2 : CALL cp_fm_to_fm(rotated_mo_coeffs, mo_occ(ispin))
512 :
513 2 : nstates = SIZE(evects, DIM=2)
514 8 : DO istate = 1, nstates
515 2 : ASSOCIATE (current_evect => evects(ispin, istate))
516 6 : CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, current_evect, emat, 0.0_dp, smo)
517 6 : CALL cp_fm_to_fm(smo, current_evect)
518 : END ASSOCIATE
519 : END DO
520 :
521 2 : CALL cp_fm_struct_release(emat_struct)
522 2 : CALL cp_fm_release(smo)
523 2 : CALL cp_fm_release(emat)
524 10 : CALL cp_fm_release(rotated_mo_coeffs)
525 :
526 : END DO ! ispin
527 :
528 2 : END SUBROUTINE rotate_vectors
529 :
530 : !**************************************************************************************************
531 : !> \brief Print RIXS spectrum.
532 : !> \param donor_state ...
533 : !> \param mu_i0 ...
534 : !> \param mu_if ...
535 : !> \param w_i0 ...
536 : !> \param w_if ...
537 : !> \param rixs_env ...
538 : !> \param rixs_section ...
539 : !> \param rixs_control ...
540 : ! **************************************************************************************************
541 16 : SUBROUTINE print_rixs_to_file(donor_state, mu_i0, mu_if, w_i0, w_if, &
542 : rixs_env, rixs_section, rixs_control)
543 :
544 : TYPE(donor_state_type), POINTER :: donor_state
545 : REAL(dp), DIMENSION(:, :) :: mu_i0
546 : REAL(dp), DIMENSION(:, :, :) :: mu_if
547 : REAL(dp), DIMENSION(:) :: w_i0, w_if
548 : TYPE(rixs_env_type), POINTER :: rixs_env
549 : TYPE(section_vals_type), POINTER :: rixs_section
550 : TYPE(rixs_control_type), POINTER :: rixs_control
551 :
552 : INTEGER :: f, i, output_unit, rixs_unit
553 : TYPE(cp_logger_type), POINTER :: logger
554 :
555 16 : NULLIFY (logger)
556 16 : logger => cp_get_default_logger()
557 :
558 : rixs_unit = cp_print_key_unit_nr(logger, rixs_section, "PRINT%SPECTRUM", &
559 : extension=".rixs", file_position="APPEND", &
560 16 : file_action="WRITE", file_form="FORMATTED")
561 :
562 16 : output_unit = cp_logger_get_default_io_unit()
563 :
564 16 : IF (rixs_unit > 0) THEN
565 :
566 : WRITE (rixs_unit, FMT="(A,/,T2,A,A,A,I3,A,A,A/,A)") &
567 8 : "=====================================================================================", &
568 8 : "Excitation from ground-state (", &
569 8 : rixs_env%core_state%state_type_char(donor_state%state_type), " for atom ", &
570 8 : donor_state%at_index, " of kind ", TRIM(donor_state%at_symbol), &
571 8 : ") to core-excited state i ", &
572 16 : "====================================================================================="
573 :
574 : WRITE (rixs_unit, FMT="(T3,A)") &
575 8 : "w_0i (eV) mu^x_0i (a.u.) mu^y_0i (a.u.) mu^z_0i (a.u.) mu^2_0i (a.u.)"
576 91 : DO i = 1, rixs_control%core_states
577 : WRITE (rixs_unit, FMT="(T2,F10.4,T26,E12.5,T42,E12.5,T58,E12.5,T74,E12.5)") &
578 91 : w_i0(i), mu_i0(1, i), mu_i0(2, i), mu_i0(3, i), mu_i0(4, i)
579 : END DO
580 :
581 : WRITE (rixs_unit, FMT="(A,/,T2,A,/,A)") &
582 8 : "=====================================================================================", &
583 8 : "Emission from core-excited state i to valence-excited state f ", &
584 16 : "====================================================================================="
585 :
586 : WRITE (rixs_unit, FMT="(T3,A)") &
587 8 : "w_0i (eV) w_if (eV) mu^x_if (a.u.) mu^y_if (a.u.) mu^z_if (a.u.) mu^2_if (a.u.)"
588 :
589 91 : DO i = 1, rixs_control%core_states
590 363 : DO f = 1, rixs_control%valence_states
591 : WRITE (rixs_unit, FMT="(T2,F10.4,T14,F8.4,T26,E12.5,T42,E12.5,T58,E12.5,T74,E12.5)") &
592 355 : w_i0(i), w_if(f), mu_if(1, i, f), mu_if(2, i, f), mu_if(3, i, f), mu_if(4, i, f)
593 : END DO
594 : END DO
595 :
596 : END IF
597 :
598 16 : CALL cp_print_key_finished_output(rixs_unit, logger, rixs_section, "PRINT%SPECTRUM")
599 :
600 16 : END SUBROUTINE print_rixs_to_file
601 :
602 : END MODULE rixs_methods
|
