Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_read_unformatted, &
16 : cp_fm_write_unformatted, &
17 : cp_fm_write_formatted
18 : USE cp_cfm_types, ONLY: cp_cfm_type, &
19 : cp_fm_to_cfm, &
20 : cp_cfm_to_fm
21 : USE kinds, ONLY: dp, &
22 : default_path_length
23 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace, &
24 : cp_fm_transpose
25 : USE cp_log_handling, ONLY: cp_logger_type, &
26 : cp_get_default_logger
27 : USE cp_output_handling, ONLY: cp_print_key_unit_nr, &
28 : cp_print_key_finished_output, &
29 : cp_print_key_generate_filename, &
30 : low_print_level, &
31 : medium_print_level
32 : USE input_section_types, ONLY: section_vals_type
33 : USE rt_bse_types, ONLY: rtbse_env_type, &
34 : multiply_cfm_fm, &
35 : multiply_fm_cfm
36 : USE cp_files, ONLY: open_file, &
37 : file_exists, &
38 : close_file
39 : USE input_constants, ONLY: do_exact, &
40 : do_bch, &
41 : evGW0, &
42 : rtp_bse_ham_gw, &
43 : rtp_bse_ham_ks, &
44 : use_rt_restart, &
45 : restart_guess
46 : USE qs_environment_types, ONLY: get_qs_env
47 : USE scf_control_types, ONLY: scf_control_type
48 : USE physcon, ONLY: evolt, femtoseconds
49 : USE rt_propagation_output, ONLY: print_moments, &
50 : print_rt_file, &
51 : rt_file_comp_real
52 :
53 : #include "../base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse_io"
60 :
61 : ! RESTART.trace on-disk format version - bump whenever the header or record layout changes
62 : INTEGER, PARAMETER, PRIVATE :: restart_trace_version = 1
63 :
64 : #:include "rt_bse_macros.fypp"
65 :
66 : PUBLIC :: output_moments, &
67 : output_field, &
68 : read_field, &
69 : output_mos_contravariant, &
70 : output_mos_covariant, &
71 : output_restart, &
72 : read_restart, &
73 : output_restart_linearized, &
74 : read_restart_info, &
75 : read_restart_trace, &
76 : read_restart_density, &
77 : read_restart_C, &
78 : check_restart_eps_consistency, &
79 : print_etrs_info_header, &
80 : print_etrs_info, &
81 : print_timestep_info, &
82 : print_rtbse_header_info
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief Writes the header and basic info to the standard output
88 : !> \param rtbse_env Entry point - rtbse environment
89 : ! **************************************************************************************************
90 14 : SUBROUTINE print_rtbse_header_info(rtbse_env)
91 : TYPE(rtbse_env_type) :: rtbse_env
92 : TYPE(cp_logger_type), POINTER :: logger
93 :
94 14 : logger => cp_get_default_logger()
95 :
96 14 : IF (rtbse_env%unit_nr > 0) THEN
97 7 : WRITE (rtbse_env%unit_nr, *) ''
98 : WRITE (rtbse_env%unit_nr, '(A)') ' /-----------------------------------------------'// &
99 7 : '------------------------------\'
100 : WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
101 7 : ' |'
102 : WRITE (rtbse_env%unit_nr, '(A)') ' | Real Time Bethe-Salpeter Propagation'// &
103 7 : ' |'
104 : WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
105 7 : ' |'
106 : WRITE (rtbse_env%unit_nr, '(A)') ' \-----------------------------------------------'// &
107 7 : '------------------------------/'
108 7 : WRITE (rtbse_env%unit_nr, *) ''
109 :
110 : ! Methods used
111 7 : WRITE (rtbse_env%unit_nr, '(A19)', advance="no") ' Exponential method'
112 14 : SELECT CASE (rtbse_env%mat_exp_method)
113 : CASE (do_bch)
114 7 : WRITE (rtbse_env%unit_nr, '(A61)') 'BCH'
115 : CASE (do_exact)
116 7 : WRITE (rtbse_env%unit_nr, '(A61)') 'EXACT'
117 : END SELECT
118 :
119 7 : WRITE (rtbse_env%unit_nr, '(A22)', advance="no") ' Reference Hamiltonian'
120 14 : SELECT CASE (rtbse_env%ham_reference_type)
121 : CASE (rtp_bse_ham_gw)
122 7 : IF (rtbse_env%bs_env%gw_flavour == evGW0) THEN
123 0 : WRITE (rtbse_env%unit_nr, '(A58)') 'evGW0'
124 : ELSE
125 7 : WRITE (rtbse_env%unit_nr, '(A58)') 'G0W0'
126 : END IF
127 : CASE (rtp_bse_ham_ks)
128 7 : WRITE (rtbse_env%unit_nr, '(A58)') 'Kohn-Sham'
129 : END SELECT
130 :
131 7 : WRITE (rtbse_env%unit_nr, '(A18,L62)') ' Apply delta pulse', &
132 14 : rtbse_env%dft_control%rtp_control%apply_delta_pulse
133 :
134 7 : WRITE (rtbse_env%unit_nr, '(A)') ''
135 : END IF
136 :
137 14 : END SUBROUTINE print_rtbse_header_info
138 :
139 : ! **************************************************************************************************
140 : !> \brief Writes the update after single etrs iteration - only for log level > medium
141 : !> \param rtbse_env Entry point - rtbse environment
142 : ! **************************************************************************************************
143 1732 : SUBROUTINE print_etrs_info(rtbse_env, step, metric)
144 : TYPE(rtbse_env_type) :: rtbse_env
145 : INTEGER :: step
146 : REAL(kind=dp) :: metric
147 : TYPE(cp_logger_type), POINTER :: logger
148 :
149 1732 : logger => cp_get_default_logger()
150 :
151 1732 : IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN
152 0 : WRITE (rtbse_env%unit_nr, '(A7,I5, E20.8E3)') ' RTBSE|', step, metric
153 : END IF
154 :
155 1732 : END SUBROUTINE print_etrs_info
156 : ! **************************************************************************************************
157 : !> \brief Writes the header for the etrs iteration updates - only for log level > medium
158 : !> \param rtbse_env Entry point - rtbse environment
159 : ! **************************************************************************************************
160 678 : SUBROUTINE print_etrs_info_header(rtbse_env)
161 : TYPE(rtbse_env_type) :: rtbse_env
162 : TYPE(cp_logger_type), POINTER :: logger
163 :
164 678 : logger => cp_get_default_logger()
165 :
166 678 : IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN
167 0 : WRITE (rtbse_env%unit_nr, '(A13, A20)') ' RTBSE| Iter.', 'Convergence'
168 : END IF
169 :
170 678 : END SUBROUTINE print_etrs_info_header
171 : ! **************************************************************************************************
172 : !> \brief Writes the summary line of a completed propagation timestep
173 : !> \param rtbse_env Entry point - rtbse environment
174 : !> \param step Index of the completed timestep
175 : !> \param electron_num_re Real part of the electron number, one entry per spin channel
176 : !> \param convergence Optional convergence metric reached by the ETRS iteration
177 : !> \param etrs_num Optional number of ETRS iterations used in this timestep
178 : !> \param step_walltime Optional wall time of this timestep, in seconds
179 : ! **************************************************************************************************
180 1508 : SUBROUTINE print_timestep_info(rtbse_env, step, electron_num_re, convergence, etrs_num, &
181 : step_walltime)
182 : TYPE(rtbse_env_type) :: rtbse_env
183 : INTEGER :: step
184 : REAL(kind=dp), DIMENSION(:), INTENT(IN) :: electron_num_re
185 : REAL(kind=dp), OPTIONAL :: convergence
186 : INTEGER, OPTIONAL :: etrs_num
187 : REAL(kind=dp), OPTIONAL :: step_walltime
188 : TYPE(cp_logger_type), POINTER :: logger
189 : LOGICAL :: flag_lrrtbse
190 : INTEGER :: nch
191 :
192 1508 : logger => cp_get_default_logger()
193 : ! one electron-number column per spin channel (1 = closed shell, 2 = open shell alpha/beta)
194 1508 : nch = SIZE(electron_num_re)
195 :
196 1508 : IF (.NOT. PRESENT(convergence) .OR. .NOT. PRESENT(etrs_num)) THEN
197 : flag_lrrtbse = .TRUE.
198 : ELSE
199 678 : flag_lrrtbse = .FALSE.
200 : END IF
201 :
202 1508 : IF (logger%iter_info%print_level > low_print_level .AND. rtbse_env%unit_nr > 0) THEN
203 754 : IF (flag_lrrtbse) THEN
204 415 : IF (step == 0) THEN
205 22 : IF (PRESENT(step_walltime)) THEN
206 : WRITE (rtbse_env%unit_nr, '(A45,T70,F11.3)') &
207 22 : " RTBSE| Estimated runtime for propagation [s]", &
208 44 : step_walltime*REAL(rtbse_env%sim_nsteps, dp)
209 : WRITE (rtbse_env%unit_nr, '(A)') &
210 22 : " RTBSE|"
211 22 : IF (nch == 1) THEN
212 : WRITE (rtbse_env%unit_nr, '(A23,T27,A13,T66,A15)') &
213 19 : " RTBSE| Simulation step", "Step time [s]", "Electron number"
214 : ELSE
215 : ! T67/T88 (not T66/T86): each α/β is 2 bytes but 1 display col, so the byte
216 : ! anchor is +1 per unicode char to right-align ')' under the value's last digit.
217 : ! T anchors carry +1 byte per α/β before the ')' (α col → +1, β col → +2),
218 : ! since each is 2 bytes but 1 display col; right-aligns ')' on the last digit.
219 : WRITE (rtbse_env%unit_nr, '(A23,T27,A13,T47,A15,T68,A15)') &
220 3 : " RTBSE| Simulation step", "Step time [s]", "El. number (α)", "El. number (β)"
221 : END IF
222 : ELSE
223 0 : IF (nch == 1) THEN
224 0 : WRITE (rtbse_env%unit_nr, '(A23,T66,A15)') " RTBSE| Simulation step", "Electron number"
225 : ELSE
226 : WRITE (rtbse_env%unit_nr, '(A23,T47,A15,T68,A15)') &
227 0 : " RTBSE| Simulation step", "El. number (α)", "El. number (β)"
228 : END IF
229 : END IF
230 : END IF
231 415 : IF (PRESENT(step_walltime)) THEN
232 415 : IF (nch == 1) THEN
233 : WRITE (rtbse_env%unit_nr, '(A7,I16,T30,F10.3,T69,E12.3E3)') &
234 355 : ' RTBSE|', step, step_walltime, electron_num_re(1)
235 : ELSE
236 : WRITE (rtbse_env%unit_nr, '(A7,I16,T30,F10.3,T49,E12.3E3,T69,E12.3E3)') &
237 60 : ' RTBSE|', step, step_walltime, electron_num_re(1), electron_num_re(2)
238 : END IF
239 : ELSE
240 0 : IF (nch == 1) THEN
241 0 : WRITE (rtbse_env%unit_nr, '(A7,I16,T61,E20.8E3)') ' RTBSE|', step, electron_num_re(1)
242 : ELSE
243 : WRITE (rtbse_env%unit_nr, '(A7,I16,T49,E12.3E3,T69,E12.3E3)') &
244 0 : ' RTBSE|', step, electron_num_re(1), electron_num_re(2)
245 : END IF
246 : END IF
247 : ELSE
248 339 : WRITE (rtbse_env%unit_nr, '(A23,A20,A20,A17)') " RTBSE| Simulation step", "Convergence", &
249 678 : "Electron number", "ETRS Iterations"
250 339 : WRITE (rtbse_env%unit_nr, '(A7,I16,E20.8E3,E20.8E3,I17)') ' RTBSE|', step, convergence, &
251 678 : electron_num_re(1), etrs_num
252 : END IF
253 : END IF
254 :
255 1508 : END SUBROUTINE print_timestep_info
256 :
257 : ! **************************************************************************************************
258 : !> \brief Outputs the matrix in MO basis for matrix coefficients corresponding to contravariant
259 : !> operator, i.e. density matrix
260 : !> \param rtbse_env Entry point - gwbse environment
261 : !> \param rho Density matrix in AO basis
262 : !> \param rtp_section RTP input section
263 : ! **************************************************************************************************
264 1508 : SUBROUTINE output_mos_contravariant(rtbse_env, rho, print_key_section)
265 : TYPE(rtbse_env_type) :: rtbse_env
266 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
267 : TYPE(section_vals_type), POINTER :: print_key_section
268 : TYPE(cp_logger_type), POINTER :: logger
269 : INTEGER :: j, rho_unit_re, rho_unit_im
270 : CHARACTER(len=14), DIMENSION(4) :: file_labels
271 :
272 1508 : file_labels(1) = "_SPIN_A_RE.dat"
273 1508 : file_labels(2) = "_SPIN_A_IM.dat"
274 1508 : file_labels(3) = "_SPIN_B_RE.dat"
275 1508 : file_labels(4) = "_SPIN_B_IM.dat"
276 1508 : logger => cp_get_default_logger()
277 :
278 : ! In the linearized RT-BSE active-MO path, rho is already in the MO basis
279 : ! restricted to the active window (sized mo_active x mo_active). Dump it
280 : ! directly without the AO-side C^T S * rho * S C transformation.
281 1508 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) THEN
282 1780 : DO j = 1, rtbse_env%n_spin
283 950 : rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1))
284 950 : rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j))
285 950 : CALL cp_cfm_to_fm(rho(j), rtbse_env%real_workspace_mo(1), rtbse_env%real_workspace_mo(2))
286 950 : CALL cp_fm_write_formatted(rtbse_env%real_workspace_mo(1), rho_unit_re)
287 950 : CALL cp_fm_write_formatted(rtbse_env%real_workspace_mo(2), rho_unit_im)
288 950 : CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section)
289 1780 : CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section)
290 : END DO
291 830 : RETURN
292 : END IF
293 : ! Start by multiplying the current density by MOS
294 1356 : DO j = 1, rtbse_env%n_spin
295 678 : rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1))
296 678 : rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j))
297 : ! Transform the density matrix into molecular orbitals basis and print it out
298 : ! S * rho
299 : CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
300 : 1.0_dp, rtbse_env%S_fm, rho(j), &
301 678 : 0.0_dp, rtbse_env%rho_workspace(1))
302 : ! C^T * S * rho
303 : CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
304 : 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), rtbse_env%rho_workspace(1), &
305 678 : 0.0_dp, rtbse_env%rho_workspace(2))
306 : ! C^T * S * rho * S
307 : CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
308 : 1.0_dp, rtbse_env%rho_workspace(2), rtbse_env%S_fm, &
309 678 : 0.0_dp, rtbse_env%rho_workspace(1))
310 : ! C^T * S * rho * S * C
311 : CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
312 : 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), &
313 678 : 0.0_dp, rtbse_env%rho_workspace(2))
314 : ! Print real and imaginary parts separately
315 : CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), &
316 678 : rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
317 678 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re)
318 678 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im)
319 678 : CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section)
320 1356 : CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section)
321 : END DO
322 : END SUBROUTINE output_mos_contravariant
323 : ! **************************************************************************************************
324 : !> \brief Outputs the matrix in MO basis for matrix components corresponding to covariant representation,
325 : !> i.e. the Hamiltonian matrix
326 : !> \param rtbse_env Entry point - gwbse environment
327 : !> \param cohsex cohsex matrix in AO basis, covariant representation
328 : !> \param rtp_section RTP input section
329 : ! **************************************************************************************************
330 0 : SUBROUTINE output_mos_covariant(rtbse_env, ham, print_key_section)
331 : TYPE(rtbse_env_type) :: rtbse_env
332 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham
333 : TYPE(section_vals_type), POINTER :: print_key_section
334 : TYPE(cp_logger_type), POINTER :: logger
335 : INTEGER :: j, rho_unit_re, rho_unit_im
336 : CHARACTER(len=21), DIMENSION(4) :: file_labels
337 :
338 0 : file_labels(1) = "_SPIN_A_RE.dat"
339 0 : file_labels(2) = "_SPIN_A_IM.dat"
340 0 : file_labels(3) = "_SPIN_B_RE.dat"
341 0 : file_labels(4) = "_SPIN_B_IM.dat"
342 0 : logger => cp_get_default_logger()
343 0 : DO j = 1, rtbse_env%n_spin
344 0 : rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1))
345 0 : rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j))
346 : ! C^T * cohsex
347 : CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
348 : 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), ham(j), &
349 0 : 0.0_dp, rtbse_env%rho_workspace(1))
350 : ! C^T * cohsex * C
351 : CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
352 : 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), &
353 0 : 0.0_dp, rtbse_env%rho_workspace(2))
354 : ! Print real and imaginary parts separately
355 : CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), &
356 0 : rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
357 0 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re)
358 0 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im)
359 0 : CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section)
360 0 : CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section)
361 : END DO
362 0 : END SUBROUTINE output_mos_covariant
363 : ! **************************************************************************************************
364 : !> \brief Prints the current field components into a file provided by input
365 : !> \param rtbse_env Entry point - gwbse environment
366 : !> \param rtp_section RTP input section
367 : ! **************************************************************************************************
368 1562 : SUBROUTINE output_field(rtbse_env, append_opt)
369 : TYPE(rtbse_env_type) :: rtbse_env
370 : LOGICAL, OPTIONAL :: append_opt
371 : TYPE(cp_logger_type), POINTER :: logger
372 : INTEGER :: field_unit, n, i
373 : LOGICAL :: append
374 :
375 : ! Figure out whether we are appending or not, true by default
376 1562 : append = .TRUE.
377 1562 : IF (PRESENT(append_opt)) append = .FALSE.
378 :
379 : ! First, write the current field to memory
380 : ! Need the absolute index
381 1562 : n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1
382 6248 : DO i = 1, 3
383 6248 : rtbse_env%field_trace(i, n) = CMPLX(rtbse_env%field(i), 0.0, kind=dp)
384 : END DO
385 1562 : rtbse_env%time_trace(n) = rtbse_env%sim_time
386 :
387 : ! Now, continue to file output
388 : ! Get logger
389 1562 : logger => cp_get_default_logger()
390 : ! Get file descriptor
391 1562 : field_unit = cp_print_key_unit_nr(logger, rtbse_env%field_section, extension=".dat")
392 1562 : IF (append) THEN
393 : CALL print_rt_file(field_unit, xvals=rtbse_env%time_trace(n:n), &
394 : yvals=rtbse_env%field_trace(:, n:n), &
395 1552 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
396 : ELSE
397 : CALL print_rt_file(field_unit, [ &
398 : "# Time [fs]", &
399 : " field x [at.u.]", &
400 : " field y [at.u.]", &
401 : " field z [at.u.]"], &
402 : rtbse_env%time_trace(n:n), rtbse_env%field_trace(:, n:n), &
403 50 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
404 : END IF
405 1562 : CALL cp_print_key_finished_output(field_unit, logger, rtbse_env%field_section)
406 :
407 1562 : END SUBROUTINE output_field
408 : ! **************************************************************************************************
409 : !> \brief Reads the field from the files provided by input - useful for the continuation run
410 : !> \param rtbse_env Entry point - gwbse environment
411 : !> \param rtp_section RTP input section
412 : ! **************************************************************************************************
413 14 : SUBROUTINE read_field(rtbse_env)
414 : TYPE(rtbse_env_type) :: rtbse_env
415 : TYPE(cp_logger_type), POINTER :: logger
416 : CHARACTER(len=default_path_length) :: save_name
417 : INTEGER :: k, n, field_unit
418 : REAL(kind=dp), DIMENSION(3) :: real_field
419 :
420 : ! Get logger
421 14 : logger => cp_get_default_logger()
422 : ! Get file name
423 14 : save_name = cp_print_key_generate_filename(logger, rtbse_env%field_section, extension=".dat", my_local=.FALSE.)
424 14 : IF (file_exists(save_name)) THEN
425 : CALL open_file(save_name, file_status="OLD", file_form="FORMATTED", file_action="READ", &
426 0 : unit_number=field_unit)
427 : ! Skip the first line - it contains headers
428 0 : READ (field_unit, '()')
429 0 : DO k = rtbse_env%sim_start_orig, rtbse_env%sim_start
430 0 : n = k - rtbse_env%sim_start_orig + 1
431 0 : READ (field_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') rtbse_env%time_trace(n), &
432 0 : real_field(1), real_field(2), real_field(3)
433 0 : rtbse_env%field_trace(:, n) = CMPLX(real_field(:), 0.0, kind=dp)
434 : ! Set the time units back to atomic units
435 0 : rtbse_env%time_trace(n) = rtbse_env%time_trace(n)/femtoseconds
436 : END DO
437 0 : CALL close_file(field_unit)
438 14 : ELSE IF (.NOT. rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. &
439 : rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
440 2 : CPWARN("Restart without RT field file - unknown field trace set to zero.")
441 : END IF
442 14 : END SUBROUTINE read_field
443 :
444 : ! **************************************************************************************************
445 : !> \brief Outputs the expectation value of moments from a given density matrix
446 : !> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3)
447 : !> \param rtbse_env Entry point - gwbse environment
448 : !> \param rho Density matrix in AO basis
449 : !> \param rtp_section RTP section of the input parameters, where moments destination may be present
450 : ! **************************************************************************************************
451 1562 : SUBROUTINE output_moments(rtbse_env, rho)
452 : TYPE(rtbse_env_type) :: rtbse_env
453 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
454 : INTEGER :: i, j, n
455 : REAL(kind=dp), DIMENSION(3) :: moments_re
456 1562 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: ws
457 :
458 1562 : n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1
459 :
460 : ! In linearized RT-BSE rho and moments are MO-active sized; otherwise AO sized.
461 1562 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) THEN
462 874 : ws => rtbse_env%real_workspace_mo
463 : ELSE
464 688 : ws => rtbse_env%real_workspace
465 : END IF
466 :
467 3250 : DO j = 1, rtbse_env%n_spin
468 : ! Need to transpose due to the definition of trace function
469 1688 : CALL cp_cfm_to_fm(msource=rho(j), mtargetr=ws(2))
470 6752 : DO i = 1, 3
471 : ! Moments should be symmetric, test without transopose?
472 5064 : CALL cp_fm_transpose(rtbse_env%moments(i, j), ws(1))
473 5064 : CALL cp_fm_trace(ws(1), ws(2), moments_re(i))
474 : ! Scale by spin degeneracy and electron charge
475 5064 : moments_re(i) = -moments_re(i)*rtbse_env%spin_degeneracy
476 6752 : rtbse_env%moments_trace(j, i, n) = CMPLX(moments_re(i), 0.0, kind=dp)
477 : END DO
478 : ! Same for imaginary part
479 1688 : CALL cp_cfm_to_fm(msource=rho(j), mtargeti=ws(2))
480 8314 : DO i = 1, 3
481 5064 : CALL cp_fm_transpose(rtbse_env%moments(i, j), ws(1))
482 5064 : CALL cp_fm_trace(ws(1), ws(2), moments_re(i))
483 : ! Scale by spin degeneracy and electron charge
484 5064 : moments_re(i) = -moments_re(i)*rtbse_env%spin_degeneracy
485 6752 : rtbse_env%moments_trace(j, i, n) = rtbse_env%moments_trace(j, i, n) + CMPLX(0.0, moments_re(i), kind=dp)
486 : END DO
487 : END DO
488 : ! Output to the file
489 : CALL print_moments(rtbse_env%moments_section, rtbse_env%unit_nr, rtbse_env%moments_trace(:, :, n), &
490 1562 : rtbse_env%sim_time, .TRUE., append_opt=(rtbse_env%sim_step /= rtbse_env%sim_start_orig))
491 1562 : END SUBROUTINE output_moments
492 : ! **************************************************************************************************
493 : !> \brief Outputs the restart info (last finished iteration step) + restard density matrix
494 : !> \param restart_section Print key section for the restart files
495 : !> \param rho Density matrix in AO basis
496 : !> \param time_index Time index to be written into the info file
497 : ! **************************************************************************************************
498 678 : SUBROUTINE output_restart(rtbse_env, rho, time_index)
499 : TYPE(rtbse_env_type), POINTER :: rtbse_env
500 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
501 : INTEGER :: time_index
502 678 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: workspace
503 : CHARACTER(len=17), DIMENSION(4) :: file_labels
504 : TYPE(cp_logger_type), POINTER :: logger
505 : INTEGER :: rho_unit_nr, i
506 :
507 : ! Default labels distinguishing up to two spin species and real/imaginary parts
508 678 : file_labels(1) = "_SPIN_A_RE.matrix"
509 678 : file_labels(2) = "_SPIN_A_IM.matrix"
510 678 : file_labels(3) = "_SPIN_B_RE.matrix"
511 678 : file_labels(4) = "_SPIN_B_IM.matrix"
512 :
513 1356 : logger => cp_get_default_logger()
514 :
515 678 : workspace => rtbse_env%real_workspace
516 :
517 1356 : DO i = 1, rtbse_env%n_spin
518 678 : CALL cp_cfm_to_fm(rho(i), workspace(1), workspace(2))
519 : ! Real part
520 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i - 1), &
521 678 : file_form="UNFORMATTED", file_position="REWIND")
522 678 : CALL cp_fm_write_unformatted(workspace(1), rho_unit_nr)
523 678 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
524 : ! Imag part
525 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i), &
526 678 : file_form="UNFORMATTED", file_position="REWIND")
527 678 : CALL cp_fm_write_unformatted(workspace(2), rho_unit_nr)
528 678 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
529 : ! Info
530 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=".info", &
531 678 : file_form="UNFORMATTED", file_position="REWIND")
532 678 : IF (rho_unit_nr > 0) WRITE (rho_unit_nr) time_index
533 1356 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
534 : END DO
535 678 : END SUBROUTINE output_restart
536 : ! **************************************************************************************************
537 : !> \brief Reads the density matrix from restart files and updates the starting time
538 : !> \param restart_section Print key section for the restart files
539 : !> \param rho Density matrix in AO basis
540 : !> \param time_index Time index to be written into the info file
541 : ! **************************************************************************************************
542 6 : SUBROUTINE read_restart(rtbse_env)
543 : TYPE(rtbse_env_type), POINTER :: rtbse_env
544 : TYPE(cp_logger_type), POINTER :: logger
545 : CHARACTER(len=default_path_length) :: save_name, save_name_2
546 : INTEGER :: rho_unit_nr, j
547 : CHARACTER(len=17), DIMENSION(4) :: file_labels
548 :
549 : ! This allows the delta kick and output of moment at time 0 in all cases
550 : ! except the case when both imaginary and real parts of the density are read
551 6 : rtbse_env%restart_extracted = .FALSE.
552 6 : logger => cp_get_default_logger()
553 : ! Start by probing/loading info file
554 6 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".info", my_local=.FALSE.)
555 6 : IF (file_exists(save_name)) THEN
556 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
557 4 : unit_number=rho_unit_nr)
558 4 : READ (rho_unit_nr) rtbse_env%sim_start
559 4 : CALL close_file(rho_unit_nr)
560 6 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A31,I25,A24)') " RTBSE| Starting from timestep ", &
561 4 : rtbse_env%sim_start, ", delta kick NOT applied"
562 : ELSE
563 2 : CPWARN("Restart required but no info file found - starting from sim_step given in input")
564 : END IF
565 :
566 : ! Default labels distinguishing up to two spin species and real/imaginary parts
567 6 : file_labels(1) = "_SPIN_A_RE.matrix"
568 6 : file_labels(2) = "_SPIN_A_IM.matrix"
569 6 : file_labels(3) = "_SPIN_B_RE.matrix"
570 6 : file_labels(4) = "_SPIN_B_IM.matrix"
571 12 : DO j = 1, rtbse_env%n_spin
572 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
573 6 : extension=file_labels(2*j - 1), my_local=.FALSE.)
574 : save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
575 6 : extension=file_labels(2*j), my_local=.FALSE.)
576 12 : IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN
577 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
578 4 : unit_number=rho_unit_nr)
579 4 : CALL cp_fm_read_unformatted(rtbse_env%real_workspace(1), rho_unit_nr)
580 4 : CALL close_file(rho_unit_nr)
581 : CALL open_file(save_name_2, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
582 4 : unit_number=rho_unit_nr)
583 4 : CALL cp_fm_read_unformatted(rtbse_env%real_workspace(2), rho_unit_nr)
584 4 : CALL close_file(rho_unit_nr)
585 : CALL cp_fm_to_cfm(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), &
586 4 : rtbse_env%rho(j))
587 4 : rtbse_env%restart_extracted = .TRUE.
588 : ELSE
589 2 : CPWARN("Restart without some restart matrices - starting from SCF density.")
590 : END IF
591 : END DO
592 6 : END SUBROUTINE read_restart
593 : ! **************************************************************************************************
594 : !> \brief Linearized RT-BSE restart writer. Writes the restart set: lab-frame MO-active density
595 : !> matrices, the .info step index (sim_step = steps completed = the resume step), the
596 : !> once-per-run C_active gauge reference, and the appended RESTART.trace record feeding the
597 : !> FT prefix on continuation. All indices derive from sim_step so the .info/.trace
598 : !> bookkeeping cannot drift apart. (Full RTBSE uses the upstream output_restart above.)
599 : !> \param rtbse_env RT-BSE environment
600 : !> \param rho Density matrix (lab frame) to store
601 : ! **************************************************************************************************
602 830 : SUBROUTINE output_restart_linearized(rtbse_env, rho)
603 : TYPE(rtbse_env_type), POINTER :: rtbse_env
604 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
605 830 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: workspace
606 : CHARACTER(len=17), DIMENSION(4) :: file_labels
607 : CHARACTER(len=16), DIMENSION(2) :: c_file_labels
608 : TYPE(cp_logger_type), POINTER :: logger
609 : INTEGER :: rho_unit_nr, i
610 :
611 : ! Default labels distinguishing up to two spin species and real/imaginary parts
612 830 : file_labels(1) = "_SPIN_A_RE.matrix"
613 830 : file_labels(2) = "_SPIN_A_IM.matrix"
614 830 : file_labels(3) = "_SPIN_B_RE.matrix"
615 830 : file_labels(4) = "_SPIN_B_IM.matrix"
616 830 : c_file_labels(1) = "_SPIN_A_C.matrix"
617 830 : c_file_labels(2) = "_SPIN_B_C.matrix"
618 :
619 1660 : logger => cp_get_default_logger()
620 :
621 : ! In linearized RT-BSE rho is MO-active sized; otherwise AO sized.
622 830 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) THEN
623 830 : workspace => rtbse_env%real_workspace_mo
624 : ELSE
625 0 : workspace => rtbse_env%real_workspace
626 : END IF
627 :
628 1780 : DO i = 1, rtbse_env%n_spin
629 950 : CALL cp_cfm_to_fm(rho(i), workspace(1), workspace(2))
630 : ! Real part
631 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i - 1), &
632 950 : file_form="UNFORMATTED", file_position="REWIND")
633 950 : CALL cp_fm_write_unformatted(workspace(1), rho_unit_nr)
634 950 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
635 : ! Imag part
636 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i), &
637 950 : file_form="UNFORMATTED", file_position="REWIND")
638 950 : CALL cp_fm_write_unformatted(workspace(2), rho_unit_nr)
639 950 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
640 : ! Info
641 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=".info", &
642 950 : file_form="UNFORMATTED", file_position="REWIND")
643 950 : IF (rho_unit_nr > 0) WRITE (rho_unit_nr) rtbse_env%sim_step
644 1780 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
645 : END DO
646 :
647 : ! Once-per-run C_active dump (linearized only): gauge reference for the restart basis bridge
648 830 : IF (ASSOCIATED(rtbse_env%real_workspace_mo) .AND. .NOT. rtbse_env%restart_C_written) THEN
649 124 : DO i = 1, rtbse_env%n_spin
650 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=c_file_labels(i), &
651 66 : file_form="UNFORMATTED", file_position="REWIND")
652 66 : CALL cp_fm_write_unformatted(rtbse_env%C_active(i), rho_unit_nr)
653 124 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
654 : END DO
655 58 : rtbse_env%restart_C_written = .TRUE.
656 : END IF
657 :
658 830 : CALL write_restart_trace(rtbse_env)
659 830 : END SUBROUTINE output_restart_linearized
660 : ! **************************************************************************************************
661 : !> \brief Appends the current observable-trace record to RESTART.trace. Records are keyed by the
662 : !> observable slot n_slot = sim_step - sim_start_orig + 1 - the SAME index output_moments/
663 : !> output_field write (both drivers bump sim_step inside the propagation call before the
664 : !> output calls run). On the first call of a run the file is rewritten from memory (header +
665 : !> records 1..n_slot), truncating leftovers from an aborted run; subsequent calls append one
666 : !> record. Ionode writes; layout matches read_restart_trace verbatim.
667 : !> \param rtbse_env RT-BSE environment
668 : ! **************************************************************************************************
669 830 : SUBROUTINE write_restart_trace(rtbse_env)
670 : TYPE(rtbse_env_type), POINTER :: rtbse_env
671 : TYPE(cp_logger_type), POINTER :: logger
672 : CHARACTER(len=default_path_length) :: save_name
673 : INTEGER :: trace_unit, n, n_eps, n_slot
674 :
675 830 : logger => cp_get_default_logger()
676 830 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".trace", my_local=.FALSE.)
677 :
678 830 : IF (rtbse_env%unit_nr > 0) THEN
679 415 : n_slot = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1
680 415 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) THEN
681 : ! scalar row count is safe: determine_active_mo_window takes the union window across
682 : ! spins, so eps_active is rectangular (mo_active, n_spin) by construction
683 415 : n_eps = SIZE(rtbse_env%eps_active, 1)
684 : ELSE
685 0 : n_eps = 0
686 : END IF
687 415 : IF (.NOT. rtbse_env%restart_trace_written) THEN
688 : CALL open_file(save_name, file_status="UNKNOWN", file_form="UNFORMATTED", file_action="WRITE", &
689 29 : file_position="REWIND", unit_number=trace_unit)
690 29 : WRITE (trace_unit) restart_trace_version
691 29 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) THEN
692 29 : WRITE (trace_unit) rtbse_env%n_spin, rtbse_env%mo_active, rtbse_env%n_ao
693 : ELSE
694 0 : WRITE (trace_unit) rtbse_env%n_spin, rtbse_env%n_ao, rtbse_env%n_ao
695 : END IF
696 29 : WRITE (trace_unit) rtbse_env%sim_dt
697 116 : WRITE (trace_unit) REAL(rtbse_env%dft_control%rtp_control%delta_pulse_direction, dp), &
698 58 : rtbse_env%dft_control%rtp_control%delta_pulse_scale
699 29 : WRITE (trace_unit) n_eps
700 452 : IF (n_eps > 0) WRITE (trace_unit) rtbse_env%eps_active
701 164 : DO n = 1, n_slot
702 540 : WRITE (trace_unit) n, rtbse_env%time_trace(n), rtbse_env%field_trace(:, n), &
703 1163 : rtbse_env%moments_trace(:, :, n)
704 : END DO
705 : ELSE
706 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="WRITE", &
707 386 : file_position="APPEND", unit_number=trace_unit)
708 1544 : WRITE (trace_unit) n_slot, rtbse_env%time_trace(n_slot), rtbse_env%field_trace(:, n_slot), &
709 3256 : rtbse_env%moments_trace(:, :, n_slot)
710 : END IF
711 415 : CALL close_file(trace_unit)
712 : END IF
713 830 : rtbse_env%restart_trace_written = .TRUE.
714 830 : END SUBROUTINE write_restart_trace
715 : ! **************************************************************************************************
716 : !> \brief Early phase of the restart read: the starting step index from the .info file, the original
717 : !> run's dt peeked from the RESTART.trace header (so ENFORCE_MAX_DT can inherit it), and an
718 : !> SCF_GUESS hygiene check. Runs BEFORE initialize_maximum_timestep; the trace prefix records
719 : !> are loaded separately by read_restart_trace once the trace arrays are sized.
720 : !> \param rtbse_env RT-BSE environment
721 : ! **************************************************************************************************
722 14 : SUBROUTINE read_restart_info(rtbse_env)
723 : TYPE(rtbse_env_type), POINTER :: rtbse_env
724 : TYPE(cp_logger_type), POINTER :: logger
725 : CHARACTER(len=default_path_length) :: save_name
726 : INTEGER :: info_unit, trace_unit, version
727 : TYPE(scf_control_type), POINTER :: scf_control
728 :
729 14 : logger => cp_get_default_logger()
730 :
731 14 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".info", my_local=.FALSE.)
732 14 : IF (file_exists(save_name)) THEN
733 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
734 14 : unit_number=info_unit)
735 14 : READ (info_unit) rtbse_env%sim_start
736 14 : CALL close_file(info_unit)
737 21 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A31,I25,A24)') " RTBSE| Starting from timestep ", &
738 14 : rtbse_env%sim_start, ", delta kick NOT applied"
739 : ELSE
740 0 : CPWARN("Restart required but no info file found - starting from sim_step given in input")
741 : END IF
742 :
743 : ! Peek the original run's dt from the trace header (version, dims record, dt) so ENFORCE_MAX_DT
744 : ! can inherit it rather than recompute a window-dependent dt. The dims record is skipped here;
745 : ! read_restart_trace re-reads the full header and validates it once the trace arrays are sized.
746 14 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".trace", my_local=.FALSE.)
747 14 : IF (file_exists(save_name)) THEN
748 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
749 14 : unit_number=trace_unit)
750 14 : READ (trace_unit) version
751 14 : IF (version == restart_trace_version) THEN
752 14 : READ (trace_unit)
753 14 : READ (trace_unit) rtbse_env%sim_dt_restart
754 : END IF
755 14 : CALL close_file(trace_unit)
756 : END IF
757 :
758 : ! Hygiene nudge only - correctness is protected by the restart basis bridge (linearized path)
759 14 : NULLIFY (scf_control)
760 14 : CALL get_qs_env(rtbse_env%qs_env, scf_control=scf_control)
761 14 : IF (scf_control%density_guess /= restart_guess) THEN
762 0 : CPWARN("RT_RESTART without SCF_GUESS RESTART - SCF may reconverge to a gauge-rotated MO basis.")
763 : END IF
764 14 : END SUBROUTINE read_restart_info
765 : ! **************************************************************************************************
766 : !> \brief Reads the RESTART.trace prefix (records 1..sim_start) into the in-memory moment/field/
767 : !> time traces so the continuation FT covers the full history. Header guards: version,
768 : !> n_spin, dims, dt (abort); kick params (warn); eps_active > 0.1 meV (warn). All ranks read
769 : !> (the traces are replicated). Record layout matches write_restart_trace verbatim.
770 : !> \param rtbse_env RT-BSE environment
771 : ! **************************************************************************************************
772 14 : SUBROUTINE read_restart_trace(rtbse_env)
773 : TYPE(rtbse_env_type), POINTER :: rtbse_env
774 : TYPE(cp_logger_type), POINTER :: logger
775 : CHARACTER(len=default_path_length) :: save_name, err_msg
776 : INTEGER :: trace_unit, version, n_spin_file, &
777 : n_basis_file, n_ao_file, n_eps, n, &
778 : n_read_max, ios
779 : REAL(kind=dp) :: dt_file, kick_scale_file, t_rec
780 : REAL(kind=dp), DIMENSION(3) :: kick_dir_file
781 14 : REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: eps_file
782 : COMPLEX(kind=dp), DIMENSION(3) :: field_rec
783 14 : COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: mom_rec
784 :
785 14 : logger => cp_get_default_logger()
786 14 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".trace", my_local=.FALSE.)
787 : ! No trace + real prior history: the FT would run on a zero prefix - either a guaranteed
788 : ! multi_fft abort AFTER the full propagation (default FT%START_TIME=0, delta_t=0) or a
789 : ! silently wrong tail-only spectrum (START_TIME>0). Fail fast instead.
790 14 : IF (.NOT. file_exists(save_name)) THEN
791 0 : IF (rtbse_env%sim_start > 0) THEN
792 : CALL cp_abort(__LOCATION__, &
793 : "RT_RESTART without RESTART.trace - the continuation FT would miss the pre-restart "// &
794 0 : "history. Restore the original run's RESTART.trace next to the density restart files.")
795 : END IF
796 : RETURN
797 : END IF
798 :
799 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
800 14 : unit_number=trace_unit)
801 14 : READ (trace_unit) version
802 14 : IF (version /= restart_trace_version) THEN
803 0 : WRITE (err_msg, '(A,I0,A,I0)') "RESTART.trace: format version ", version, &
804 0 : " does not match this binary's version ", restart_trace_version
805 0 : CALL cp_abort(__LOCATION__, TRIM(err_msg))
806 : END IF
807 14 : READ (trace_unit) n_spin_file, n_basis_file, n_ao_file
808 14 : IF (n_spin_file /= rtbse_env%n_spin) CPABORT("RESTART.trace: n_spin mismatch")
809 14 : IF (n_ao_file /= rtbse_env%n_ao) CPABORT("RESTART.trace: n_ao mismatch")
810 : ! v1 is strict same-window: the linearized active-MO count must match (D6)
811 14 : IF (ASSOCIATED(rtbse_env%real_workspace_mo) .AND. n_basis_file /= rtbse_env%mo_active) THEN
812 : CALL cp_abort(__LOCATION__, &
813 0 : "RESTART.trace: active-MO count differs from the original run - same-window continuation only")
814 : END IF
815 14 : READ (trace_unit) dt_file
816 : ! ENFORCE_MAX_DT inherits dt_file (read_restart_info -> initialize_maximum_timestep), so this
817 : ! fires only on the manual path (ENFORCE off + user dt /= original); name the exact fix.
818 14 : IF (ABS(dt_file - rtbse_env%sim_dt) > 1.0e-12_dp*MAX(1.0_dp, ABS(rtbse_env%sim_dt))) THEN
819 : WRITE (err_msg, '(A,ES16.9,A)') &
820 : "RESTART.trace: TIMESTEP differs from the original run - continuation undefined. "// &
821 0 : "Set MD%TIMESTEP [fs] ", dt_file*femtoseconds, &
822 0 : " or enable RTBSE%ENFORCE_MAX_DT to inherit it automatically."
823 0 : CALL cp_abort(__LOCATION__, TRIM(err_msg))
824 : END IF
825 14 : READ (trace_unit) kick_dir_file, kick_scale_file
826 56 : IF (MAXVAL(ABS(kick_dir_file - REAL(rtbse_env%dft_control%rtp_control%delta_pulse_direction, dp))) > 1.0e-12_dp .OR. &
827 : ABS(kick_scale_file - rtbse_env%dft_control%rtp_control%delta_pulse_scale) > 1.0e-12_dp) THEN
828 : CALL cp_warn(__LOCATION__, &
829 0 : "RESTART.trace: delta-kick parameters differ from the original run - FT normalization inconsistent.")
830 : END IF
831 14 : READ (trace_unit) n_eps
832 14 : IF (n_eps > 0) THEN
833 56 : ALLOCATE (eps_file(n_eps, rtbse_env%n_spin))
834 14 : READ (trace_unit) eps_file
835 14 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) THEN
836 14 : IF (n_eps /= SIZE(rtbse_env%eps_active, 1)) THEN
837 0 : CALL cp_abort(__LOCATION__, "RESTART.trace: active-window size mismatch")
838 : END IF
839 : ! eps_active is not populated yet (built later in the Hamiltonian init); stash for the
840 : ! post-Hamiltonian consistency check in check_restart_eps_consistency.
841 56 : ALLOCATE (rtbse_env%eps_active_restart(n_eps, rtbse_env%n_spin))
842 222 : rtbse_env%eps_active_restart(:, :) = eps_file
843 : END IF
844 14 : DEALLOCATE (eps_file)
845 : END IF
846 :
847 42 : ALLOCATE (mom_rec(rtbse_env%n_spin, 3))
848 14 : n_read_max = 0
849 168 : DO
850 182 : READ (trace_unit, IOSTAT=ios) n, t_rec, field_rec, mom_rec
851 182 : IF (ios /= 0) EXIT
852 : ! Fill through slot sim_start+1: both drivers bump sim_step inside the propagation call
853 : ! (etrs_scf_loop / solve_rk4_timestep), so the live loop's first step writes slot
854 : ! sim_start+2 and slot sim_start+1 must be reloaded here.
855 168 : IF (n > rtbse_env%sim_start + 1) EXIT
856 168 : IF (n > SIZE(rtbse_env%time_trace)) THEN
857 0 : CPABORT("RESTART.trace: record index exceeds trace size - increase MOTION%MD%STEPS")
858 : END IF
859 168 : rtbse_env%time_trace(n) = t_rec
860 672 : rtbse_env%field_trace(:, n) = field_rec
861 1242 : rtbse_env%moments_trace(:, :, n) = mom_rec
862 182 : n_read_max = MAX(n_read_max, n)
863 : END DO
864 14 : DEALLOCATE (mom_rec)
865 14 : CALL close_file(trace_unit)
866 14 : IF (n_read_max < rtbse_env%sim_start) THEN
867 0 : CPWARN("RESTART.trace: fewer records than restart step - trace prefix incomplete.")
868 : END IF
869 14 : END SUBROUTINE read_restart_trace
870 : ! **************************************************************************************************
871 : !> \brief Compares the original run's active eigenvalues (stashed by read_restart_trace) against
872 : !> the recomputed eps_active, once the Hamiltonian is built. Prints the max deviation and
873 : !> warns above 0.1 meV (GW analytic continuation gives run-to-run QP noise above FP, so the
874 : !> threshold is deliberately loose). Frees the stash. No-op if nothing was stashed.
875 : !> \param rtbse_env RT-BSE environment
876 : ! **************************************************************************************************
877 58 : SUBROUTINE check_restart_eps_consistency(rtbse_env)
878 : TYPE(rtbse_env_type), POINTER :: rtbse_env
879 : REAL(kind=dp) :: eps_dev
880 :
881 58 : IF (.NOT. ASSOCIATED(rtbse_env%eps_active_restart)) RETURN
882 :
883 222 : eps_dev = MAXVAL(ABS(rtbse_env%eps_active_restart - rtbse_env%eps_active))
884 14 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A,ES12.3,A)') &
885 7 : " RTBSE| Restart eps_active max deviation vs original run ", eps_dev, " Ha"
886 : ! 0.1 meV = 1.0e-4 eV, converted to Ha via evolt (Ha -> eV factor)
887 14 : IF (eps_dev > 1.0e-4_dp/evolt) THEN
888 : CALL cp_warn(__LOCATION__, &
889 : "RESTART.trace: active eigenvalues deviate beyond 0.1 meV - "// &
890 0 : "Hamiltonian changed; continuation is physically inconsistent.")
891 : END IF
892 :
893 14 : DEALLOCATE (rtbse_env%eps_active_restart)
894 14 : NULLIFY (rtbse_env%eps_active_restart)
895 : END SUBROUTINE check_restart_eps_consistency
896 : ! **************************************************************************************************
897 : !> \brief Late phase of the restart read: overwrites rho from the lab-frame restart matrices and
898 : !> sets restart_extracted. MO-active for the linearized path, AO otherwise
899 : !> (cp_fm_read_unformatted aborts on a size mismatch, guarding a changed active window).
900 : !> \param rtbse_env RT-BSE environment
901 : ! **************************************************************************************************
902 14 : SUBROUTINE read_restart_density(rtbse_env)
903 : TYPE(rtbse_env_type), POINTER :: rtbse_env
904 : TYPE(cp_logger_type), POINTER :: logger
905 : CHARACTER(len=default_path_length) :: save_name, save_name_2
906 : INTEGER :: rho_unit_nr, j
907 : CHARACTER(len=17), DIMENSION(4) :: file_labels
908 14 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: ws
909 :
910 : ! This allows the delta kick and output of moment at time 0 in all cases
911 : ! except the case when both imaginary and real parts of the density are read
912 14 : rtbse_env%restart_extracted = .FALSE.
913 28 : logger => cp_get_default_logger()
914 :
915 14 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) THEN
916 14 : ws => rtbse_env%real_workspace_mo
917 : ELSE
918 0 : ws => rtbse_env%real_workspace
919 : END IF
920 :
921 : ! Default labels distinguishing up to two spin species and real/imaginary parts
922 14 : file_labels(1) = "_SPIN_A_RE.matrix"
923 14 : file_labels(2) = "_SPIN_A_IM.matrix"
924 14 : file_labels(3) = "_SPIN_B_RE.matrix"
925 14 : file_labels(4) = "_SPIN_B_IM.matrix"
926 30 : DO j = 1, rtbse_env%n_spin
927 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
928 16 : extension=file_labels(2*j - 1), my_local=.FALSE.)
929 : save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
930 16 : extension=file_labels(2*j), my_local=.FALSE.)
931 30 : IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN
932 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
933 16 : unit_number=rho_unit_nr)
934 16 : CALL cp_fm_read_unformatted(ws(1), rho_unit_nr)
935 16 : CALL close_file(rho_unit_nr)
936 : CALL open_file(save_name_2, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
937 16 : unit_number=rho_unit_nr)
938 16 : CALL cp_fm_read_unformatted(ws(2), rho_unit_nr)
939 16 : CALL close_file(rho_unit_nr)
940 : CALL cp_fm_to_cfm(ws(1), ws(2), &
941 16 : rtbse_env%rho(j))
942 16 : rtbse_env%restart_extracted = .TRUE.
943 : ELSE
944 0 : CPWARN("Restart without some restart matrices - starting from SCF density.")
945 : END IF
946 : END DO
947 14 : END SUBROUTINE read_restart_density
948 : ! **************************************************************************************************
949 : !> \brief Reads the previous run's C_active slabs (gauge reference for the restart basis bridge).
950 : !> \param rtbse_env RT-BSE environment
951 : !> \param C_old Caller-created fm array (n_spin) on fm_struct_ao_mo_active, filled on success
952 : !> \param found .TRUE. iff all per-spin C files were present and read
953 : ! **************************************************************************************************
954 14 : SUBROUTINE read_restart_C(rtbse_env, C_old, found)
955 : TYPE(rtbse_env_type), POINTER :: rtbse_env
956 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: C_old
957 : LOGICAL, INTENT(OUT) :: found
958 : TYPE(cp_logger_type), POINTER :: logger
959 : CHARACTER(len=default_path_length) :: save_name
960 : CHARACTER(len=16), DIMENSION(2) :: c_file_labels
961 : INTEGER :: c_unit, j
962 :
963 14 : c_file_labels(1) = "_SPIN_A_C.matrix"
964 14 : c_file_labels(2) = "_SPIN_B_C.matrix"
965 14 : logger => cp_get_default_logger()
966 14 : found = .FALSE.
967 30 : DO j = 1, rtbse_env%n_spin
968 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
969 16 : extension=c_file_labels(j), my_local=.FALSE.)
970 16 : IF (.NOT. file_exists(save_name)) THEN
971 0 : CPWARN("Restart without C_active file - assuming identical MO gauge (no basis bridge).")
972 0 : RETURN
973 : END IF
974 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
975 16 : unit_number=c_unit)
976 16 : CALL cp_fm_read_unformatted(C_old(j), c_unit)
977 30 : CALL close_file(c_unit)
978 : END DO
979 14 : found = .TRUE.
980 : END SUBROUTINE read_restart_C
981 : END MODULE rt_bse_io
|
