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