Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Input/output from the propagation via RT-BSE method.
10 : !> \author Stepan Marek (08.24)
11 : ! **************************************************************************************************
12 :
13 : MODULE rt_bse_io
14 : USE cp_fm_types, ONLY: cp_fm_type, &
15 : cp_fm_p_type, &
16 : cp_fm_create, &
17 : cp_fm_release, &
18 : cp_fm_read_unformatted, &
19 : cp_fm_write_unformatted, &
20 : cp_fm_write_formatted
21 : USE cp_cfm_types, ONLY: cp_cfm_type, &
22 : cp_cfm_p_type, &
23 : cp_fm_to_cfm, &
24 : cp_cfm_to_fm
25 : USE kinds, ONLY: dp, &
26 : default_path_length
27 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace, &
28 : cp_fm_transpose, &
29 : cp_fm_norm
30 : USE parallel_gemm_api, ONLY: parallel_gemm
31 : USE cp_log_handling, ONLY: cp_logger_type, &
32 : cp_get_default_logger
33 : USE cp_output_handling, ONLY: cp_print_key_unit_nr, &
34 : cp_print_key_finished_output, &
35 : cp_print_key_generate_filename, &
36 : low_print_level, &
37 : medium_print_level
38 : USE input_section_types, ONLY: section_vals_type
39 : USE rt_bse_types, ONLY: rtbse_env_type, &
40 : multiply_cfm_fm, &
41 : multiply_fm_cfm
42 : USE cp_files, ONLY: open_file, &
43 : file_exists, &
44 : close_file
45 : USE input_constants, ONLY: do_exact, &
46 : do_bch, &
47 : rtp_bse_ham_g0w0, &
48 : rtp_bse_ham_ks, &
49 : use_rt_restart
50 : USE physcon, ONLY: femtoseconds, &
51 : evolt
52 : USE mathconstants, ONLY: twopi
53 :
54 : #include "../base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse_io"
61 :
62 : #:include "rt_bse_macros.fypp"
63 :
64 : PUBLIC :: output_moments, &
65 : read_moments, &
66 : output_moments_ft, &
67 : output_polarizability, &
68 : output_field, &
69 : read_field, &
70 : output_mos_contravariant, &
71 : output_mos_covariant, &
72 : output_restart, &
73 : read_restart, &
74 : print_etrs_info_header, &
75 : print_etrs_info, &
76 : print_timestep_info, &
77 : print_rtbse_header_info
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief Writes the header and basic info to the standard output
83 : !> \param rtbse_env Entry point - rtbse environment
84 : ! **************************************************************************************************
85 6 : SUBROUTINE print_rtbse_header_info(rtbse_env)
86 : TYPE(rtbse_env_type) :: rtbse_env
87 :
88 6 : WRITE (rtbse_env%unit_nr, *) ''
89 : WRITE (rtbse_env%unit_nr, '(A)') ' /-----------------------------------------------'// &
90 6 : '------------------------------\'
91 : WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
92 6 : ' |'
93 : WRITE (rtbse_env%unit_nr, '(A)') ' | Real Time Bethe-Salpeter Propagation'// &
94 6 : ' |'
95 : WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
96 6 : ' |'
97 : WRITE (rtbse_env%unit_nr, '(A)') ' \-----------------------------------------------'// &
98 6 : '------------------------------/'
99 6 : WRITE (rtbse_env%unit_nr, *) ''
100 :
101 : ! Methods used
102 6 : WRITE (rtbse_env%unit_nr, '(A19)', advance="no") ' Exponential method'
103 12 : SELECT CASE (rtbse_env%mat_exp_method)
104 : CASE (do_bch)
105 6 : WRITE (rtbse_env%unit_nr, '(A61)') 'BCH'
106 : CASE (do_exact)
107 6 : WRITE (rtbse_env%unit_nr, '(A61)') 'EXACT'
108 : END SELECT
109 :
110 6 : WRITE (rtbse_env%unit_nr, '(A22)', advance="no") ' Reference Hamiltonian'
111 12 : SELECT CASE (rtbse_env%ham_reference_type)
112 : CASE (rtp_bse_ham_g0w0)
113 6 : WRITE (rtbse_env%unit_nr, '(A58)') 'G0W0'
114 : CASE (rtp_bse_ham_ks)
115 6 : WRITE (rtbse_env%unit_nr, '(A58)') 'Kohn-Sham'
116 : END SELECT
117 :
118 6 : WRITE (rtbse_env%unit_nr, '(A18,L62)') ' Apply delta pulse', &
119 12 : rtbse_env%dft_control%rtp_control%apply_delta_pulse
120 :
121 6 : WRITE (rtbse_env%unit_nr, '(A)') ''
122 :
123 6 : END SUBROUTINE
124 :
125 : ! **************************************************************************************************
126 : !> \brief Writes the update after single etrs iteration - only for log level > medium
127 : !> \param rtbse_env Entry point - rtbse environment
128 : ! **************************************************************************************************
129 1824 : SUBROUTINE print_etrs_info(rtbse_env, step, metric)
130 : TYPE(rtbse_env_type) :: rtbse_env
131 : INTEGER :: step
132 : REAL(kind=dp) :: metric
133 : TYPE(cp_logger_type), POINTER :: logger
134 :
135 1824 : logger => cp_get_default_logger()
136 :
137 1824 : IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN
138 0 : WRITE (rtbse_env%unit_nr, '(A7,I5, E20.8E3)') ' RTBSE|', step, metric
139 : END IF
140 :
141 1824 : END SUBROUTINE
142 : ! **************************************************************************************************
143 : !> \brief Writes the header for the etrs iteration updates - only for log level > medium
144 : !> \param rtbse_env Entry point - rtbse environment
145 : ! **************************************************************************************************
146 608 : SUBROUTINE print_etrs_info_header(rtbse_env)
147 : TYPE(rtbse_env_type) :: rtbse_env
148 : TYPE(cp_logger_type), POINTER :: logger
149 :
150 608 : logger => cp_get_default_logger()
151 :
152 608 : IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN
153 0 : WRITE (rtbse_env%unit_nr, '(A13, A20)') ' RTBSE| Iter.', 'Convergence'
154 : END IF
155 :
156 608 : END SUBROUTINE
157 : ! **************************************************************************************************
158 : !> \brief Writes the header for the etrs iteration updates - only for log level > medium
159 : !> \param rtbse_env Entry point - rtbse environment
160 : ! **************************************************************************************************
161 608 : SUBROUTINE print_timestep_info(rtbse_env, step, convergence, electron_num_re, etrs_num)
162 : TYPE(rtbse_env_type) :: rtbse_env
163 : INTEGER :: step
164 : REAL(kind=dp) :: convergence
165 : REAL(kind=dp) :: electron_num_re
166 : INTEGER :: etrs_num
167 : TYPE(cp_logger_type), POINTER :: logger
168 :
169 608 : logger => cp_get_default_logger()
170 :
171 608 : IF (logger%iter_info%print_level > low_print_level .AND. rtbse_env%unit_nr > 0) THEN
172 304 : WRITE (rtbse_env%unit_nr, '(A23,A20,A20,A17)') " RTBSE| Simulation step", "Convergence", &
173 608 : "Electron number", "ETRS Iterations"
174 304 : WRITE (rtbse_env%unit_nr, '(A7,I16,E20.8E3,E20.8E3,I17)') ' RTBSE|', step, convergence, &
175 608 : electron_num_re, etrs_num
176 : END IF
177 :
178 608 : END SUBROUTINE
179 :
180 : ! **************************************************************************************************
181 : !> \brief Outputs the matrix in MO basis for matrix coefficients corresponding to contravariant
182 : !> operator, i.e. density matrix
183 : !> \param rtbse_env Entry point - gwbse environment
184 : !> \param rho Density matrix in AO basis
185 : !> \param rtp_section RTP input section
186 : ! **************************************************************************************************
187 608 : SUBROUTINE output_mos_contravariant(rtbse_env, rho, print_key_section)
188 : TYPE(rtbse_env_type) :: rtbse_env
189 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
190 : TYPE(section_vals_type), POINTER :: print_key_section
191 : TYPE(cp_logger_type), POINTER :: logger
192 : INTEGER :: j, rho_unit_re, rho_unit_im
193 : CHARACTER(len=14), DIMENSION(4) :: file_labels
194 :
195 608 : file_labels(1) = "_SPIN_A_RE.dat"
196 608 : file_labels(2) = "_SPIN_A_IM.dat"
197 608 : file_labels(3) = "_SPIN_B_RE.dat"
198 608 : file_labels(4) = "_SPIN_B_IM.dat"
199 608 : logger => cp_get_default_logger()
200 : ! Start by multiplying the current density by MOS
201 1216 : DO j = 1, rtbse_env%n_spin
202 608 : rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1))
203 608 : rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j))
204 : ! Transform the density matrix into molecular orbitals basis and print it out
205 : ! S * rho
206 : CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
207 : 1.0_dp, rtbse_env%S_fm, rho(j), &
208 608 : 0.0_dp, rtbse_env%rho_workspace(1))
209 : ! C^T * S * rho
210 : CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
211 : 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), rtbse_env%rho_workspace(1), &
212 608 : 0.0_dp, rtbse_env%rho_workspace(2))
213 : ! C^T * S * rho * S
214 : CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
215 : 1.0_dp, rtbse_env%rho_workspace(2), rtbse_env%S_fm, &
216 608 : 0.0_dp, rtbse_env%rho_workspace(1))
217 : ! C^T * S * rho * S * C
218 : CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
219 : 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), &
220 608 : 0.0_dp, rtbse_env%rho_workspace(2))
221 : ! Print real and imaginary parts separately
222 : CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), &
223 608 : rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
224 608 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re)
225 608 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im)
226 608 : CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section)
227 1216 : CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section)
228 : END DO
229 608 : END SUBROUTINE
230 : ! **************************************************************************************************
231 : !> \brief Outputs the matrix in MO basis for matrix components corresponding to covariant representation,
232 : !> i.e. the Hamiltonian matrix
233 : !> \param rtbse_env Entry point - gwbse environment
234 : !> \param cohsex cohsex matrix in AO basis, covariant representation
235 : !> \param rtp_section RTP input section
236 : ! **************************************************************************************************
237 0 : SUBROUTINE output_mos_covariant(rtbse_env, ham, print_key_section)
238 : TYPE(rtbse_env_type) :: rtbse_env
239 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham
240 : TYPE(section_vals_type), POINTER :: print_key_section
241 : TYPE(cp_logger_type), POINTER :: logger
242 : INTEGER :: j, rho_unit_re, rho_unit_im
243 : CHARACTER(len=21), DIMENSION(4) :: file_labels
244 :
245 0 : file_labels(1) = "_SPIN_A_RE.dat"
246 0 : file_labels(2) = "_SPIN_A_IM.dat"
247 0 : file_labels(3) = "_SPIN_B_RE.dat"
248 0 : file_labels(4) = "_SPIN_B_IM.dat"
249 0 : logger => cp_get_default_logger()
250 0 : DO j = 1, rtbse_env%n_spin
251 0 : rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1))
252 0 : rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j))
253 : ! C^T * cohsex
254 : CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
255 : 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), ham(j), &
256 0 : 0.0_dp, rtbse_env%rho_workspace(1))
257 : ! C^T * cohsex * C
258 : CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
259 : 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), &
260 0 : 0.0_dp, rtbse_env%rho_workspace(2))
261 : ! Print real and imaginary parts separately
262 : CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), &
263 0 : rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
264 0 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re)
265 0 : CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im)
266 0 : CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section)
267 0 : CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section)
268 : END DO
269 0 : END SUBROUTINE
270 : ! **************************************************************************************************
271 : !> \brief Prints the current field components into a file provided by input
272 : !> \param rtbse_env Entry point - gwbse environment
273 : !> \param rtp_section RTP input section
274 : ! **************************************************************************************************
275 616 : SUBROUTINE output_field(rtbse_env)
276 : TYPE(rtbse_env_type) :: rtbse_env
277 : TYPE(cp_logger_type), POINTER :: logger
278 : INTEGER :: field_unit, n, i
279 :
280 : ! Get logger
281 616 : logger => cp_get_default_logger()
282 : ! Get file descriptor
283 616 : field_unit = cp_print_key_unit_nr(logger, rtbse_env%field_section, extension=".dat")
284 : ! If the file descriptor is non-zero, output field
285 : ! TODO : Output also in SI
286 616 : IF (field_unit /= -1) THEN
287 0 : WRITE (field_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') rtbse_env%sim_time*femtoseconds, &
288 0 : rtbse_env%field(1), rtbse_env%field(2), rtbse_env%field(3)
289 : END IF
290 : ! Write the output to memory for FT
291 : ! Need the absolute index
292 616 : n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1
293 2464 : DO i = 1, 3
294 2464 : rtbse_env%field_trace(i)%series(n) = rtbse_env%field(i)
295 : END DO
296 616 : rtbse_env%time_trace%series(n) = rtbse_env%sim_time
297 616 : CALL cp_print_key_finished_output(field_unit, logger, rtbse_env%field_section)
298 :
299 616 : END SUBROUTINE
300 : ! **************************************************************************************************
301 : !> \brief Reads the field from the files provided by input - useful for the continuation run
302 : !> \param rtbse_env Entry point - gwbse environment
303 : !> \param rtp_section RTP input section
304 : ! **************************************************************************************************
305 12 : SUBROUTINE read_field(rtbse_env)
306 : TYPE(rtbse_env_type) :: rtbse_env
307 : TYPE(cp_logger_type), POINTER :: logger
308 : CHARACTER(len=default_path_length) :: save_name
309 : INTEGER :: k, n, field_unit
310 :
311 : ! Get logger
312 12 : logger => cp_get_default_logger()
313 : ! Get file name
314 12 : save_name = cp_print_key_generate_filename(logger, rtbse_env%field_section, extension=".dat", my_local=.FALSE.)
315 12 : IF (file_exists(save_name)) THEN
316 : CALL open_file(save_name, file_status="OLD", file_form="FORMATTED", file_action="READ", &
317 0 : unit_number=field_unit)
318 0 : DO k = rtbse_env%sim_start_orig, rtbse_env%sim_start
319 0 : n = k - rtbse_env%sim_start_orig + 1
320 0 : READ (field_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') rtbse_env%time_trace%series(n), &
321 0 : rtbse_env%field_trace(1)%series(n), rtbse_env%field_trace(2)%series(n), rtbse_env%field_trace(3)%series(n)
322 : ! Set the time units back to atomic units
323 0 : rtbse_env%time_trace%series(n) = rtbse_env%time_trace%series(n)/femtoseconds
324 : END DO
325 0 : CALL close_file(field_unit)
326 12 : ELSE IF (.NOT. rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. &
327 : rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
328 2 : CPWARN("Restart without RT field file - unknown field trace set to zero.")
329 : END IF
330 12 : END SUBROUTINE read_field
331 :
332 : ! **************************************************************************************************
333 : !> \brief Outputs the expectation value of moments from a given density matrix
334 : !> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3)
335 : !> \param rtbse_env Entry point - gwbse environment
336 : !> \param rho Density matrix in AO basis
337 : !> \param rtp_section RTP section of the input parameters, where moments destination may be present
338 : ! **************************************************************************************************
339 616 : SUBROUTINE output_moments(rtbse_env, rho)
340 : TYPE(rtbse_env_type) :: rtbse_env
341 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
342 : TYPE(cp_logger_type), POINTER :: logger
343 : INTEGER :: i, j, n, moments_unit_re, moments_unit_im
344 : CHARACTER(len=14), DIMENSION(4) :: file_labels
345 : REAL(kind=dp), DIMENSION(3) :: moments
346 :
347 : ! Start by getting the relevant file unit
348 616 : file_labels(1) = "_SPIN_A_RE.dat"
349 616 : file_labels(2) = "_SPIN_A_IM.dat"
350 616 : file_labels(3) = "_SPIN_B_RE.dat"
351 616 : file_labels(4) = "_SPIN_B_IM.dat"
352 616 : logger => cp_get_default_logger()
353 1232 : DO j = 1, rtbse_env%n_spin
354 616 : moments_unit_re = cp_print_key_unit_nr(logger, rtbse_env%moments_section, extension=file_labels(2*j - 1))
355 616 : moments_unit_im = cp_print_key_unit_nr(logger, rtbse_env%moments_section, extension=file_labels(2*j))
356 : ! If, for any reason, the file unit is not provided, skip to next cycle immediately
357 : ! TODO : Specify output units in config
358 : ! Need to transpose due to the definition of trace function
359 616 : CALL cp_cfm_to_fm(msource=rho(j), mtargetr=rtbse_env%real_workspace(2))
360 2464 : DO i = 1, 3
361 : ! Moments should be symmetric, test without transopose?
362 1848 : CALL cp_fm_transpose(rtbse_env%moments(i), rtbse_env%real_workspace(1))
363 1848 : CALL cp_fm_trace(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), moments(i))
364 : ! Scale by spin degeneracy and electron charge
365 2464 : moments(i) = -moments(i)*rtbse_env%spin_degeneracy
366 : END DO
367 : ! Output to the file
368 616 : IF (moments_unit_re > 0) THEN
369 : ! If outputting to standard output, also prepend identifying 'MOMENTS_TRACE|' string
370 256 : IF (rtbse_env%unit_nr == moments_unit_re) THEN
371 : WRITE (moments_unit_re, '(A18,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
372 256 : " MOMENTS_TRACE_RE|", rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
373 : ELSE
374 : WRITE (moments_unit_re, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
375 0 : rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
376 0 : CALL cp_print_key_finished_output(moments_unit_re, logger, rtbse_env%moments_section)
377 : END IF
378 : END IF
379 : ! Save to memory for FT - real part
380 616 : n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1
381 2464 : DO i = 1, 3
382 2464 : rtbse_env%moments_trace(i)%series(n) = CMPLX(moments(i), 0.0, kind=dp)
383 : END DO
384 : ! Same for imaginary part
385 616 : CALL cp_cfm_to_fm(msource=rho(j), mtargeti=rtbse_env%real_workspace(2))
386 2464 : DO i = 1, 3
387 1848 : CALL cp_fm_transpose(rtbse_env%moments(i), rtbse_env%real_workspace(1))
388 1848 : CALL cp_fm_trace(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), moments(i))
389 : ! Scale by spin degeneracy and electron charge
390 2464 : moments(i) = -moments(i)*rtbse_env%spin_degeneracy
391 : END DO
392 : ! Output to the file
393 616 : IF (moments_unit_im > 0) THEN
394 : ! If outputting to standard output, also prepend identifying 'MOMENTS_TRACE|' string
395 256 : IF (rtbse_env%unit_nr == moments_unit_im) THEN
396 : WRITE (moments_unit_im, '(A18,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
397 256 : " MOMENTS_TRACE_IM|", rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
398 : ELSE
399 : WRITE (moments_unit_im, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
400 0 : rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3)
401 : ! Close the files
402 0 : CALL cp_print_key_finished_output(moments_unit_im, logger, rtbse_env%moments_section)
403 : END IF
404 : END IF
405 : ! Save to memory for FT - imag part
406 3080 : DO i = 1, 3
407 2464 : rtbse_env%moments_trace(i)%series(n) = rtbse_env%moments_trace(i)%series(n) + CMPLX(0.0, moments(i), kind=dp)
408 : END DO
409 : END DO
410 616 : END SUBROUTINE
411 : ! **************************************************************************************************
412 : !> \brief Outputs the restart info (last finished iteration step) + restard density matrix
413 : !> \param restart_section Print key section for the restart files
414 : !> \param rho Density matrix in AO basis
415 : !> \param time_index Time index to be written into the info file
416 : ! **************************************************************************************************
417 608 : SUBROUTINE output_restart(rtbse_env, rho, time_index)
418 : TYPE(rtbse_env_type), POINTER :: rtbse_env
419 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
420 : INTEGER :: time_index
421 608 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: workspace
422 : CHARACTER(len=17), DIMENSION(4) :: file_labels
423 : TYPE(cp_logger_type), POINTER :: logger
424 : INTEGER :: rho_unit_nr, i
425 :
426 : ! Default labels distinguishing up to two spin species and real/imaginary parts
427 608 : file_labels(1) = "_SPIN_A_RE.matrix"
428 608 : file_labels(2) = "_SPIN_A_IM.matrix"
429 608 : file_labels(3) = "_SPIN_B_RE.matrix"
430 608 : file_labels(4) = "_SPIN_B_IM.matrix"
431 :
432 1216 : logger => cp_get_default_logger()
433 :
434 608 : workspace => rtbse_env%real_workspace
435 :
436 1216 : DO i = 1, rtbse_env%n_spin
437 608 : CALL cp_cfm_to_fm(rho(i), workspace(1), workspace(2))
438 : ! Real part
439 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i - 1), &
440 608 : file_form="UNFORMATTED", file_position="REWIND")
441 608 : CALL cp_fm_write_unformatted(workspace(1), rho_unit_nr)
442 608 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
443 : ! Imag part
444 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i), &
445 608 : file_form="UNFORMATTED", file_position="REWIND")
446 608 : CALL cp_fm_write_unformatted(workspace(2), rho_unit_nr)
447 608 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
448 : ! Info
449 : rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=".info", &
450 608 : file_form="UNFORMATTED", file_position="REWIND")
451 608 : IF (rho_unit_nr > 0) WRITE (rho_unit_nr) time_index
452 1216 : CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section)
453 : END DO
454 608 : END SUBROUTINE output_restart
455 : ! **************************************************************************************************
456 : !> \brief Reads the density matrix from restart files and updates the starting time
457 : !> \param restart_section Print key section for the restart files
458 : !> \param rho Density matrix in AO basis
459 : !> \param time_index Time index to be written into the info file
460 : ! **************************************************************************************************
461 4 : SUBROUTINE read_restart(rtbse_env)
462 : TYPE(rtbse_env_type), POINTER :: rtbse_env
463 : TYPE(cp_logger_type), POINTER :: logger
464 : CHARACTER(len=default_path_length) :: save_name, save_name_2
465 : INTEGER :: rho_unit_nr, j
466 : CHARACTER(len=17), DIMENSION(4) :: file_labels
467 :
468 : ! This allows the delta kick and output of moment at time 0 in all cases
469 : ! except the case when both imaginary and real parts of the density are read
470 4 : rtbse_env%restart_extracted = .FALSE.
471 4 : logger => cp_get_default_logger()
472 : ! Start by probing/loading info file
473 4 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".info", my_local=.FALSE.)
474 4 : IF (file_exists(save_name)) THEN
475 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
476 4 : unit_number=rho_unit_nr)
477 4 : READ (rho_unit_nr) rtbse_env%sim_start
478 4 : CALL close_file(rho_unit_nr)
479 6 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A31,I25,A24)') " RTBSE| Starting from timestep ", &
480 4 : rtbse_env%sim_start, ", delta kick NOT applied"
481 : ELSE
482 0 : CPWARN("Restart required but no info file found - starting from sim_step given in input")
483 : END IF
484 :
485 : ! Default labels distinguishing up to two spin species and real/imaginary parts
486 4 : file_labels(1) = "_SPIN_A_RE.matrix"
487 4 : file_labels(2) = "_SPIN_A_IM.matrix"
488 4 : file_labels(3) = "_SPIN_B_RE.matrix"
489 4 : file_labels(4) = "_SPIN_B_IM.matrix"
490 8 : DO j = 1, rtbse_env%n_spin
491 : save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
492 4 : extension=file_labels(2*j - 1), my_local=.FALSE.)
493 : save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%restart_section, &
494 4 : extension=file_labels(2*j), my_local=.FALSE.)
495 8 : IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN
496 : CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
497 4 : unit_number=rho_unit_nr)
498 4 : CALL cp_fm_read_unformatted(rtbse_env%real_workspace(1), rho_unit_nr)
499 4 : CALL close_file(rho_unit_nr)
500 : CALL open_file(save_name_2, file_status="OLD", file_form="UNFORMATTED", file_action="READ", &
501 4 : unit_number=rho_unit_nr)
502 4 : CALL cp_fm_read_unformatted(rtbse_env%real_workspace(2), rho_unit_nr)
503 4 : CALL close_file(rho_unit_nr)
504 : CALL cp_fm_to_cfm(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), &
505 4 : rtbse_env%rho(j))
506 4 : rtbse_env%restart_extracted = .TRUE.
507 : ELSE
508 0 : CPWARN("Restart without some restart matrices - starting from SCF density.")
509 : END IF
510 : END DO
511 4 : END SUBROUTINE read_restart
512 : ! **************************************************************************************************
513 : !> \brief Reads the moments and time traces from the save files
514 : !> \param rtbse_env GW-BSE environment (assumes consistent setup, i.e. a continuation run).
515 : !> Otherwise, the traces are set at zero
516 : ! **************************************************************************************************
517 12 : SUBROUTINE read_moments(rtbse_env)
518 : TYPE(rtbse_env_type), POINTER :: rtbse_env
519 : TYPE(cp_logger_type), POINTER :: logger
520 : CHARACTER(len=default_path_length) :: save_name, save_name_2
521 : INTEGER :: i, j, k, moments_unit_re, moments_unit_im, n
522 : CHARACTER(len=17), DIMENSION(4) :: file_labels
523 : REAL(kind=dp), DIMENSION(3) :: moments_re, moments_im
524 : REAL(kind=dp) :: timestamp
525 :
526 12 : logger => cp_get_default_logger()
527 :
528 : ! Read moments from the previous run
529 : ! Default labels distinguishing up to two spin species and real/imaginary parts
530 12 : file_labels(1) = "_SPIN_A_RE.dat"
531 12 : file_labels(2) = "_SPIN_A_IM.dat"
532 12 : file_labels(3) = "_SPIN_B_RE.dat"
533 12 : file_labels(4) = "_SPIN_B_IM.dat"
534 24 : DO j = 1, rtbse_env%n_spin
535 : save_name = cp_print_key_generate_filename(logger, rtbse_env%moments_section, &
536 12 : extension=file_labels(2*j - 1), my_local=.FALSE.)
537 : save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%moments_section, &
538 12 : extension=file_labels(2*j), my_local=.FALSE.)
539 24 : IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN
540 : CALL open_file(save_name, file_status="OLD", file_form="FORMATTED", file_action="READ", &
541 0 : unit_number=moments_unit_re)
542 : CALL open_file(save_name_2, file_status="OLD", file_form="FORMATTED", file_action="READ", &
543 0 : unit_number=moments_unit_im)
544 : ! Extra time step for the initial one
545 0 : DO k = rtbse_env%sim_start_orig, rtbse_env%sim_start
546 : ! Determine the absolute time step - offset in memory
547 0 : n = k - rtbse_env%sim_start_orig + 1
548 0 : READ (moments_unit_re, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') timestamp, &
549 0 : moments_re(1), moments_re(2), moments_re(3)
550 0 : READ (moments_unit_im, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') timestamp, &
551 0 : moments_im(1), moments_im(2), moments_im(3)
552 0 : DO i = 1, 3
553 0 : rtbse_env%moments_trace(i)%series(n) = CMPLX(moments_re(i), moments_im(i), kind=dp)
554 : END DO
555 : END DO
556 : ! Change back to atomic units in the trace
557 0 : rtbse_env%time_trace%series(:) = rtbse_env%time_trace%series(:)/femtoseconds
558 0 : CALL close_file(moments_unit_re)
559 0 : CALL close_file(moments_unit_im)
560 12 : ELSE IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
561 4 : CPWARN("Restart without previous moments - unknown moments trace set to zero.")
562 : END IF
563 : END DO
564 12 : END SUBROUTINE read_moments
565 : ! **************************************************************************************************
566 : !> \brief Outputs the Fourier transform of moments stored in the environment memory to the configured file
567 : !> \param rtbse_env GW-BSE environment
568 : ! **************************************************************************************************
569 12 : SUBROUTINE output_moments_ft(rtbse_env)
570 : TYPE(rtbse_env_type), POINTER :: rtbse_env
571 : TYPE(cp_logger_type), POINTER :: logger
572 : REAL(kind=dp), DIMENSION(:), POINTER :: omega_series, &
573 : ft_real_series, &
574 : ft_imag_series, &
575 : value_series
576 : ! Stores the data in ready for output format
577 : ! - first dimension is 6 - 1 - real part along x, 2 - imag part along x, 3 - real part along y, ...
578 12 : REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: ft_full_series
579 : INTEGER :: i, n, ft_unit
580 :
581 12 : logger => cp_get_default_logger()
582 :
583 12 : n = rtbse_env%sim_nsteps + 2
584 : NULLIFY (omega_series)
585 860 : ALLOCATE (omega_series(n), source=0.0_dp)
586 : NULLIFY (ft_real_series)
587 848 : ALLOCATE (ft_real_series(n), source=0.0_dp)
588 : NULLIFY (ft_imag_series)
589 848 : ALLOCATE (ft_imag_series(n), source=0.0_dp)
590 : NULLIFY (value_series)
591 848 : ALLOCATE (value_series(n), source=0.0_dp)
592 36 : ALLOCATE (ft_full_series(6, n))
593 : ! Carry out for each direction independently and real and imaginary parts also independently
594 48 : DO i = 1, 3
595 : ! Real part of the value first
596 2508 : value_series(:) = REAL(rtbse_env%moments_trace(i)%series(:))
597 : CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, &
598 36 : rtbse_env%ft_damping, rtbse_env%ft_start)
599 2508 : ft_full_series(2*i - 1, :) = ft_real_series(:)
600 2508 : ft_full_series(2*i, :) = ft_imag_series(:)
601 : ! Now imaginary part
602 2508 : value_series(:) = AIMAG(rtbse_env%moments_trace(i)%series(:))
603 : CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, &
604 36 : rtbse_env%ft_damping, rtbse_env%ft_start)
605 2508 : ft_full_series(2*i - 1, :) = ft_full_series(2*i - 1, :) - ft_imag_series
606 2520 : ft_full_series(2*i, :) = ft_full_series(2*i, :) + ft_real_series
607 : END DO
608 12 : DEALLOCATE (ft_real_series)
609 12 : DEALLOCATE (ft_imag_series)
610 12 : DEALLOCATE (value_series)
611 : ! Now, write these to file
612 : ft_unit = cp_print_key_unit_nr(logger, rtbse_env%ft_section, extension=".dat", &
613 12 : file_form="FORMATTED", file_position="REWIND")
614 12 : IF (ft_unit > 0) THEN
615 418 : DO i = 1, n
616 : WRITE (ft_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
617 412 : omega_series(i), ft_full_series(1, i), ft_full_series(2, i), &
618 412 : ft_full_series(3, i), ft_full_series(4, i), &
619 830 : ft_full_series(5, i), ft_full_series(6, i)
620 : END DO
621 : END IF
622 12 : CALL cp_print_key_finished_output(ft_unit, logger, rtbse_env%ft_section)
623 12 : DEALLOCATE (omega_series)
624 12 : DEALLOCATE (ft_full_series)
625 12 : END SUBROUTINE output_moments_ft
626 : ! **************************************************************************************************
627 : !> \brief Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega),
628 : !> where i and j are provided by the configuration. The tensor element is energy dependent and
629 : !> has real and imaginary parts
630 : !> \param rtbse_env GW-BSE environment
631 : ! **************************************************************************************************
632 12 : SUBROUTINE output_polarizability(rtbse_env)
633 : TYPE(rtbse_env_type), POINTER :: rtbse_env
634 : TYPE(cp_logger_type), POINTER :: logger
635 : REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_series, &
636 : ft_real_series, &
637 : ft_imag_series, &
638 : value_series
639 12 : COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: moment_series, &
640 12 : field_series, &
641 12 : polarizability_series
642 : INTEGER :: pol_unit, &
643 : i, n
644 :
645 12 : logger => cp_get_default_logger()
646 :
647 12 : n = rtbse_env%sim_nsteps + 2
648 : ! All allocations together, although could save some memory, if required by consequent deallocations
649 860 : ALLOCATE (omega_series(n), source=0.0_dp)
650 848 : ALLOCATE (ft_real_series(n), source=0.0_dp)
651 848 : ALLOCATE (ft_imag_series(n), source=0.0_dp)
652 848 : ALLOCATE (value_series(n), source=0.0_dp)
653 872 : ALLOCATE (moment_series(n), source=CMPLX(0.0, 0.0, kind=dp))
654 860 : ALLOCATE (field_series(n), source=CMPLX(0.0, 0.0, kind=dp))
655 860 : ALLOCATE (polarizability_series(n), source=CMPLX(0.0, 0.0, kind=dp))
656 :
657 : ! The moment ft
658 : ! Real part
659 836 : value_series(:) = REAL(rtbse_env%moments_trace(rtbse_env%pol_element(1))%series(:))
660 : CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, &
661 12 : rtbse_env%ft_damping, rtbse_env%ft_start)
662 836 : moment_series(:) = CMPLX(ft_real_series(:), ft_imag_series(:), kind=dp)
663 : ! Imaginary part
664 836 : value_series(:) = AIMAG(rtbse_env%moments_trace(rtbse_env%pol_element(1))%series(:))
665 : CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, &
666 12 : rtbse_env%ft_damping, rtbse_env%ft_start)
667 836 : moment_series(:) = moment_series(:) + CMPLX(-ft_imag_series(:), ft_real_series(:), kind=dp)
668 :
669 : ! Calculate the field transform - store it in ft_real_series
670 12 : IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse) THEN
671 : ! Only divide by constant magnitude in atomic units
672 418 : field_series(:) = CMPLX(rtbse_env%dft_control%rtp_control%delta_pulse_scale, 0.0, kind=dp)
673 : ELSE
674 : ! Calculate the transform of the field as well and divide by it
675 : ! The field FT is not damped - assume field is localised in time
676 : ! The field is strictly real
677 : CALL ft_simple(rtbse_env%time_trace%series, rtbse_env%field_trace(rtbse_env%pol_element(2))%series, &
678 6 : omega_series, ft_real_series, ft_imag_series, 0.0_dp, rtbse_env%ft_start)
679 : ! Regularization for the case when ft_series becomes identically zero - TODO : Set in config
680 418 : field_series(:) = CMPLX(ft_real_series(:), ft_imag_series(:) + 1.0e-20, kind=dp)
681 : END IF
682 : ! Divide to get the polarizability series
683 : ! Real part
684 836 : polarizability_series(:) = moment_series(:)/field_series(:)
685 :
686 : ! Change units to eV for energy
687 : ! use value_series for energy and moment_series for polarizability
688 836 : value_series(:) = omega_series(:)*evolt
689 : ! Use below for printing of field FT in case of problems
690 : ! PRINT *, "=== Field FT"
691 : ! DO i=1,n
692 : ! PRINT '(E20.8E3,E20.8E3,E20.8E3)', value_series(i), REAL(field_series(i)), AIMAG(field_series(i))
693 : ! END DO
694 : ! PRINT *, "=== Field FT"
695 : ! Print out the polarizability to a file
696 : pol_unit = cp_print_key_unit_nr(logger, rtbse_env%pol_section, extension=".dat", &
697 12 : file_form="FORMATTED", file_position="REWIND")
698 12 : IF (pol_unit > 0) THEN
699 6 : IF (pol_unit == rtbse_env%unit_nr) THEN
700 : WRITE (pol_unit, '(A16,A20,A22,A22)') &
701 1 : " POLARIZABILITY|", "Energy [eV]", "Real polarizability", "Imag polarizability"
702 53 : DO i = 1, n
703 : WRITE (pol_unit, '(A16,E20.8E3,E22.8E3,E22.8E3)') &
704 53 : " POLARIZABILITY|", value_series(i), REAL(polarizability_series(i)), AIMAG(polarizability_series(i))
705 : END DO
706 : ELSE
707 : WRITE (pol_unit, '(A1,A19,A20,A20,A20)') &
708 5 : "#", "omega [a.u.]", "Energy [eV]", "Real polarizability", "Imag polarizability"
709 365 : DO i = 1, n
710 : WRITE (pol_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
711 365 : omega_series(i), value_series(i), REAL(polarizability_series(i)), AIMAG(polarizability_series(i))
712 : END DO
713 5 : CALL cp_print_key_finished_output(pol_unit, logger, rtbse_env%ft_section)
714 : END IF
715 : END IF
716 :
717 12 : DEALLOCATE (value_series)
718 12 : DEALLOCATE (ft_real_series)
719 12 : DEALLOCATE (ft_imag_series)
720 :
721 12 : DEALLOCATE (field_series)
722 12 : DEALLOCATE (moment_series)
723 :
724 12 : DEALLOCATE (omega_series)
725 12 : DEALLOCATE (polarizability_series)
726 12 : END SUBROUTINE output_polarizability
727 : ! **************************************************************************************************
728 : !> \brief Naively calculates the Fourier transform - it is not the bottleneck of this calculation
729 : !> \param time_series Timestamps in atomic units of time
730 : !> \param value_series Values to be Fourier transformed - moments, field etc. So far only real.
731 : !> \param omega_series Array to be filled by sampled values of frequency
732 : !> \param result_series FT of the value series - real values (cosines)
733 : !> \param iresult_series FT of the value series - imaginary values (sines)
734 : !> \param damping_opt Supply custom exponential damping - default is 4.0/totalTime, i.e. ratio
735 : !> of last and first element in windowed value series is reduced by e^(-4)
736 : !> \param t0_opt Carry the FT only starting from certain time - allows for exclusion of trace before
737 : !> the pulse application etc.
738 : !> \author Stepan Marek
739 : !> \date 09.2024
740 : ! **************************************************************************************************
741 : ! So far only for defined one dimensional series
742 102 : SUBROUTINE ft_simple(time_series, value_series, omega_series, result_series, iresult_series, &
743 : damping_opt, t0_opt)
744 : REAL(kind=dp), DIMENSION(:) :: time_series
745 : REAL(kind=dp), DIMENSION(:) :: value_series
746 : REAL(kind=dp), DIMENSION(:) :: omega_series
747 : REAL(kind=dp), DIMENSION(:) :: result_series
748 : REAL(kind=dp), DIMENSION(:) :: iresult_series
749 : REAL(kind=dp), OPTIONAL :: damping_opt
750 : REAL(kind=dp), OPTIONAL :: t0_opt
751 : CHARACTER(len=*), PARAMETER :: routineN = "ft_simple"
752 : INTEGER :: N, i, j, t0_i, j_wrap, handle
753 : REAL(kind=dp) :: t0, delta_t, delta_omega, damping
754 :
755 102 : CALL timeset(routineN, handle)
756 :
757 102 : N = SIZE(time_series)
758 :
759 102 : t0_i = 1
760 102 : IF (PRESENT(t0_opt)) THEN
761 : ! Determine the index at which we start applying the damping
762 : ! Also the index around which we fold around
763 102 : DO i = 1, N
764 : ! Increase until we break or reach the end of the time series
765 102 : t0_i = i
766 102 : IF (time_series(i) >= t0_opt) THEN
767 : EXIT
768 : END IF
769 : END DO
770 : END IF
771 :
772 102 : t0 = time_series(t0_i)
773 :
774 : ! Default damping so that at the end of the time series, divide value by e^-4
775 102 : damping = 4.0_dp/(time_series(N) - time_series(t0_i))
776 102 : IF (PRESENT(damping_opt)) THEN
777 102 : IF (damping_opt >= 0.0_dp) damping = damping_opt
778 : END IF
779 :
780 : ! Construct the grid
781 : ! Assume series equidistant
782 102 : delta_t = time_series(2) - time_series(1)
783 102 : delta_omega = twopi/(time_series(N) - time_series(1))
784 7106 : DO i = 1, N
785 7004 : result_series(i) = 0.0_dp
786 7004 : iresult_series(i) = 0.0_dp
787 7004 : omega_series(i) = delta_omega*(i - 1)
788 544714 : DO j = 1, N
789 537608 : j_wrap = MODULO(j + t0_i - 2, N) + 1
790 : result_series(i) = result_series(i) + COS(twopi*(i - 1)*(j - 1)/N)* &
791 537608 : EXP(-damping*delta_t*(j - 1))*value_series(j_wrap)
792 : iresult_series(i) = iresult_series(i) + SIN(twopi*(i - 1)*(j - 1)/N)* &
793 544612 : EXP(-damping*delta_t*(j - 1))*value_series(j_wrap)
794 : END DO
795 : END DO
796 7106 : delta_omega = twopi/(time_series(N) - time_series(1))
797 7106 : result_series(:) = delta_t*result_series(:)
798 7106 : iresult_series(:) = delta_t*iresult_series(:)
799 :
800 102 : CALL timestop(handle)
801 :
802 102 : END SUBROUTINE
803 : END MODULE rt_bse_io
|
