(git:419edc0)
Loading...
Searching...
No Matches
rt_bse_io.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Input/output from the propagation via RT-BSE method.
10!> \author Stepan Marek (08.24)
11! **************************************************************************************************
12
13MODULE rt_bse_io
14 USE cp_fm_types, ONLY: cp_fm_type, &
15 cp_fm_p_type, &
16 cp_fm_create, &
17 cp_fm_release, &
18 cp_fm_read_unformatted, &
19 cp_fm_write_unformatted, &
20 cp_fm_write_formatted
21 USE cp_cfm_types, ONLY: cp_cfm_type, &
22 cp_cfm_p_type, &
23 cp_fm_to_cfm, &
24 cp_cfm_to_fm
25 USE kinds, ONLY: dp, &
26 default_path_length
27 USE cp_fm_basic_linalg, ONLY: cp_fm_trace, &
28 cp_fm_transpose, &
29 cp_fm_norm
30 USE parallel_gemm_api, ONLY: parallel_gemm
31 USE cp_log_handling, ONLY: cp_logger_type, &
32 cp_get_default_logger
33 USE cp_output_handling, ONLY: cp_print_key_unit_nr, &
34 cp_print_key_finished_output, &
35 cp_print_key_generate_filename, &
36 low_print_level, &
37 medium_print_level
38 USE input_section_types, ONLY: section_vals_type
39 USE rt_bse_types, ONLY: rtbse_env_type, &
40 multiply_cfm_fm, &
41 multiply_fm_cfm
42 USE cp_files, ONLY: open_file, &
43 file_exists, &
44 close_file
45 USE input_constants, ONLY: do_exact, &
46 do_bch, &
47 rtp_bse_ham_g0w0, &
48 rtp_bse_ham_ks, &
49 use_rt_restart
50 USE physcon, ONLY: femtoseconds, &
51 evolt
52 USE mathconstants, ONLY: twopi
53#if defined (__GREENX)
54 USE gx_ac, ONLY: create_thiele_pade, &
55 evaluate_thiele_pade_at, &
56 free_params, &
57 params
58#endif
59
60#include "../base/base_uses.f90"
61
62 IMPLICIT NONE
63
64 PRIVATE
65
66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse_io"
67
68
69
70
71 PUBLIC :: output_moments, &
72 read_moments, &
73 output_moments_ft, &
74 output_polarizability, &
75 output_field, &
76 read_field, &
77 output_mos_contravariant, &
78 output_mos_covariant, &
79 output_restart, &
80 read_restart, &
81 print_etrs_info_header, &
82 print_etrs_info, &
83 print_timestep_info, &
84 print_rtbse_header_info
85
86CONTAINS
87
88! **************************************************************************************************
89!> \brief Writes the header and basic info to the standard output
90!> \param rtbse_env Entry point - rtbse environment
91! **************************************************************************************************
92 SUBROUTINE print_rtbse_header_info(rtbse_env)
93 TYPE(rtbse_env_type) :: rtbse_env
94
95 WRITE (rtbse_env%unit_nr, *) ''
96 WRITE (rtbse_env%unit_nr, '(A)') ' /-----------------------------------------------'// &
97 '------------------------------\'
98 WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
99 ' |'
100 WRITE (rtbse_env%unit_nr, '(A)') ' | Real Time Bethe-Salpeter Propagation'// &
101 ' |'
102 WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
103 ' |'
104 WRITE (rtbse_env%unit_nr, '(A)') ' \-----------------------------------------------'// &
105 '------------------------------/'
106 WRITE (rtbse_env%unit_nr, *) ''
107
108 ! Methods used
109 WRITE (rtbse_env%unit_nr, '(A19)', advance="no") ' Exponential method'
110 SELECT CASE (rtbse_env%mat_exp_method)
111 CASE (do_bch)
112 WRITE (rtbse_env%unit_nr, '(A61)') 'BCH'
113 CASE (do_exact)
114 WRITE (rtbse_env%unit_nr, '(A61)') 'EXACT'
115 END SELECT
116
117 WRITE (rtbse_env%unit_nr, '(A22)', advance="no") ' Reference Hamiltonian'
118 SELECT CASE (rtbse_env%ham_reference_type)
119 CASE (rtp_bse_ham_g0w0)
120 WRITE (rtbse_env%unit_nr, '(A58)') 'G0W0'
121 CASE (rtp_bse_ham_ks)
122 WRITE (rtbse_env%unit_nr, '(A58)') 'Kohn-Sham'
123 END SELECT
124
125 WRITE (rtbse_env%unit_nr, '(A18,L62)') ' Apply delta pulse', &
126 rtbse_env%dft_control%rtp_control%apply_delta_pulse
127
128 WRITE (rtbse_env%unit_nr, '(A)') ''
129
130 END SUBROUTINE
131
132! **************************************************************************************************
133!> \brief Writes the update after single etrs iteration - only for log level > medium
134!> \param rtbse_env Entry point - rtbse environment
135! **************************************************************************************************
136 SUBROUTINE print_etrs_info(rtbse_env, step, metric)
137 TYPE(rtbse_env_type) :: rtbse_env
138 INTEGER :: step
139 REAL(kind=dp) :: metric
140 TYPE(cp_logger_type), POINTER :: logger
141
142 logger => cp_get_default_logger()
143
144 IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN
145 WRITE (rtbse_env%unit_nr, '(A7,I5, E20.8E3)') ' RTBSE|', step, metric
146 END IF
147
148 END SUBROUTINE
149! **************************************************************************************************
150!> \brief Writes the header for the etrs iteration updates - only for log level > medium
151!> \param rtbse_env Entry point - rtbse environment
152! **************************************************************************************************
153 SUBROUTINE print_etrs_info_header(rtbse_env)
154 TYPE(rtbse_env_type) :: rtbse_env
155 TYPE(cp_logger_type), POINTER :: logger
156
157 logger => cp_get_default_logger()
158
159 IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN
160 WRITE (rtbse_env%unit_nr, '(A13, A20)') ' RTBSE| Iter.', 'Convergence'
161 END IF
162
163 END SUBROUTINE
164! **************************************************************************************************
165!> \brief Writes the header for the etrs iteration updates - only for log level > medium
166!> \param rtbse_env Entry point - rtbse environment
167! **************************************************************************************************
168 SUBROUTINE print_timestep_info(rtbse_env, step, convergence, electron_num_re, etrs_num)
169 TYPE(rtbse_env_type) :: rtbse_env
170 INTEGER :: step
171 REAL(kind=dp) :: convergence
172 REAL(kind=dp) :: electron_num_re
173 INTEGER :: etrs_num
174 TYPE(cp_logger_type), POINTER :: logger
175
176 logger => cp_get_default_logger()
177
178 IF (logger%iter_info%print_level > low_print_level .AND. rtbse_env%unit_nr > 0) THEN
179 WRITE (rtbse_env%unit_nr, '(A23,A20,A20,A17)') " RTBSE| Simulation step", "Convergence", &
180 "Electron number", "ETRS Iterations"
181 WRITE (rtbse_env%unit_nr, '(A7,I16,E20.8E3,E20.8E3,I17)') ' RTBSE|', step, convergence, &
182 electron_num_re, etrs_num
183 END IF
184
185 END SUBROUTINE
186
187! **************************************************************************************************
188!> \brief Outputs the matrix in MO basis for matrix coefficients corresponding to contravariant
189!> operator, i.e. density matrix
190!> \param rtbse_env Entry point - gwbse environment
191!> \param rho Density matrix in AO basis
192!> \param rtp_section RTP input section
193! **************************************************************************************************
194 SUBROUTINE output_mos_contravariant(rtbse_env, rho, print_key_section)
195 TYPE(rtbse_env_type) :: rtbse_env
196 TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
197 TYPE(section_vals_type), POINTER :: print_key_section
198 TYPE(cp_logger_type), POINTER :: logger
199 INTEGER :: j, rho_unit_re, rho_unit_im
200 CHARACTER(len=14), DIMENSION(4) :: file_labels
201
202 file_labels(1) = "_SPIN_A_RE.dat"
203 file_labels(2) = "_SPIN_A_IM.dat"
204 file_labels(3) = "_SPIN_B_RE.dat"
205 file_labels(4) = "_SPIN_B_IM.dat"
206 logger => cp_get_default_logger()
207 ! Start by multiplying the current density by MOS
208 DO j = 1, rtbse_env%n_spin
209 rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1))
210 rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j))
211 ! Transform the density matrix into molecular orbitals basis and print it out
212 ! S * rho
213 CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
214 1.0_dp, rtbse_env%S_fm, rho(j), &
215 0.0_dp, rtbse_env%rho_workspace(1))
216 ! C^T * S * rho
217 CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
218 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), rtbse_env%rho_workspace(1), &
219 0.0_dp, rtbse_env%rho_workspace(2))
220 ! C^T * S * rho * S
221 CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
222 1.0_dp, rtbse_env%rho_workspace(2), rtbse_env%S_fm, &
223 0.0_dp, rtbse_env%rho_workspace(1))
224 ! C^T * S * rho * S * C
225 CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
226 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), &
227 0.0_dp, rtbse_env%rho_workspace(2))
228 ! Print real and imaginary parts separately
229 CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), &
230 rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
231 CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re)
232 CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im)
233 CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section)
234 CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section)
235 END DO
236 END SUBROUTINE
237! **************************************************************************************************
238!> \brief Outputs the matrix in MO basis for matrix components corresponding to covariant representation,
239!> i.e. the Hamiltonian matrix
240!> \param rtbse_env Entry point - gwbse environment
241!> \param cohsex cohsex matrix in AO basis, covariant representation
242!> \param rtp_section RTP input section
243! **************************************************************************************************
244 SUBROUTINE output_mos_covariant(rtbse_env, ham, print_key_section)
245 TYPE(rtbse_env_type) :: rtbse_env
246 TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham
247 TYPE(section_vals_type), POINTER :: print_key_section
248 TYPE(cp_logger_type), POINTER :: logger
249 INTEGER :: j, rho_unit_re, rho_unit_im
250 CHARACTER(len=21), DIMENSION(4) :: file_labels
251
252 file_labels(1) = "_SPIN_A_RE.dat"
253 file_labels(2) = "_SPIN_A_IM.dat"
254 file_labels(3) = "_SPIN_B_RE.dat"
255 file_labels(4) = "_SPIN_B_IM.dat"
256 logger => cp_get_default_logger()
257 DO j = 1, rtbse_env%n_spin
258 rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1))
259 rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j))
260 ! C^T * cohsex
261 CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
262 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), ham(j), &
263 0.0_dp, rtbse_env%rho_workspace(1))
264 ! C^T * cohsex * C
265 CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
266 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), &
267 0.0_dp, rtbse_env%rho_workspace(2))
268 ! Print real and imaginary parts separately
269 CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), &
270 rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
271 CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re)
272 CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im)
273 CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section)
274 CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section)
275 END DO
276 END SUBROUTINE
277! **************************************************************************************************
278!> \brief Prints the current field components into a file provided by input
279!> \param rtbse_env Entry point - gwbse environment
280!> \param rtp_section RTP input section
281! **************************************************************************************************
282 SUBROUTINE output_field(rtbse_env)
283 TYPE(rtbse_env_type) :: rtbse_env
284 TYPE(cp_logger_type), POINTER :: logger
285 INTEGER :: field_unit, n, i
286
287 ! Get logger
288 logger => cp_get_default_logger()
289 ! Get file descriptor
290 field_unit = cp_print_key_unit_nr(logger, rtbse_env%field_section, extension=".dat")
291 ! If the file descriptor is non-zero, output field
292 ! TODO : Output also in SI
293 IF (field_unit /= -1) THEN
294 WRITE (field_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') rtbse_env%sim_time*femtoseconds, &
295 rtbse_env%field(1), rtbse_env%field(2), rtbse_env%field(3)
296 END IF
297 ! Write the output to memory for FT
298 ! Need the absolute index
299 n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1
300 DO i = 1, 3
301 rtbse_env%field_trace(i)%series(n) = rtbse_env%field(i)
302 END DO
303 rtbse_env%time_trace%series(n) = rtbse_env%sim_time
304 CALL cp_print_key_finished_output(field_unit, logger, rtbse_env%field_section)
305
306 END SUBROUTINE
307! **************************************************************************************************
308!> \brief Reads the field from the files provided by input - useful for the continuation run
309!> \param rtbse_env Entry point - gwbse environment
310!> \param rtp_section RTP input section
311! **************************************************************************************************
312 SUBROUTINE read_field(rtbse_env)
313 TYPE(rtbse_env_type) :: rtbse_env
314 TYPE(cp_logger_type), POINTER :: logger
315 CHARACTER(len=default_path_length) :: save_name
316 INTEGER :: k, n, field_unit
317
318 ! Get logger
319 logger => cp_get_default_logger()
320 ! Get file name
321 save_name = cp_print_key_generate_filename(logger, rtbse_env%field_section, extension=".dat", my_local=.false.)
322 IF (file_exists(save_name)) THEN
323 CALL open_file(save_name, file_status="OLD", file_form="FORMATTED", file_action="READ", &
324 unit_number=field_unit)
325 DO k = rtbse_env%sim_start_orig, rtbse_env%sim_start
326 n = k - rtbse_env%sim_start_orig + 1
327 READ (field_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') rtbse_env%time_trace%series(n), &
328 rtbse_env%field_trace(1)%series(n), rtbse_env%field_trace(2)%series(n), rtbse_env%field_trace(3)%series(n)
329 ! Set the time units back to atomic units
330 rtbse_env%time_trace%series(n) = rtbse_env%time_trace%series(n)/femtoseconds
331 END DO
332 CALL close_file(field_unit)
333 ELSE IF (.NOT. rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. &
334 rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
335 cpwarn("Restart without RT field file - unknown field trace set to zero.")
336 END IF
337 END SUBROUTINE read_field
338
339! **************************************************************************************************
340!> \brief Outputs the expectation value of moments from a given density matrix
341!> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3)
342!> \param rtbse_env Entry point - gwbse environment
343!> \param rho Density matrix in AO basis
344!> \param rtp_section RTP section of the input parameters, where moments destination may be present
345! **************************************************************************************************
346 SUBROUTINE output_moments(rtbse_env, rho)
347 TYPE(rtbse_env_type) :: rtbse_env
348 TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
349 TYPE(cp_logger_type), POINTER :: logger
350 INTEGER :: i, j, n, moments_unit_re, moments_unit_im
351 CHARACTER(len=14), DIMENSION(4) :: file_labels
352 REAL(kind=dp), DIMENSION(3) :: moments
353
354 ! Start by getting the relevant file unit
355 file_labels(1) = "_SPIN_A_RE.dat"
356 file_labels(2) = "_SPIN_A_IM.dat"
357 file_labels(3) = "_SPIN_B_RE.dat"
358 file_labels(4) = "_SPIN_B_IM.dat"
359 logger => cp_get_default_logger()
360 DO j = 1, rtbse_env%n_spin
361 moments_unit_re = cp_print_key_unit_nr(logger, rtbse_env%moments_section, extension=file_labels(2*j - 1))
362 moments_unit_im = cp_print_key_unit_nr(logger, rtbse_env%moments_section, extension=file_labels(2*j))
363 ! If, for any reason, the file unit is not provided, skip to next cycle immediately
364 ! TODO : Specify output units in config
365 ! Need to transpose due to the definition of trace function
366 CALL cp_cfm_to_fm(msource=rho(j), mtargetr=rtbse_env%real_workspace(2))
367 DO i = 1, 3
368 ! Moments should be symmetric, test without transopose?
369 CALL cp_fm_transpose(rtbse_env%moments(i), rtbse_env%real_workspace(1))
370 CALL cp_fm_trace(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), moments(i))
371 ! Scale by spin degeneracy and electron charge
372 moments(i) = -moments(i)*rtbse_env%spin_degeneracy
373 END DO
374 ! Output to the file
375 IF (moments_unit_re > 0) THEN
376 ! If outputting to standard output, also prepend identifying 'MOMENTS_TRACE|' string
377 IF (rtbse_env%unit_nr == moments_unit_re) THEN
378 WRITE (moments_unit_re, '(A18,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
379 " MOMENTS_TRACE_RE|", rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
380 ELSE
381 WRITE (moments_unit_re, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
382 rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
383 CALL cp_print_key_finished_output(moments_unit_re, logger, rtbse_env%moments_section)
384 END IF
385 END IF
386 ! Save to memory for FT - real part
387 n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1
388 DO i = 1, 3
389 rtbse_env%moments_trace(i)%series(n) = cmplx(moments(i), 0.0, kind=dp)
390 END DO
391 ! Same for imaginary part
392 CALL cp_cfm_to_fm(msource=rho(j), mtargeti=rtbse_env%real_workspace(2))
393 DO i = 1, 3
394 CALL cp_fm_transpose(rtbse_env%moments(i), rtbse_env%real_workspace(1))
395 CALL cp_fm_trace(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), moments(i))
396 ! Scale by spin degeneracy and electron charge
397 moments(i) = -moments(i)*rtbse_env%spin_degeneracy
398 END DO
399 ! Output to the file
400 IF (moments_unit_im > 0) THEN
401 ! If outputting to standard output, also prepend identifying 'MOMENTS_TRACE|' string
402 IF (rtbse_env%unit_nr == moments_unit_im) THEN
403 WRITE (moments_unit_im, '(A18,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
404 " MOMENTS_TRACE_IM|", rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
405 ELSE
406 WRITE (moments_unit_im, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
407 rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
408 ! Close the files
409 CALL cp_print_key_finished_output(moments_unit_im, logger, rtbse_env%moments_section)
410 END IF
411 END IF
412 ! Save to memory for FT - imag part
413 DO i = 1, 3
414 rtbse_env%moments_trace(i)%series(n) = rtbse_env%moments_trace(i)%series(n) + cmplx(0.0, moments(i), kind=dp)
415 END DO
416 END DO
417 END SUBROUTINE
418! **************************************************************************************************
419!> \brief Outputs the restart info (last finished iteration step) + restard density matrix
420!> \param restart_section Print key section for the restart files
421!> \param rho Density matrix in AO basis
422!> \param time_index Time index to be written into the info file
423! **************************************************************************************************
424 SUBROUTINE output_restart(rtbse_env, rho, time_index)
425 TYPE(rtbse_env_type), POINTER :: rtbse_env
426 TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
427 INTEGER :: time_index
428 TYPE(cp_fm_type), DIMENSION(:), POINTER :: workspace
429 CHARACTER(len=17), DIMENSION(4) :: file_labels
430 TYPE(cp_logger_type), POINTER :: logger
431 INTEGER :: rho_unit_nr, i
432
433 ! Default labels distinguishing up to two spin species and real/imaginary parts
434 file_labels(1) = "_SPIN_A_RE.matrix"
435 file_labels(2) = "_SPIN_A_IM.matrix"
436 file_labels(3) = "_SPIN_B_RE.matrix"
437 file_labels(4) = "_SPIN_B_IM.matrix"
438
439 logger => cp_get_default_logger()
440
441 workspace => rtbse_env%real_workspace
442
443 DO i = 1, rtbse_env%n_spin
444 CALL cp_cfm_to_fm(rho(i), workspace(1), workspace(2))
445 ! Real part
446 rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i - 1), &
447 file_form="UNFORMATTED", file_position="REWIND")
448 CALL cp_fm_write_unformatted(workspace(1), rho_unit_nr)
449 CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
450 ! Imag part
451 rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i), &
452 file_form="UNFORMATTED", file_position="REWIND")
453 CALL cp_fm_write_unformatted(workspace(2), rho_unit_nr)
454 CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
455 ! Info
456 rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=".info", &
457 file_form="UNFORMATTED", file_position="REWIND")
458 IF (rho_unit_nr > 0) WRITE (rho_unit_nr) time_index
459 CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
460 END DO
461 END SUBROUTINE output_restart
462! **************************************************************************************************
463!> \brief Reads the density matrix from restart files and updates the starting time
464!> \param restart_section Print key section for the restart files
465!> \param rho Density matrix in AO basis
466!> \param time_index Time index to be written into the info file
467! **************************************************************************************************
468 SUBROUTINE read_restart(rtbse_env)
469 TYPE(rtbse_env_type), POINTER :: rtbse_env
470 TYPE(cp_logger_type), POINTER :: logger
471 CHARACTER(len=default_path_length) :: save_name, save_name_2
472 INTEGER :: rho_unit_nr, j
473 CHARACTER(len=17), DIMENSION(4) :: file_labels
474
475 ! This allows the delta kick and output of moment at time 0 in all cases
476 ! except the case when both imaginary and real parts of the density are read
477 rtbse_env%restart_extracted = .false.
478 logger => cp_get_default_logger()
479 ! Start by probing/loading info file
480 save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".info", my_local=.false.)
481 IF (file_exists(save_name)) THEN
482 CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
483 unit_number=rho_unit_nr)
484 READ (rho_unit_nr) rtbse_env%sim_start
485 CALL close_file(rho_unit_nr)
486 IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A31,I25,A24)') " RTBSE| Starting from timestep ", &
487 rtbse_env%sim_start, ", delta kick NOT applied"
488 ELSE
489 cpwarn("Restart required but no info file found - starting from sim_step given in input")
490 END IF
491
492 ! Default labels distinguishing up to two spin species and real/imaginary parts
493 file_labels(1) = "_SPIN_A_RE.matrix"
494 file_labels(2) = "_SPIN_A_IM.matrix"
495 file_labels(3) = "_SPIN_B_RE.matrix"
496 file_labels(4) = "_SPIN_B_IM.matrix"
497 DO j = 1, rtbse_env%n_spin
498 save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
499 extension=file_labels(2*j - 1), my_local=.false.)
500 save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
501 extension=file_labels(2*j), my_local=.false.)
502 IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN
503 CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
504 unit_number=rho_unit_nr)
505 CALL cp_fm_read_unformatted(rtbse_env%real_workspace(1), rho_unit_nr)
506 CALL close_file(rho_unit_nr)
507 CALL open_file(save_name_2, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
508 unit_number=rho_unit_nr)
509 CALL cp_fm_read_unformatted(rtbse_env%real_workspace(2), rho_unit_nr)
510 CALL close_file(rho_unit_nr)
511 CALL cp_fm_to_cfm(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), &
512 rtbse_env%rho(j))
513 rtbse_env%restart_extracted = .true.
514 ELSE
515 cpwarn("Restart without some restart matrices - starting from SCF density.")
516 END IF
517 END DO
518 END SUBROUTINE read_restart
519! **************************************************************************************************
520!> \brief Reads the moments and time traces from the save files
521!> \param rtbse_env GW-BSE environment (assumes consistent setup, i.e. a continuation run).
522!> Otherwise, the traces are set at zero
523! **************************************************************************************************
524 SUBROUTINE read_moments(rtbse_env)
525 TYPE(rtbse_env_type), POINTER :: rtbse_env
526 TYPE(cp_logger_type), POINTER :: logger
527 CHARACTER(len=default_path_length) :: save_name, save_name_2
528 INTEGER :: i, j, k, moments_unit_re, moments_unit_im, n
529 CHARACTER(len=17), DIMENSION(4) :: file_labels
530 REAL(kind=dp), DIMENSION(3) :: moments_re, moments_im
531 REAL(kind=dp) :: timestamp
532
533 logger => cp_get_default_logger()
534
535 ! Read moments from the previous run
536 ! Default labels distinguishing up to two spin species and real/imaginary parts
537 file_labels(1) = "_SPIN_A_RE.dat"
538 file_labels(2) = "_SPIN_A_IM.dat"
539 file_labels(3) = "_SPIN_B_RE.dat"
540 file_labels(4) = "_SPIN_B_IM.dat"
541 DO j = 1, rtbse_env%n_spin
542 save_name = cp_print_key_generate_filename(logger, rtbse_env%moments_section, &
543 extension=file_labels(2*j - 1), my_local=.false.)
544 save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%moments_section, &
545 extension=file_labels(2*j), my_local=.false.)
546 IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN
547 CALL open_file(save_name, file_status="OLD", file_form="FORMATTED", file_action="READ", &
548 unit_number=moments_unit_re)
549 CALL open_file(save_name_2, file_status="OLD", file_form="FORMATTED", file_action="READ", &
550 unit_number=moments_unit_im)
551 ! Extra time step for the initial one
552 DO k = rtbse_env%sim_start_orig, rtbse_env%sim_start
553 ! Determine the absolute time step - offset in memory
554 n = k - rtbse_env%sim_start_orig + 1
555 READ (moments_unit_re, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') timestamp, &
556 moments_re(1), moments_re(2), moments_re(3)
557 READ (moments_unit_im, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') timestamp, &
558 moments_im(1), moments_im(2), moments_im(3)
559 DO i = 1, 3
560 rtbse_env%moments_trace(i)%series(n) = cmplx(moments_re(i), moments_im(i), kind=dp)
561 END DO
562 rtbse_env%time_trace%series(n) = timestamp
563 END DO
564 ! Change back to atomic units in the trace
565 rtbse_env%time_trace%series(:) = rtbse_env%time_trace%series(:)/femtoseconds
566 CALL close_file(moments_unit_re)
567 CALL close_file(moments_unit_im)
568 ELSE IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
569 cpwarn("Restart without previous moments - unknown moments trace set to zero.")
570 END IF
571 END DO
572 END SUBROUTINE read_moments
573! **************************************************************************************************
574!> \brief Outputs the Fourier transform of moments stored in the environment memory to the configured file
575!> \param rtbse_env GW-BSE environment
576! **************************************************************************************************
577 SUBROUTINE output_moments_ft(rtbse_env)
578 TYPE(rtbse_env_type), POINTER :: rtbse_env
579 TYPE(cp_logger_type), POINTER :: logger
580 REAL(kind=dp), DIMENSION(:), POINTER :: omega_series, &
581 ft_real_series, &
582 ft_imag_series, &
583 value_series
584 ! Stores the data in ready for output format
585 ! - first dimension is 6 - 1 - real part along x, 2 - imag part along x, 3 - real part along y, ...
586 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: ft_full_series
587 INTEGER :: i, n, ft_unit
588#if defined (__GREENX)
589 COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_complex, &
590 moments_ft_complex
591 COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: moments_eval_complex
592#endif
593
594 logger => cp_get_default_logger()
595
596 n = rtbse_env%sim_nsteps + 2
597 NULLIFY (omega_series)
598 ALLOCATE (omega_series(n), source=0.0_dp)
599 NULLIFY (ft_real_series)
600 ALLOCATE (ft_real_series(n), source=0.0_dp)
601 NULLIFY (ft_imag_series)
602 ALLOCATE (ft_imag_series(n), source=0.0_dp)
603 NULLIFY (value_series)
604 ALLOCATE (value_series(n), source=0.0_dp)
605 ALLOCATE (ft_full_series(6, n))
606 ! Carry out for each direction independently and real and imaginary parts also independently
607 DO i = 1, 3
608 ! Real part of the value first
609 value_series(:) = real(rtbse_env%moments_trace(i)%series(:))
610 CALL ft_simple(rtbse_env%time_trace%series, value_series, ft_real_series, ft_imag_series, &
611 damping_opt=rtbse_env%ft_damping, t0_opt=rtbse_env%ft_start, subtract_initial_opt=.true.)
612 ft_full_series(2*i - 1, :) = ft_real_series(:)
613 ft_full_series(2*i, :) = ft_imag_series(:)
614 ! Now imaginary part
615 value_series(:) = aimag(rtbse_env%moments_trace(i)%series(:))
616 CALL ft_simple(rtbse_env%time_trace%series, value_series, ft_real_series, ft_imag_series, omega_series, &
617 damping_opt=rtbse_env%ft_damping, t0_opt=rtbse_env%ft_start, subtract_initial_opt=.true.)
618 ft_full_series(2*i - 1, :) = ft_full_series(2*i - 1, :) - ft_imag_series
619 ft_full_series(2*i, :) = ft_full_series(2*i, :) + ft_real_series
620 END DO
621 DEALLOCATE (value_series)
622 ! Now, write these to file
623 ft_unit = cp_print_key_unit_nr(logger, rtbse_env%ft_section, extension=".dat", &
624 file_form="FORMATTED", file_position="REWIND")
625 IF (ft_unit > 0) THEN
626 DO i = 1, n
627 WRITE (ft_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
628 omega_series(i), ft_full_series(1, i), ft_full_series(2, i), &
629 ft_full_series(3, i), ft_full_series(4, i), &
630 ft_full_series(5, i), ft_full_series(6, i)
631 END DO
632 END IF
633 CALL cp_print_key_finished_output(ft_unit, logger, rtbse_env%ft_section)
634 DEALLOCATE (ft_real_series)
635 DEALLOCATE (ft_imag_series)
636#if defined (__GREENX)
637 IF (rtbse_env%pade_requested) THEN
638 ! Report Padé refinement
639 IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A10,A27,E23.8E3,E20.8E3)') &
640 " PADE_FT| ", "Evaluation grid bounds [eV]", rtbse_env%pade_e_min, rtbse_env%pade_e_max
641 ALLOCATE (omega_complex(n))
642 ALLOCATE (moments_ft_complex(n))
643 ALLOCATE (moments_eval_complex(3, rtbse_env%pade_npoints))
644 omega_complex(:) = cmplx(omega_series(:), 0.0, kind=dp)
645 DO i = 1, 3
646 moments_ft_complex(:) = cmplx(ft_full_series(2*i - 1, :), &
647 ft_full_series(2*i, :), &
648 kind=dp)
649 ! Copy the fitting parameters
650 ! TODO : Optional direct setting of parameters?
651 CALL refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, &
652 omega_complex, moments_ft_complex, rtbse_env%pade_x_eval, moments_eval_complex(i, :))
653 END DO
654 ! Write into alternative file
655 ft_unit = cp_print_key_unit_nr(logger, rtbse_env%ft_section, extension="_PADE.dat", &
656 file_form="FORMATTED", file_position="REWIND")
657 IF (ft_unit > 0) THEN
658 DO i = 1, rtbse_env%pade_npoints
659 WRITE (ft_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
660 REAL(rtbse_env%pade_x_eval(i)), REAL(moments_eval_complex(1, i)), aimag(moments_eval_complex(1, i)), &
661 REAL(moments_eval_complex(2, i)), aimag(moments_eval_complex(2, i)), &
662 REAL(moments_eval_complex(3, i)), aimag(moments_eval_complex(3, i))
663 END DO
664 END IF
665 CALL cp_print_key_finished_output(ft_unit, logger, rtbse_env%ft_section)
666 DEALLOCATE (omega_complex)
667 DEALLOCATE (moments_ft_complex)
668 DEALLOCATE (moments_eval_complex)
669 END IF
670#endif
671 DEALLOCATE (omega_series)
672 DEALLOCATE (ft_full_series)
673 ! Perform control with gx_ac
674 END SUBROUTINE output_moments_ft
675! **************************************************************************************************
676!> \brief Refines the FT grid using Padé approximants
677!> \param x_fit Input x-variables
678!> \param y_fit Input y-variables
679!> \param x_eval Refined x-variables
680!> \param y_eval Refined y-variables
681! **************************************************************************************************
682 SUBROUTINE refine_ft(fit_e_min, fit_e_max, x_fit, y_fit, x_eval, y_eval, n_pade_opt)
683 REAL(kind=dp) :: fit_e_min, &
684 fit_e_max
685 COMPLEX(kind=dp), DIMENSION(:) :: x_fit, &
686 y_fit, &
687 x_eval, &
688 y_eval
689 INTEGER, OPTIONAL :: n_pade_opt
690#if defined (__GREENX)
691 CHARACTER(len=*), PARAMETER :: routinen = "refine_ft"
692 INTEGER :: handle, &
693 fit_start, &
694 fit_end, &
695 max_fit, &
696 n_fit, &
697 n_pade, &
698 n_eval, &
699 i
700 TYPE(params) :: pade_params
701 CALL timeset(routinen, handle)
702
703 ! Get the sizes from arrays
704 max_fit = SIZE(x_fit)
705 n_eval = SIZE(x_eval)
706
707 ! Search for the fit start and end indices
708 fit_start = -1
709 fit_end = -1
710 ! Search for the subset of FT points which is within energy limits given by
711 ! the input
712 ! Do not search when automatic request of highest energy is made
713 IF (fit_e_max < 0) fit_end = max_fit
714 DO i = 1, max_fit
715 IF (fit_start == -1 .AND. real(x_fit(i)) >= fit_e_min) fit_start = i
716 IF (fit_end == -1 .AND. real(x_fit(i)) > fit_e_max) fit_end = i - 1
717 IF (fit_start > 0 .AND. fit_end > 0) EXIT
718 END DO
719 IF (fit_start == -1) fit_start = 1
720 IF (fit_end == -1) fit_end = max_fit
721 n_fit = fit_end - fit_start + 1
722
723 n_pade = n_fit/2
724 IF (PRESENT(n_pade_opt)) n_pade = n_pade_opt
725
726 ! Warn about a large number of Padé parameters
727 IF (n_pade > 1000) THEN
728 cpwarn(é"More then 1000 Pad parameters requested - may reduce with FIT_E_MIN/FIT_E_MAX.")
729 END IF
730 ! TODO : Symmetry mode settable?
731 ! Here, we assume that ft corresponds to transform of real trace
732 pade_params = create_thiele_pade(n_pade, x_fit(fit_start:fit_end), y_fit(fit_start:fit_end), &
733 enforce_symmetry="conjugate")
734
735 ! Check whetner the splice is needed or not
736 y_eval(1:n_eval) = evaluate_thiele_pade_at(pade_params, x_eval)
737
738 CALL free_params(pade_params)
739
740 CALL timestop(handle)
741#else
742 ! Mark used
743 mark_used(fit_e_min)
744 mark_used(fit_e_max)
745 mark_used(x_fit)
746 mark_used(y_fit)
747 mark_used(x_eval)
748 mark_used(y_eval)
749 mark_used(n_pade_opt)
750 cpabort("refine_ft called but CP2K compiled without GreenX - GreenX needed.")
751#endif
752 END SUBROUTINE refine_ft
753! **************************************************************************************************
754!> \brief Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega),
755!> where i and j are provided by the configuration. The tensor element is energy dependent and
756!> has real and imaginary parts
757!> \param rtbse_env GW-BSE environment
758! **************************************************************************************************
759 SUBROUTINE output_polarizability(rtbse_env)
760 TYPE(rtbse_env_type), POINTER :: rtbse_env
761 TYPE(cp_logger_type), POINTER :: logger
762 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_series, &
763 ft_real_series, &
764 ft_imag_series, &
765 value_series
766 COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: moment_series, &
767 field_series, &
768 moment_refined, &
769 field_refined, &
770 omega_complex
771 COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: polarizability_series, &
772 polarizability_refined
773 INTEGER :: pol_unit, &
774 i, k, n, n_elems
775
776 logger => cp_get_default_logger()
777
778 n = rtbse_env%sim_nsteps + 2
779 n_elems = SIZE(rtbse_env%pol_elements, 1)
780 ! All allocations together, although could save some memory, if required by consequent deallocations
781 ALLOCATE (omega_series(n), source=0.0_dp)
782 ALLOCATE (ft_real_series(n), source=0.0_dp)
783 ALLOCATE (ft_imag_series(n), source=0.0_dp)
784 ALLOCATE (value_series(n), source=0.0_dp)
785 ALLOCATE (moment_series(n), source=cmplx(0.0, 0.0, kind=dp))
786 ALLOCATE (field_series(n), source=cmplx(0.0, 0.0, kind=dp))
787 ALLOCATE (polarizability_series(n_elems, n), source=cmplx(0.0, 0.0, kind=dp))
788
789 ! Allocations for Padé refinement
790 IF (rtbse_env%pade_requested) THEN
791 ALLOCATE (omega_complex(n), source=cmplx(0.0, 0.0, kind=dp))
792 ALLOCATE (field_refined(rtbse_env%pade_npoints), source=cmplx(0.0, 0.0, kind=dp))
793 ALLOCATE (moment_refined(rtbse_env%pade_npoints), source=cmplx(0.0, 0.0, kind=dp))
794 ALLOCATE (polarizability_refined(n_elems, rtbse_env%pade_npoints), source=cmplx(0.0, 0.0, kind=dp))
795 ELSE
796 ALLOCATE (omega_complex(0), source=cmplx(0.0, 0.0, kind=dp))
797 ALLOCATE (field_refined(0), source=cmplx(0.0, 0.0, kind=dp))
798 ALLOCATE (moment_refined(0), source=cmplx(0.0, 0.0, kind=dp))
799 ALLOCATE (polarizability_refined(0, 0), source=cmplx(0.0, 0.0, kind=dp))
800 END IF
801
802 DO k = 1, n_elems
803 ! The moment ft
804 ! Real part
805 value_series(:) = real(rtbse_env%moments_trace(rtbse_env%pol_elements(k, 1))%series(:))
806 CALL ft_simple(rtbse_env%time_trace%series, value_series, ft_real_series, ft_imag_series, &
807 damping_opt=rtbse_env%ft_damping, t0_opt=rtbse_env%ft_start, subtract_initial_opt=.true.)
808 moment_series(:) = cmplx(ft_real_series(:), ft_imag_series(:), kind=dp)
809 ! Imaginary part
810 value_series(:) = aimag(rtbse_env%moments_trace(rtbse_env%pol_elements(k, 1))%series(:))
811 CALL ft_simple(rtbse_env%time_trace%series, value_series, ft_real_series, ft_imag_series, omega_series, &
812 damping_opt=rtbse_env%ft_damping, t0_opt=rtbse_env%ft_start, subtract_initial_opt=.true.)
813 moment_series(:) = moment_series(:) + cmplx(-ft_imag_series(:), ft_real_series(:), kind=dp)
814
815 ! Calculate the field transform - store it in ft_real_series
816 IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse) THEN
817 ! Only divide by constant magnitude in atomic units
818 ! TODO : Fix for field with more than one direction
819 field_series(:) = cmplx( &
820 rtbse_env%dft_control%rtp_control%delta_pulse_direction(rtbse_env%pol_elements(k, 2))* &
821 rtbse_env%dft_control%rtp_control%delta_pulse_scale, 0.0, kind=dp)
822 ELSE
823 ! Calculate the transform of the field as well and divide by it
824 ! The field FT is not damped - assume field is localised in time
825 ! The field is strictly real
826 CALL ft_simple(rtbse_env%time_trace%series, rtbse_env%field_trace(rtbse_env%pol_elements(k, 2))%series, &
827 ft_real_series, ft_imag_series, omega_series, &
828 0.0_dp, rtbse_env%ft_start, .false.)
829 ! Regularization for the case when ft_series becomes identically zero - TODO : Set in config
830 field_series(:) = cmplx(ft_real_series(:), ft_imag_series(:) + 1.0e-20, kind=dp)
831 END IF
832 ! Divide to get the polarizability series
833 ! Real part
834 polarizability_series(k, :) = moment_series(:)/field_series(:)
835
836 ! Refinement of the grid via Padé approximants, if requested
837 IF (rtbse_env%pade_requested) THEN
838 omega_complex(:) = cmplx(omega_series(:), 0.0, kind=dp)
839 ! No need to refine the field if it is known exactly
840 ! Also GreenX has trouble fitting constant functions
841 IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse) THEN
842 field_refined(:) = cmplx( &
843 rtbse_env%dft_control%rtp_control%delta_pulse_direction(rtbse_env%pol_elements(k, 2))* &
844 rtbse_env%dft_control%rtp_control%delta_pulse_scale, 0.0, kind=dp)
845 ELSE
846 CALL refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, omega_complex, field_series, &
847 rtbse_env%pade_x_eval, field_refined)
848 END IF
849 CALL refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, omega_complex, moment_series, &
850 rtbse_env%pade_x_eval, moment_refined)
851 polarizability_refined(k, :) = moment_refined(:)/field_refined(:)
852 END IF
853 END DO
854
855 ! Change units to eV for energy
856 ! use value_series for energy and moment_series for polarizability
857 value_series(:) = omega_series(:)*evolt
858 ! Print out the polarizability to a file
859 pol_unit = cp_print_key_unit_nr(logger, rtbse_env%pol_section, extension=".dat", &
860 file_form="FORMATTED", file_position="REWIND")
861 ! Printing for both the stdout and separate file
862 IF (pol_unit > 0) THEN
863 IF (pol_unit == rtbse_env%unit_nr) THEN
864 ! Print the stdout preline
865 WRITE (pol_unit, '(A16)', advance="no") " POLARIZABILITY|"
866 ELSE
867 ! Print also the energy in atomic units
868 WRITE (pol_unit, '(A1,A19)', advance="no") "#", "omega [a.u.]"
869 END IF
870 ! Common - print the energy in eV
871 WRITE (pol_unit, '(A20)', advance="no") "Energy [eV]"
872 ! Print a header for each polarizability element
873 DO k = 1, n_elems - 1
874 WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)', advance="no") &
875 "Real pol.", rtbse_env%pol_elements(k, 1), rtbse_env%pol_elements(k, 2), &
876 "Imag pol.", rtbse_env%pol_elements(k, 1), rtbse_env%pol_elements(k, 2)
877 END DO
878 WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)') &
879 "Real pol.", rtbse_env%pol_elements(n_elems, 1), rtbse_env%pol_elements(n_elems, 2), &
880 "Imag pol.", rtbse_env%pol_elements(n_elems, 1), rtbse_env%pol_elements(n_elems, 2)
881 DO i = 1, n
882 IF (pol_unit == rtbse_env%unit_nr) THEN
883 ! Print the stdout preline
884 WRITE (pol_unit, '(A16)', advance="no") " POLARIZABILITY|"
885 ELSE
886 ! omega in a.u.
887 WRITE (pol_unit, '(E20.8E3)', advance="no") omega_series(i)
888 END IF
889 ! Common values
890 WRITE (pol_unit, '(E20.8E3)', advance="no") value_series(i)
891 DO k = 1, n_elems - 1
892 WRITE (pol_unit, '(E20.8E3,E20.8E3)', advance="no") &
893 REAL(polarizability_series(k, i)), aimag(polarizability_series(k, i))
894 END DO
895 ! Print the final value and advance
896 WRITE (pol_unit, '(E20.8E3,E20.8E3)') &
897 REAL(polarizability_series(n_elems, i)), aimag(polarizability_series(n_elems, i))
898 END DO
899 CALL cp_print_key_finished_output(pol_unit, logger, rtbse_env%pol_section)
900 END IF
901#if defined(__GREENX)
902 ! Print out the refined polarizability to a file
903 pol_unit = cp_print_key_unit_nr(logger, rtbse_env%pol_section, extension="_PADE.dat", &
904 file_form="FORMATTED", file_position="REWIND")
905 ! Printing for both the stdout and separate file
906 IF (pol_unit > 0 .AND. rtbse_env%pade_requested) THEN
907 IF (pol_unit == rtbse_env%unit_nr) THEN
908 ! Print the stdout preline
909 WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
910 ELSE
911 ! Print also the energy in atomic units
912 WRITE (pol_unit, '(A1,A19)', advance="no") "#", "omega [a.u.]"
913 END IF
914 ! Common - print the energy in eV
915 WRITE (pol_unit, '(A20)', advance="no") "Energy [eV]"
916 ! Print a header for each polarizability element
917 DO k = 1, n_elems - 1
918 WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)', advance="no") &
919 "Real pol.", rtbse_env%pol_elements(k, 1), rtbse_env%pol_elements(k, 2), &
920 "Imag pol.", rtbse_env%pol_elements(k, 1), rtbse_env%pol_elements(k, 2)
921 END DO
922 WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)') &
923 "Real pol.", rtbse_env%pol_elements(n_elems, 1), rtbse_env%pol_elements(n_elems, 2), &
924 "Imag pol.", rtbse_env%pol_elements(n_elems, 1), rtbse_env%pol_elements(n_elems, 2)
925 DO i = 1, SIZE(rtbse_env%pade_x_eval)
926 IF (pol_unit == rtbse_env%unit_nr) THEN
927 ! Print the stdout preline
928 WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
929 ELSE
930 ! omega in a.u.
931 WRITE (pol_unit, '(E20.8E3)', advance="no") real(rtbse_env%pade_x_eval(i), kind=dp)
932 END IF
933 ! Common values
934 WRITE (pol_unit, '(E20.8E3)', advance="no") real(rtbse_env%pade_x_eval(i), kind=dp)*evolt
935 DO k = 1, n_elems - 1
936 WRITE (pol_unit, '(E20.8E3,E20.8E3)', advance="no") &
937 REAL(polarizability_refined(k, i)), aimag(polarizability_refined(k, i))
938 END DO
939 ! Print the final value and advance
940 WRITE (pol_unit, '(E20.8E3,E20.8E3)') &
941 REAL(polarizability_refined(n_elems, i)), aimag(polarizability_refined(n_elems, i))
942 END DO
943 CALL cp_print_key_finished_output(pol_unit, logger, rtbse_env%pol_section)
944 END IF
945#endif
946
947 DEALLOCATE (value_series)
948 DEALLOCATE (ft_real_series)
949 DEALLOCATE (ft_imag_series)
950
951 DEALLOCATE (field_series)
952 DEALLOCATE (moment_series)
953
954 DEALLOCATE (omega_series)
955 DEALLOCATE (polarizability_series)
956
957 DEALLOCATE (omega_complex)
958 DEALLOCATE (field_refined)
959 DEALLOCATE (moment_refined)
960 DEALLOCATE (polarizability_refined)
961 END SUBROUTINE output_polarizability
962! **************************************************************************************************
963!> \brief Naively calculates the Fourier transform - it is not the bottleneck of this calculation
964!> \param time_series Timestamps in atomic units of time
965!> \param value_series Values to be Fourier transformed - moments, field etc. So far only real.
966!> \param omega_series Array to be filled by sampled values of frequency
967!> \param result_series FT of the value series - real values (cosines)
968!> \param iresult_series FT of the value series - imaginary values (sines)
969!> \param damping_opt Supply custom exponential damping - default is 4.0/totalTime, i.e. ratio
970!> of last and first element in windowed value series is reduced by e^(-4)
971!> \param t0_opt Carry the FT only starting from certain time - allows for exclusion of trace before
972!> the pulse application etc.
973!> \author Stepan Marek
974!> \date 09.2024
975! **************************************************************************************************
976 ! So far only for defined one dimensional series
977 SUBROUTINE ft_simple(time_series, value_series, result_series, iresult_series, omega_series, &
978 damping_opt, t0_opt, subtract_initial_opt)
979 REAL(kind=dp), DIMENSION(:) :: time_series
980 REAL(kind=dp), DIMENSION(:) :: value_series
981 REAL(kind=dp), DIMENSION(:) :: result_series
982 REAL(kind=dp), DIMENSION(:) :: iresult_series
983 REAL(kind=dp), DIMENSION(:), OPTIONAL :: omega_series
984 REAL(kind=dp), OPTIONAL :: damping_opt
985 REAL(kind=dp), OPTIONAL :: t0_opt
986 LOGICAL, OPTIONAL :: subtract_initial_opt
987 CHARACTER(len=*), PARAMETER :: routinen = "ft_simple"
988 INTEGER :: n, i, j, t0_i, j_wrap, handle
989 REAL(kind=dp) :: t0, delta_t, delta_omega, damping, &
990 value_subtract
991 LOGICAL :: subtract_initial
992
993 CALL timeset(routinen, handle)
994
995 n = SIZE(time_series)
996
997 t0_i = 1
998 IF (PRESENT(t0_opt)) THEN
999 ! Determine the index at which we start applying the damping
1000 ! Also the index around which we fold around
1001 DO i = 1, n
1002 ! Increase until we break or reach the end of the time series
1003 t0_i = i
1004 IF (time_series(i) >= t0_opt) THEN
1005 EXIT
1006 END IF
1007 END DO
1008 END IF
1009
1010 t0 = time_series(t0_i)
1011
1012 ! Default damping so that at the end of the time series, divide value by e^-4
1013 damping = 4.0_dp/(time_series(n) - time_series(t0_i))
1014 IF (PRESENT(damping_opt)) THEN
1015 ! Damping is given a time in which the moments reduce by factor of 1/e
1016 IF (damping_opt > 0.0_dp) damping = 1.0_dp/damping_opt
1017 ! Special case - zero damping
1018 IF (damping_opt == 0.0_dp) damping = 0.0_dp
1019 END IF
1020
1021 IF (PRESENT(subtract_initial_opt)) THEN
1022 subtract_initial = subtract_initial_opt
1023 ELSE
1024 subtract_initial = .true.
1025 END IF
1026
1027 ! Construct the grid
1028 ! Assume series equidistant
1029 delta_t = time_series(2) - time_series(1)
1030 delta_omega = twopi/(time_series(n) - time_series(1))
1031 ! Subtract initial value, if requested (default is to subtract the value)
1032 value_subtract = 0.0_dp
1033 IF (subtract_initial) value_subtract = value_series(1)
1034 DO i = 1, n
1035 result_series(i) = 0.0_dp
1036 iresult_series(i) = 0.0_dp
1037 IF (PRESENT(omega_series)) omega_series(i) = delta_omega*(i - 1)
1038 DO j = 1, n
1039 j_wrap = modulo(j + t0_i - 2, n) + 1
1040 result_series(i) = result_series(i) + cos(twopi*(i - 1)*(j - 1)/n)* &
1041 exp(-damping*delta_t*(j - 1))*(value_series(j_wrap) - value_subtract)
1042 iresult_series(i) = iresult_series(i) + sin(twopi*(i - 1)*(j - 1)/n)* &
1043 exp(-damping*delta_t*(j - 1))*(value_series(j_wrap) - value_subtract)
1044 END DO
1045 END DO
1046 result_series(:) = delta_t*result_series(:)
1047 iresult_series(:) = delta_t*iresult_series(:)
1048
1049 CALL timestop(handle)
1050
1051 END SUBROUTINE
1052END MODULE rt_bse_io
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:75
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:813
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:761
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:501
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition cp_fm_basic_linalg.F:1695
cp_fm_basic_linalg::cp_fm_norm
real(kind=dp) function, public cp_fm_norm(matrix, mode)
norm of matrix using (p)dlange
Definition cp_fm_basic_linalg.F:2595
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2146
cp_fm_types::cp_fm_read_unformatted
subroutine, public cp_fm_read_unformatted(fm, unit)
...
Definition cp_fm_types.F:2377
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
cp_fm_types::cp_fm_write_formatted
subroutine, public cp_fm_write_formatted(fm, unit, header, value_format)
Write out a full matrix in plain text.
Definition cp_fm_types.F:2250
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:73
cp_output_handling::medium_print_level
integer, parameter, public medium_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:754
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_bch
integer, parameter, public do_bch
Definition input_constants.F:925
input_constants::use_rt_restart
integer, parameter, public use_rt_restart
Definition input_constants.F:935
input_constants::do_exact
integer, parameter, public do_exact
Definition input_constants.F:925
input_constants::rtp_bse_ham_g0w0
integer, parameter, public rtp_bse_ham_g0w0
Definition input_constants.F:947
input_constants::rtp_bse_ham_ks
integer, parameter, public rtp_bse_ham_ks
Definition input_constants.F:947
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::femtoseconds
real(kind=dp), parameter, public femtoseconds
Definition physcon.F:153
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
rt_bse_io
Input/output from the propagation via RT-BSE method.
Definition rt_bse_io.F:13
rt_bse_io::print_rtbse_header_info
subroutine, public print_rtbse_header_info(rtbse_env)
Writes the header and basic info to the standard output.
Definition rt_bse_io.F:93
rt_bse_io::output_mos_covariant
subroutine, public output_mos_covariant(rtbse_env, ham, print_key_section)
Outputs the matrix in MO basis for matrix components corresponding to covariant representation,...
Definition rt_bse_io.F:245
rt_bse_io::output_restart
subroutine, public output_restart(rtbse_env, rho, time_index)
Outputs the restart info (last finished iteration step) + restard density matrix.
Definition rt_bse_io.F:425
rt_bse_io::print_etrs_info
subroutine, public print_etrs_info(rtbse_env, step, metric)
Writes the update after single etrs iteration - only for log level > medium.
Definition rt_bse_io.F:137
rt_bse_io::read_field
subroutine, public read_field(rtbse_env)
Reads the field from the files provided by input - useful for the continuation run.
Definition rt_bse_io.F:313
rt_bse_io::output_field
subroutine, public output_field(rtbse_env)
Prints the current field components into a file provided by input.
Definition rt_bse_io.F:283
rt_bse_io::print_etrs_info_header
subroutine, public print_etrs_info_header(rtbse_env)
Writes the header for the etrs iteration updates - only for log level > medium.
Definition rt_bse_io.F:154
rt_bse_io::output_mos_contravariant
subroutine, public output_mos_contravariant(rtbse_env, rho, print_key_section)
Outputs the matrix in MO basis for matrix coefficients corresponding to contravariant operator,...
Definition rt_bse_io.F:195
rt_bse_io::read_moments
subroutine, public read_moments(rtbse_env)
Reads the moments and time traces from the save files.
Definition rt_bse_io.F:525
rt_bse_io::output_moments_ft
subroutine, public output_moments_ft(rtbse_env)
Outputs the Fourier transform of moments stored in the environment memory to the configured file.
Definition rt_bse_io.F:578
rt_bse_io::read_restart
subroutine, public read_restart(rtbse_env)
Reads the density matrix from restart files and updates the starting time.
Definition rt_bse_io.F:469
rt_bse_io::print_timestep_info
subroutine, public print_timestep_info(rtbse_env, step, convergence, electron_num_re, etrs_num)
Writes the header for the etrs iteration updates - only for log level > medium.
Definition rt_bse_io.F:169
rt_bse_io::output_polarizability
subroutine, public output_polarizability(rtbse_env)
Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega),...
Definition rt_bse_io.F:760
rt_bse_io::output_moments
subroutine, public output_moments(rtbse_env, rho)
Outputs the expectation value of moments from a given density matrix.
Definition rt_bse_io.F:347
rt_bse_types
Data storage and other types for propagation via RT-BSE method.
Definition rt_bse_types.F:13
rt_bse_types::multiply_cfm_fm
subroutine, public multiply_cfm_fm(trans_c, trans_r, na, nb, nc, alpha, matrix_c, matrix_r, beta, res)
Multiplies complex matrix by a real matrix from the right.
Definition rt_bse_types.F:941
rt_bse_types::multiply_fm_cfm
subroutine, public multiply_fm_cfm(trans_r, trans_c, na, nb, nc, alpha, matrix_r, matrix_c, beta, res)
Multiplies real matrix by a complex matrix from the right.
Definition rt_bse_types.F:892
cp_cfm_types::cp_cfm_p_type
Just to build arrays of pointers to matrices.
Definition cp_cfm_types.F:78
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_fm_types::cp_fm_p_type
just to build arrays of pointers to matrices
Definition cp_fm_types.F:129
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
rt_bse_types::rtbse_env_type
Definition rt_bse_types.F:181