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