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 Routine for the real time propagation output.
10 : !> \author Florian Schiffmann (02.09)
11 : ! **************************************************************************************************
12 :
13 : MODULE rt_propagation_output
14 : USE atomic_kind_types, ONLY: atomic_kind_type
15 : USE cell_types, ONLY: cell_type
16 : USE cp_control_types, ONLY: dft_control_type,&
17 : rtp_control_type
18 : USE cp_dbcsr_api, ONLY: &
19 : dbcsr_add, dbcsr_binary_write, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
20 : dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_filter, &
21 : dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, dbcsr_iterator_blocks_left, &
22 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
23 : dbcsr_multiply, dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type, &
24 : dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
25 : USE cp_dbcsr_contrib, ONLY: dbcsr_checksum,&
26 : dbcsr_trace
27 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
28 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
29 : dbcsr_allocate_matrix_set,&
30 : dbcsr_deallocate_matrix_set
31 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_double,&
34 : cp_fm_struct_release,&
35 : cp_fm_struct_type
36 : USE cp_fm_types, ONLY: cp_fm_create,&
37 : cp_fm_get_info,&
38 : cp_fm_release,&
39 : cp_fm_set_all,&
40 : cp_fm_type
41 : USE cp_log_handling, ONLY: cp_get_default_logger,&
42 : cp_logger_get_default_io_unit,&
43 : cp_logger_get_default_unit_nr,&
44 : cp_logger_type,&
45 : cp_to_string
46 : USE cp_output_handling, ONLY: cp_iter_string,&
47 : cp_p_file,&
48 : cp_print_key_finished_output,&
49 : cp_print_key_should_output,&
50 : cp_print_key_unit_nr,&
51 : cp_printkey_is_on
52 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
53 : USE efield_utils, ONLY: make_field
54 : USE greenx_interface, ONLY: greenx_refine_ft
55 : USE input_constants, ONLY: ehrenfest,&
56 : real_time_propagation
57 : USE input_section_types, ONLY: section_get_ivals,&
58 : section_vals_get_subs_vals,&
59 : section_vals_type
60 : USE kahan_sum, ONLY: accurate_sum
61 : USE kinds, ONLY: default_path_length,&
62 : dp
63 : USE machine, ONLY: m_flush
64 : USE mathconstants, ONLY: twopi
65 : USE message_passing, ONLY: mp_comm_type
66 : USE parallel_gemm_api, ONLY: parallel_gemm
67 : USE particle_list_types, ONLY: particle_list_type
68 : USE particle_methods, ONLY: get_particle_set
69 : USE particle_types, ONLY: particle_type
70 : USE physcon, ONLY: evolt,&
71 : femtoseconds
72 : USE pw_env_types, ONLY: pw_env_get,&
73 : pw_env_type
74 : USE pw_methods, ONLY: pw_zero
75 : USE pw_pool_types, ONLY: pw_pool_type
76 : USE pw_types, ONLY: pw_c1d_gs_type,&
77 : pw_r3d_rs_type
78 : USE qs_energy_types, ONLY: qs_energy_type
79 : USE qs_environment_types, ONLY: get_qs_env,&
80 : qs_environment_type
81 : USE qs_kind_types, ONLY: get_qs_kind_set,&
82 : qs_kind_type
83 : USE qs_linres_current, ONLY: calculate_jrho_resp
84 : USE qs_linres_types, ONLY: current_env_type
85 : USE qs_mo_io, ONLY: write_rt_mos_to_restart
86 : USE qs_moments, ONLY: build_local_moment_matrix
87 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
88 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
89 : USE qs_rho_types, ONLY: qs_rho_get,&
90 : qs_rho_type
91 : USE qs_scf_post_gpw, ONLY: qs_scf_post_moments,&
92 : write_mo_dependent_results,&
93 : write_mo_free_results
94 : USE qs_scf_post_tb, ONLY: scf_post_calculation_tb
95 : USE qs_scf_types, ONLY: qs_scf_env_type
96 : USE qs_subsys_types, ONLY: qs_subsys_get,&
97 : qs_subsys_type
98 : USE rt_projection_mo_utils, ONLY: compute_and_write_proj_mo
99 : USE rt_propagation_ft, ONLY: multi_fft
100 : USE rt_propagation_types, ONLY: get_rtp,&
101 : rt_prop_type
102 : USE rt_propagation_utils, ONLY: write_rtp_mo_cubes,&
103 : write_rtp_mos_to_output_unit
104 : USE rt_propagation_velocity_gauge, ONLY: velocity_gauge_nl_commutator
105 : #include "../base/base_uses.f90"
106 :
107 : IMPLICIT NONE
108 :
109 : PRIVATE
110 :
111 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_output'
112 :
113 : PUBLIC :: rt_prop_output, &
114 : rt_convergence, &
115 : rt_convergence_density, &
116 : report_density_occupation, &
117 : print_moments, &
118 : calc_local_moment, &
119 : print_ft, &
120 : print_rt_file
121 :
122 : INTEGER, PARAMETER, PUBLIC :: rt_file_comp_both = 0, &
123 : rt_file_comp_real = 1, &
124 : rt_file_comp_imag = 2
125 :
126 : CONTAINS
127 :
128 : ! **************************************************************************************************
129 : !> \brief ...
130 : !> \param qs_env ...
131 : !> \param run_type ...
132 : !> \param delta_iter ...
133 : !> \param used_time ...
134 : ! **************************************************************************************************
135 2408 : SUBROUTINE rt_prop_output(qs_env, run_type, delta_iter, used_time)
136 : TYPE(qs_environment_type), POINTER :: qs_env
137 : INTEGER, INTENT(in) :: run_type
138 : REAL(dp), INTENT(in), OPTIONAL :: delta_iter, used_time
139 :
140 : INTEGER :: i, n_electrons, n_proj, natom, nspin, &
141 : output_unit, spin, unit_nr
142 2408 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
143 : INTEGER, DIMENSION(2) :: nelectron_spin
144 2408 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
145 : LOGICAL :: new_file
146 : REAL(dp) :: orthonormality, strace, tot_rho_r, trace
147 : REAL(kind=dp), DIMENSION(3) :: field, reference_point, vec_pot
148 2408 : REAL(KIND=dp), DIMENSION(:), POINTER :: qs_tot_rho_r
149 2408 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: j_int
150 2408 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
151 : TYPE(cell_type), POINTER :: cell
152 2408 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
153 : TYPE(cp_logger_type), POINTER :: logger
154 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
155 2408 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_cim, matrix_cre, matrix_s, P_im, &
156 2408 : P_re, p_xyz, rho_new
157 : TYPE(dbcsr_type), POINTER :: tmp_ao
158 : TYPE(dft_control_type), POINTER :: dft_control
159 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
160 2408 : POINTER :: sab_all, sab_orb, sap_ppnl
161 2408 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
162 2408 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
163 : TYPE(qs_rho_type), POINTER :: rho
164 : TYPE(rt_prop_type), POINTER :: rtp
165 : TYPE(section_vals_type), POINTER :: dft_section, input, rtp_section
166 :
167 2408 : NULLIFY (logger, dft_control)
168 :
169 4816 : logger => cp_get_default_logger()
170 : CALL get_qs_env(qs_env, &
171 : rtp=rtp, &
172 : matrix_s=matrix_s, &
173 : input=input, &
174 : rho=rho, &
175 : particle_set=particle_set, &
176 : atomic_kind_set=atomic_kind_set, &
177 : qs_kind_set=qs_kind_set, &
178 : dft_control=dft_control, sab_all=sab_all, sab_orb=sab_orb, &
179 : sap_ppnl=sap_ppnl, cell=cell, &
180 2408 : dbcsr_dist=dbcsr_dist, nelectron_spin=nelectron_spin)
181 :
182 2408 : rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
183 :
184 2408 : CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
185 2408 : n_electrons = n_electrons - dft_control%charge
186 :
187 2408 : CALL qs_rho_get(rho_struct=rho, tot_rho_r=qs_tot_rho_r)
188 :
189 2408 : tot_rho_r = accurate_sum(qs_tot_rho_r)
190 :
191 : output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
192 2408 : extension=".scfLog")
193 :
194 2408 : IF (output_unit > 0) THEN
195 : WRITE (output_unit, FMT="(/,(T3,A,T40,I5))") &
196 1204 : "Information at iteration step:", rtp%iter
197 : WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
198 1204 : "Total electronic density (r-space): ", &
199 1204 : tot_rho_r, &
200 : tot_rho_r + &
201 2408 : REAL(n_electrons, dp)
202 : WRITE (UNIT=output_unit, FMT="((T3,A,T59,F22.14))") &
203 1204 : "Total energy:", rtp%energy_new
204 1204 : IF (run_type == ehrenfest) THEN
205 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
206 576 : "Energy difference to previous iteration step:", rtp%energy_new - rtp%energy_old
207 : END IF
208 1204 : IF (run_type == real_time_propagation) THEN
209 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
210 628 : "Energy difference to initial state:", rtp%energy_new - rtp%energy_old
211 : END IF
212 1204 : IF (PRESENT(delta_iter)) THEN
213 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,E20.6))") &
214 1204 : "Convergence:", delta_iter
215 : END IF
216 1204 : IF (rtp%converged) THEN
217 326 : IF (run_type == real_time_propagation) THEN
218 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F12.2))") &
219 188 : "Time needed for propagation:", used_time
220 : END IF
221 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,3X,F16.14))") &
222 326 : "CONVERGENCE REACHED", rtp%energy_new - rtp%energy_old
223 : END IF
224 : END IF
225 :
226 2408 : IF (rtp%converged) THEN
227 652 : IF (.NOT. rtp%linear_scaling) THEN
228 468 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
229 : CALL rt_calculate_orthonormality(orthonormality, &
230 468 : mos_new, matrix_s(1)%matrix)
231 468 : IF (output_unit > 0) THEN
232 : WRITE (output_unit, FMT="(/,(T3,A,T60,F20.10))") &
233 234 : "Max deviation from orthonormalization:", orthonormality
234 : END IF
235 : END IF
236 : END IF
237 :
238 2408 : IF (output_unit > 0) THEN
239 1204 : CALL m_flush(output_unit)
240 : END IF
241 : CALL cp_print_key_finished_output(output_unit, logger, rtp_section, &
242 2408 : "PRINT%PROGRAM_RUN_INFO")
243 :
244 2408 : IF (rtp%converged) THEN
245 652 : dft_section => section_vals_get_subs_vals(input, "DFT")
246 652 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
247 : dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) THEN
248 30 : CALL print_field_applied(qs_env, dft_section)
249 : END IF
250 652 : CALL make_moment(qs_env)
251 652 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
252 : dft_section, "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS"), cp_p_file)) THEN
253 4 : CALL print_rtp_energy_components(qs_env, dft_section)
254 : END IF
255 652 : IF (.NOT. dft_control%qs_control%dftb) THEN
256 532 : CALL write_available_results(qs_env=qs_env, rtp=rtp)
257 : END IF
258 :
259 652 : IF (rtp%linear_scaling) THEN
260 184 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
261 :
262 : ! Probably have to rebuild the moment matrix, since atoms can also move, in principle
263 184 : IF (dft_control%rtp_control%save_local_moments) THEN
264 : ! Save the field value
265 36 : IF (dft_control%apply_efield_field) THEN
266 0 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
267 0 : rtp%fields(:, rtp%istep + rtp%i_start + 1) = CMPLX(field(:), 0.0, kind=dp)
268 : END IF
269 36 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
270 0 : CALL build_local_moment_matrix(qs_env, rtp%local_moments, 1, reference_point)
271 : END IF
272 : ! TODO : Is symmetric rho possible?
273 : ! Spin + complex parts
274 : ! Extensions setup
275 : CALL calc_local_moment(rtp%local_moments, rho_new, &
276 36 : rtp%local_moments_work, rtp%moments(:, :, rtp%istep + rtp%i_start + 1))
277 : ! Time 1 is zero (start) time
278 36 : rtp%times(rtp%istep + rtp%i_start + 1) = qs_env%sim_time
279 36 : output_unit = cp_logger_get_default_io_unit(logger)
280 : CALL print_moments(section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS"), output_unit, &
281 36 : rtp%moments(:, :, rtp%istep + rtp%i_start + 1), qs_env%sim_time, rtp%track_imag_density)
282 : END IF
283 :
284 184 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
285 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART"), cp_p_file)) THEN
286 98 : CALL write_rt_p_to_restart(rho_new, .FALSE.)
287 : END IF
288 184 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
289 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY"), cp_p_file)) THEN
290 2 : CALL write_rt_p_to_restart(rho_new, .TRUE.)
291 : END IF
292 184 : IF (.NOT. dft_control%qs_control%dftb) THEN
293 : !Not sure if these things could also work with dftb or not
294 184 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
295 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
296 64 : DO spin = 1, SIZE(rho_new)/2
297 64 : CALL rt_current(qs_env, rho_new(2*spin)%matrix, dft_section, spin, SIZE(rho_new)/2)
298 : END DO
299 : END IF
300 : END IF
301 : ELSE
302 468 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
303 468 : IF (.NOT. dft_control%qs_control%dftb .AND. .NOT. dft_control%qs_control%xtb) THEN
304 288 : IF (rtp%track_imag_density) THEN
305 210 : NULLIFY (P_im, p_xyz)
306 210 : CALL dbcsr_allocate_matrix_set(p_xyz, 3)
307 :
308 : ! Linear momentum operator
309 : ! prepare for allocation
310 210 : natom = SIZE(particle_set, 1)
311 630 : ALLOCATE (first_sgf(natom))
312 420 : ALLOCATE (last_sgf(natom))
313 : CALL get_particle_set(particle_set, qs_kind_set, &
314 : first_sgf=first_sgf, &
315 210 : last_sgf=last_sgf)
316 420 : ALLOCATE (row_blk_sizes(natom))
317 210 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
318 210 : DEALLOCATE (first_sgf)
319 210 : DEALLOCATE (last_sgf)
320 :
321 210 : ALLOCATE (p_xyz(1)%matrix)
322 : CALL dbcsr_create(matrix=p_xyz(1)%matrix, &
323 : name="p_xyz", &
324 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
325 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
326 210 : mutable_work=.TRUE.)
327 210 : CALL cp_dbcsr_alloc_block_from_nbl(p_xyz(1)%matrix, sab_orb)
328 210 : CALL dbcsr_set(p_xyz(1)%matrix, 0.0_dp)
329 630 : DO i = 2, 3
330 420 : ALLOCATE (p_xyz(i)%matrix)
331 420 : CALL dbcsr_copy(p_xyz(i)%matrix, p_xyz(1)%matrix, "p_xyz-"//TRIM(ADJUSTL(cp_to_string(i))))
332 630 : CALL dbcsr_set(p_xyz(i)%matrix, 0.0_dp)
333 : END DO
334 210 : CALL build_lin_mom_matrix(qs_env, p_xyz)
335 210 : DEALLOCATE (row_blk_sizes)
336 :
337 210 : nspin = SIZE(mos_new)/2
338 210 : CALL qs_rho_get(rho, rho_ao_im=P_im)
339 630 : ALLOCATE (j_int(nspin, 3))
340 1566 : j_int = 0.0_dp
341 :
342 210 : NULLIFY (tmp_ao)
343 210 : CALL dbcsr_init_p(tmp_ao)
344 210 : CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
345 210 : CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
346 210 : CALL dbcsr_set(tmp_ao, 0.0_dp)
347 :
348 840 : DO i = 1, 3
349 630 : strace = 0.0_dp
350 1566 : DO spin = 1, nspin
351 726 : CALL dbcsr_set(tmp_ao, 0.0_dp)
352 : CALL dbcsr_multiply("T", "N", 1.0_dp, P_im(spin)%matrix, p_xyz(i)%matrix, &
353 726 : 0.0_dp, tmp_ao)
354 726 : CALL dbcsr_trace(tmp_ao, trace)
355 726 : strace = strace + trace
356 : ! dft_control%rtp_control%vec_pot(1)
357 1356 : j_int(spin, i) = -trace + dft_control%rtp_control%vec_pot(i)*nelectron_spin(spin)
358 : !! j_int(spin, i) = strace
359 : END DO
360 : END DO
361 :
362 : ! Nonlocal pseudopotential correction to the integrated current:
363 : ! the velocity operator contains i*[Vnl, r], which with C_d = [Vnl~, r_d]
364 : ! contributes j_nl(d) = -Tr[P_re^T C_im(d)] - Tr[P_im^T C_re(d)].
365 : ! Absent in all-electron calculations (no nonlocal projectors).
366 210 : IF (ASSOCIATED(sap_ppnl)) THEN
367 : ! the commutator is taken with the gauge-transformed Vnl only if
368 : ! the propagation Hamiltonian uses it
369 92 : vec_pot = 0.0_dp
370 92 : IF (dft_control%rtp_control%velocity_gauge .AND. &
371 : dft_control%rtp_control%nl_gauge_transform) THEN
372 72 : vec_pot = dft_control%rtp_control%vec_pot
373 : END IF
374 :
375 92 : NULLIFY (matrix_cre, matrix_cim)
376 92 : CALL dbcsr_allocate_matrix_set(matrix_cre, 3)
377 92 : CALL dbcsr_allocate_matrix_set(matrix_cim, 3)
378 368 : DO i = 1, 3
379 276 : CALL dbcsr_init_p(matrix_cre(i)%matrix)
380 : CALL dbcsr_create(matrix_cre(i)%matrix, template=matrix_s(1)%matrix, &
381 : matrix_type=dbcsr_type_antisymmetric, &
382 276 : name="nl commutator real part")
383 276 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_cre(i)%matrix, sab_orb)
384 276 : CALL dbcsr_set(matrix_cre(i)%matrix, 0.0_dp)
385 276 : CALL dbcsr_init_p(matrix_cim(i)%matrix)
386 : CALL dbcsr_create(matrix_cim(i)%matrix, template=matrix_s(1)%matrix, &
387 : matrix_type=dbcsr_type_symmetric, &
388 276 : name="nl commutator imaginary part")
389 276 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_cim(i)%matrix, sab_orb)
390 368 : CALL dbcsr_set(matrix_cim(i)%matrix, 0.0_dp)
391 : END DO
392 :
393 92 : CALL velocity_gauge_nl_commutator(qs_env, matrix_cre, matrix_cim, vec_pot)
394 :
395 92 : CALL qs_rho_get(rho, rho_ao=P_re)
396 368 : DO i = 1, 3
397 692 : DO spin = 1, nspin
398 324 : CALL dbcsr_set(tmp_ao, 0.0_dp)
399 : CALL dbcsr_multiply("T", "N", 1.0_dp, P_re(spin)%matrix, matrix_cim(i)%matrix, &
400 324 : 0.0_dp, tmp_ao)
401 324 : CALL dbcsr_trace(tmp_ao, trace)
402 324 : j_int(spin, i) = j_int(spin, i) - trace
403 324 : CALL dbcsr_set(tmp_ao, 0.0_dp)
404 : CALL dbcsr_multiply("T", "N", 1.0_dp, P_im(spin)%matrix, matrix_cre(i)%matrix, &
405 324 : 0.0_dp, tmp_ao)
406 324 : CALL dbcsr_trace(tmp_ao, trace)
407 924 : j_int(spin, i) = j_int(spin, i) - trace
408 : END DO
409 : END DO
410 :
411 92 : CALL dbcsr_deallocate_matrix_set(matrix_cre)
412 92 : CALL dbcsr_deallocate_matrix_set(matrix_cim)
413 : END IF
414 :
415 210 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
416 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT"), cp_p_file)) THEN
417 :
418 16 : output_unit = cp_logger_get_default_io_unit(logger)
419 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
420 16 : "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT", extension=".dat", is_new_file=new_file)
421 :
422 16 : IF (output_unit > 0) THEN
423 : WRITE (UNIT=output_unit, FMT="(T2,A,E24.16)") &
424 62 : "RTP_CURRENT| CheckSum j_int=", SQRT(SUM(j_int**2))
425 8 : IF (new_file) THEN
426 4 : IF (nspin == 2) THEN
427 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,4X,A,2X,A,2(10X,A),4X,A,2(10X,A))') &
428 1 : "Step Nr.", "Time[fs]", "ALPHA jint[X]", "jint[Y]", "jint[Z]", &
429 2 : "BETA jint[X]", "jint[Y]", "jint[Z]"
430 : ELSE
431 3 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,4X,A,8X,A,2(10X,A))') "Step Nr.", "Time[fs]", &
432 6 : "jint[X]", "jint[Y]", "jint[Z]"
433 : END IF
434 : END IF
435 :
436 8 : IF (nspin == 2) THEN
437 2 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
438 16 : j_int(1, 1:3), j_int(2, 1:3)
439 : ELSE
440 6 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
441 30 : j_int(1, 1:3)
442 : END IF
443 : END IF
444 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
445 16 : "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT")
446 : END IF
447 210 : DEALLOCATE (j_int)
448 :
449 210 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
450 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
451 24 : DO spin = 1, nspin
452 24 : CALL rt_current(qs_env, P_im(spin)%matrix, dft_section, spin, nspin)
453 : END DO
454 : END IF
455 210 : CALL dbcsr_deallocate_matrix(tmp_ao)
456 210 : CALL dbcsr_deallocate_matrix_set(p_xyz)
457 : END IF
458 :
459 : ! projection of molecular orbitals
460 288 : IF (dft_control%rtp_control%is_proj_mo) THEN
461 100 : DO n_proj = 1, SIZE(dft_control%rtp_control%proj_mo_list)
462 : CALL compute_and_write_proj_mo(qs_env, mos_new, &
463 100 : dft_control%rtp_control%proj_mo_list(n_proj)%proj_mo, n_proj)
464 : END DO
465 : END IF
466 : END IF
467 : CALL write_rt_mos_to_restart(qs_env%mos, mos_new, particle_set, &
468 468 : dft_section, qs_kind_set)
469 : END IF
470 : END IF
471 :
472 2408 : rtp%energy_old = rtp%energy_new
473 :
474 2408 : IF (.NOT. rtp%converged .AND. rtp%iter >= dft_control%rtp_control%max_iter) THEN
475 : CALL cp_abort(__LOCATION__, "EMD did not converge, either increase MAX_ITER "// &
476 0 : "or use a smaller TIMESTEP")
477 : END IF
478 :
479 4816 : END SUBROUTINE rt_prop_output
480 :
481 : ! **************************************************************************************************
482 : !> \brief computes the effective orthonormality of a set of mos given an s-matrix
483 : !> orthonormality is the max deviation from unity of the C^T S C
484 : !> \param orthonormality ...
485 : !> \param mos_new ...
486 : !> \param matrix_s ...
487 : !> \author Florian Schiffmann (02.09)
488 : ! **************************************************************************************************
489 468 : SUBROUTINE rt_calculate_orthonormality(orthonormality, mos_new, matrix_s)
490 : REAL(KIND=dp), INTENT(out) :: orthonormality
491 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
492 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
493 :
494 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_calculate_orthonormality'
495 :
496 : INTEGER :: handle, i, im, ispin, j, k, n, &
497 : ncol_local, nrow_local, nspin, re
498 468 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
499 : REAL(KIND=dp) :: alpha, max_alpha, max_beta
500 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
501 : TYPE(cp_fm_type) :: overlap_re, svec_im, svec_re
502 :
503 468 : NULLIFY (tmp_fm_struct)
504 :
505 468 : CALL timeset(routineN, handle)
506 :
507 468 : nspin = SIZE(mos_new)/2
508 468 : max_alpha = 0.0_dp
509 468 : max_beta = 0.0_dp
510 1046 : DO ispin = 1, nspin
511 578 : re = ispin*2 - 1
512 578 : im = ispin*2
513 : ! get S*C
514 578 : CALL cp_fm_create(svec_re, mos_new(im)%matrix_struct)
515 578 : CALL cp_fm_create(svec_im, mos_new(im)%matrix_struct)
516 : CALL cp_fm_get_info(mos_new(im), &
517 578 : nrow_global=n, ncol_global=k)
518 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(re), &
519 578 : svec_re, k)
520 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(im), &
521 578 : svec_im, k)
522 :
523 : ! get C^T (S*C)
524 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
525 : para_env=mos_new(re)%matrix_struct%para_env, &
526 578 : context=mos_new(re)%matrix_struct%context)
527 578 : CALL cp_fm_create(overlap_re, tmp_fm_struct)
528 :
529 578 : CALL cp_fm_struct_release(tmp_fm_struct)
530 :
531 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(re), &
532 578 : svec_re, 0.0_dp, overlap_re)
533 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(im), &
534 578 : svec_im, 1.0_dp, overlap_re)
535 :
536 578 : CALL cp_fm_release(svec_re)
537 578 : CALL cp_fm_release(svec_im)
538 :
539 : CALL cp_fm_get_info(overlap_re, nrow_local=nrow_local, ncol_local=ncol_local, &
540 578 : row_indices=row_indices, col_indices=col_indices)
541 1705 : DO i = 1, nrow_local
542 10720 : DO j = 1, ncol_local
543 9015 : alpha = overlap_re%local_data(i, j)
544 9015 : IF (row_indices(i) == col_indices(j)) alpha = alpha - 1.0_dp
545 10142 : max_alpha = MAX(max_alpha, ABS(alpha))
546 : END DO
547 : END DO
548 2780 : CALL cp_fm_release(overlap_re)
549 : END DO
550 468 : CALL mos_new(1)%matrix_struct%para_env%max(max_alpha)
551 468 : CALL mos_new(1)%matrix_struct%para_env%max(max_beta)
552 468 : orthonormality = max_alpha
553 :
554 468 : CALL timestop(handle)
555 :
556 468 : END SUBROUTINE rt_calculate_orthonormality
557 :
558 : ! **************************************************************************************************
559 : !> \brief computes the convergence criterion for RTP and EMD
560 : !> \param rtp ...
561 : !> \param matrix_s Overlap matrix without the derivatives
562 : !> \param delta_mos ...
563 : !> \param delta_eps ...
564 : !> \author Florian Schiffmann (02.09)
565 : ! **************************************************************************************************
566 :
567 1594 : SUBROUTINE rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
568 : TYPE(rt_prop_type), POINTER :: rtp
569 : TYPE(dbcsr_type), POINTER :: matrix_s
570 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: delta_mos
571 : REAL(dp), INTENT(out) :: delta_eps
572 :
573 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence'
574 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
575 :
576 : INTEGER :: handle, i, icol, im, ispin, j, lcol, &
577 : lrow, nao, newdim, nmo, nspin, re
578 : LOGICAL :: double_col, double_row
579 : REAL(KIND=dp) :: alpha, max_alpha
580 : TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1, tmp_fm_struct
581 : TYPE(cp_fm_type) :: work, work1, work2
582 1594 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
583 :
584 1594 : NULLIFY (tmp_fm_struct)
585 :
586 1594 : CALL timeset(routineN, handle)
587 :
588 1594 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
589 :
590 1594 : nspin = SIZE(delta_mos)/2
591 1594 : max_alpha = 0.0_dp
592 :
593 5654 : DO i = 1, SIZE(mos_new)
594 5654 : CALL cp_fm_scale_and_add(-one, delta_mos(i), one, mos_new(i))
595 : END DO
596 :
597 3624 : DO ispin = 1, nspin
598 2030 : re = ispin*2 - 1
599 2030 : im = ispin*2
600 :
601 2030 : double_col = .TRUE.
602 2030 : double_row = .FALSE.
603 : CALL cp_fm_struct_double(newstruct, &
604 : delta_mos(re)%matrix_struct, &
605 : delta_mos(re)%matrix_struct%context, &
606 : double_col, &
607 2030 : double_row)
608 :
609 2030 : CALL cp_fm_create(work, matrix_struct=newstruct)
610 2030 : CALL cp_fm_create(work1, matrix_struct=newstruct)
611 :
612 : CALL cp_fm_get_info(delta_mos(re), ncol_local=lcol, ncol_global=nmo, &
613 2030 : nrow_global=nao)
614 2030 : CALL cp_fm_get_info(work, ncol_global=newdim)
615 :
616 2030 : CALL cp_fm_set_all(work, zero, zero)
617 :
618 10088 : DO icol = 1, lcol
619 91848 : work%local_data(:, icol) = delta_mos(re)%local_data(:, icol)
620 93878 : work%local_data(:, icol + lcol) = delta_mos(im)%local_data(:, icol)
621 : END DO
622 :
623 2030 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, work, work1, ncol=newdim)
624 :
625 2030 : CALL cp_fm_release(work)
626 :
627 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, &
628 : para_env=delta_mos(re)%matrix_struct%para_env, &
629 2030 : context=delta_mos(re)%matrix_struct%context)
630 : CALL cp_fm_struct_double(newstruct1, &
631 : tmp_fm_struct, &
632 : delta_mos(re)%matrix_struct%context, &
633 : double_col, &
634 2030 : double_row)
635 :
636 2030 : CALL cp_fm_create(work, matrix_struct=newstruct1)
637 2030 : CALL cp_fm_create(work2, matrix_struct=newstruct1)
638 :
639 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(re), &
640 2030 : work1, zero, work)
641 :
642 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(im), &
643 2030 : work1, zero, work2)
644 :
645 2030 : CALL cp_fm_get_info(work, nrow_local=lrow)
646 6059 : DO i = 1, lrow
647 40044 : DO j = 1, lcol
648 : alpha = SQRT((work%local_data(i, j) + work2%local_data(i, j + lcol))**2 + &
649 33985 : (work%local_data(i, j + lcol) - work2%local_data(i, j))**2)
650 38014 : max_alpha = MAX(max_alpha, ABS(alpha))
651 : END DO
652 : END DO
653 :
654 2030 : CALL cp_fm_release(work)
655 2030 : CALL cp_fm_release(work1)
656 2030 : CALL cp_fm_release(work2)
657 2030 : CALL cp_fm_struct_release(tmp_fm_struct)
658 2030 : CALL cp_fm_struct_release(newstruct)
659 9714 : CALL cp_fm_struct_release(newstruct1)
660 :
661 : END DO
662 :
663 1594 : CALL delta_mos(1)%matrix_struct%para_env%max(max_alpha)
664 1594 : delta_eps = SQRT(max_alpha)
665 :
666 1594 : CALL timestop(handle)
667 :
668 1594 : END SUBROUTINE rt_convergence
669 :
670 : ! **************************************************************************************************
671 : !> \brief computes the convergence criterion for RTP and EMD based on the density matrix
672 : !> \param rtp ...
673 : !> \param delta_P ...
674 : !> \param delta_eps ...
675 : !> \author Samuel Andermatt (02.14)
676 : ! **************************************************************************************************
677 :
678 1628 : SUBROUTINE rt_convergence_density(rtp, delta_P, delta_eps)
679 :
680 : TYPE(rt_prop_type), POINTER :: rtp
681 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P
682 : REAL(dp), INTENT(out) :: delta_eps
683 :
684 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence_density'
685 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
686 :
687 : INTEGER :: col_atom, handle, i, ispin, row_atom
688 : REAL(dp) :: alpha, max_alpha
689 814 : REAL(dp), DIMENSION(:, :), POINTER :: block_values
690 : TYPE(dbcsr_iterator_type) :: iter
691 814 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
692 : TYPE(dbcsr_type), POINTER :: tmp
693 : TYPE(mp_comm_type) :: group
694 :
695 814 : CALL timeset(routineN, handle)
696 :
697 814 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
698 :
699 2934 : DO i = 1, SIZE(rho_new)
700 2934 : CALL dbcsr_add(delta_P(i)%matrix, rho_new(i)%matrix, one, -one)
701 : END DO
702 : !get the maximum value of delta_P
703 2934 : DO i = 1, SIZE(delta_P)
704 : !square all entries of both matrices
705 2120 : CALL dbcsr_iterator_start(iter, delta_P(i)%matrix)
706 12550 : DO WHILE (dbcsr_iterator_blocks_left(iter))
707 10430 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
708 808178 : block_values = block_values*block_values
709 : END DO
710 5054 : CALL dbcsr_iterator_stop(iter)
711 : END DO
712 : NULLIFY (tmp)
713 814 : ALLOCATE (tmp)
714 814 : CALL dbcsr_create(tmp, template=delta_P(1)%matrix, matrix_type="N")
715 1874 : DO ispin = 1, SIZE(delta_P)/2
716 1060 : CALL dbcsr_desymmetrize(delta_P(2*ispin - 1)%matrix, tmp)
717 1874 : CALL dbcsr_add(delta_P(2*ispin)%matrix, tmp, one, one)
718 : END DO
719 : !the absolute values are now in the even entries of delta_P
720 814 : max_alpha = zero
721 1874 : DO ispin = 1, SIZE(delta_P)/2
722 1060 : CALL dbcsr_iterator_start(iter, delta_P(2*ispin)%matrix)
723 6372 : DO WHILE (dbcsr_iterator_blocks_left(iter))
724 5312 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
725 408540 : alpha = MAXVAL(block_values)
726 6372 : IF (alpha > max_alpha) max_alpha = alpha
727 : END DO
728 2934 : CALL dbcsr_iterator_stop(iter)
729 : END DO
730 814 : CALL dbcsr_get_info(delta_P(1)%matrix, group=group)
731 814 : CALL group%max(max_alpha)
732 814 : delta_eps = SQRT(max_alpha)
733 814 : CALL dbcsr_deallocate_matrix(tmp)
734 814 : CALL timestop(handle)
735 :
736 814 : END SUBROUTINE rt_convergence_density
737 :
738 : ! **************************************************************************************************
739 : !> \brief interface to qs_moments. Does only work for nonperiodic dipole
740 : !> \param qs_env ...
741 : !> \author Florian Schiffmann (02.09)
742 : ! **************************************************************************************************
743 :
744 652 : SUBROUTINE make_moment(qs_env)
745 :
746 : TYPE(qs_environment_type), POINTER :: qs_env
747 :
748 : CHARACTER(len=*), PARAMETER :: routineN = 'make_moment'
749 :
750 : INTEGER :: handle, output_unit
751 : TYPE(cp_logger_type), POINTER :: logger
752 : TYPE(dft_control_type), POINTER :: dft_control
753 :
754 652 : CALL timeset(routineN, handle)
755 :
756 652 : NULLIFY (dft_control)
757 :
758 652 : logger => cp_get_default_logger()
759 652 : output_unit = cp_logger_get_default_io_unit(logger)
760 652 : CALL get_qs_env(qs_env, dft_control=dft_control)
761 652 : IF (dft_control%qs_control%dftb) THEN
762 120 : CALL scf_post_calculation_tb(qs_env, "DFTB", .FALSE.)
763 532 : ELSE IF (dft_control%qs_control%xtb) THEN
764 60 : CALL scf_post_calculation_tb(qs_env, "xTB", .FALSE.)
765 : ELSE
766 472 : CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
767 : END IF
768 652 : CALL timestop(handle)
769 :
770 652 : END SUBROUTINE make_moment
771 :
772 : ! **************************************************************************************************
773 : !> \brief Reports the sparsity pattern of the complex density matrix
774 : !> \param filter_eps ...
775 : !> \param rho ...
776 : !> \author Samuel Andermatt (09.14)
777 : ! **************************************************************************************************
778 :
779 184 : SUBROUTINE report_density_occupation(filter_eps, rho)
780 :
781 : REAL(KIND=dp) :: filter_eps
782 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho
783 :
784 : CHARACTER(len=*), PARAMETER :: routineN = 'report_density_occupation'
785 :
786 : INTEGER :: handle, i, im, ispin, re, unit_nr
787 : REAL(KIND=dp) :: eps, occ
788 : TYPE(cp_logger_type), POINTER :: logger
789 184 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmp
790 :
791 184 : CALL timeset(routineN, handle)
792 :
793 184 : logger => cp_get_default_logger()
794 184 : unit_nr = cp_logger_get_default_io_unit(logger)
795 184 : NULLIFY (tmp)
796 184 : CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
797 692 : DO i = 1, SIZE(rho)
798 508 : CALL dbcsr_init_p(tmp(i)%matrix)
799 508 : CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
800 692 : CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
801 : END DO
802 438 : DO ispin = 1, SIZE(rho)/2
803 254 : re = 2*ispin - 1
804 254 : im = 2*ispin
805 254 : eps = MAX(filter_eps, 1.0E-11_dp)
806 2354 : DO WHILE (eps < 1.1_dp)
807 2100 : CALL dbcsr_filter(tmp(re)%matrix, eps)
808 2100 : occ = dbcsr_get_occupation(tmp(re)%matrix)
809 3150 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
810 2100 : ispin, " eps ", eps, " real: ", occ
811 2100 : eps = eps*10
812 : END DO
813 254 : eps = MAX(filter_eps, 1.0E-11_dp)
814 2538 : DO WHILE (eps < 1.1_dp)
815 2100 : CALL dbcsr_filter(tmp(im)%matrix, eps)
816 2100 : occ = dbcsr_get_occupation(tmp(im)%matrix)
817 3150 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
818 2100 : ispin, " eps ", eps, " imag: ", occ
819 2100 : eps = eps*10.0_dp
820 : END DO
821 : END DO
822 184 : CALL dbcsr_deallocate_matrix_set(tmp)
823 184 : CALL timestop(handle)
824 :
825 184 : END SUBROUTINE report_density_occupation
826 :
827 : ! **************************************************************************************************
828 : !> \brief Writes the density matrix and the atomic positions to a restart file
829 : !> \param rho_new ...
830 : !> \param history ...
831 : !> \author Samuel Andermatt (09.14)
832 : ! **************************************************************************************************
833 :
834 100 : SUBROUTINE write_rt_p_to_restart(rho_new, history)
835 :
836 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
837 : LOGICAL :: history
838 :
839 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_p_to_restart'
840 :
841 : CHARACTER(LEN=default_path_length) :: file_name, project_name
842 : INTEGER :: handle, im, ispin, re, unit_nr
843 : REAL(KIND=dp) :: cs_pos
844 : TYPE(cp_logger_type), POINTER :: logger
845 :
846 100 : CALL timeset(routineN, handle)
847 100 : logger => cp_get_default_logger()
848 100 : IF (logger%para_env%is_source()) THEN
849 50 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
850 : ELSE
851 : unit_nr = -1
852 : END IF
853 :
854 100 : project_name = logger%iter_info%project_name
855 238 : DO ispin = 1, SIZE(rho_new)/2
856 138 : re = 2*ispin - 1
857 138 : im = 2*ispin
858 138 : IF (history) THEN
859 : WRITE (file_name, '(A,I0,A)') &
860 2 : TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
861 : ELSE
862 136 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
863 : END IF
864 138 : cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.TRUE.)
865 138 : IF (unit_nr > 0) THEN
866 69 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
867 : END IF
868 138 : CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
869 138 : IF (history) THEN
870 : WRITE (file_name, '(A,I0,A)') &
871 2 : TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
872 : ELSE
873 136 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
874 : END IF
875 138 : cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.TRUE.)
876 138 : IF (unit_nr > 0) THEN
877 69 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
878 : END IF
879 238 : CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
880 : END DO
881 :
882 100 : CALL timestop(handle)
883 :
884 100 : END SUBROUTINE write_rt_p_to_restart
885 :
886 : ! **************************************************************************************************
887 : !> \brief Collocation of the current and printing of it in a cube file
888 : !> \param qs_env ...
889 : !> \param P_im ...
890 : !> \param dft_section ...
891 : !> \param spin ...
892 : !> \param nspin ...
893 : !> \author Samuel Andermatt (06.15)
894 : ! **************************************************************************************************
895 48 : SUBROUTINE rt_current(qs_env, P_im, dft_section, spin, nspin)
896 : TYPE(qs_environment_type), POINTER :: qs_env
897 : TYPE(dbcsr_type), POINTER :: P_im
898 : TYPE(section_vals_type), POINTER :: dft_section
899 : INTEGER :: spin, nspin
900 :
901 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_current'
902 :
903 : CHARACTER(len=1) :: char_spin
904 : CHARACTER(len=14) :: ext
905 : CHARACTER(len=2) :: sdir
906 : INTEGER :: dir, handle, print_unit
907 48 : INTEGER, DIMENSION(:), POINTER :: stride
908 : LOGICAL :: mpi_io
909 : TYPE(cp_logger_type), POINTER :: logger
910 : TYPE(current_env_type) :: current_env
911 : TYPE(dbcsr_type), POINTER :: tmp, zero
912 : TYPE(particle_list_type), POINTER :: particles
913 : TYPE(pw_c1d_gs_type) :: gs
914 : TYPE(pw_env_type), POINTER :: pw_env
915 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
916 : TYPE(pw_r3d_rs_type) :: rs
917 : TYPE(qs_subsys_type), POINTER :: subsys
918 :
919 48 : CALL timeset(routineN, handle)
920 :
921 48 : logger => cp_get_default_logger()
922 48 : CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
923 48 : CALL qs_subsys_get(subsys, particles=particles)
924 48 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
925 :
926 48 : NULLIFY (zero, tmp)
927 48 : ALLOCATE (zero, tmp)
928 48 : CALL dbcsr_create(zero, template=P_im)
929 48 : CALL dbcsr_copy(zero, P_im)
930 48 : CALL dbcsr_set(zero, 0.0_dp)
931 48 : CALL dbcsr_create(tmp, template=P_im)
932 48 : CALL dbcsr_copy(tmp, P_im)
933 48 : IF (nspin == 1) THEN
934 32 : CALL dbcsr_scale(tmp, 0.5_dp)
935 : END IF
936 48 : current_env%gauge = -1
937 48 : current_env%gauge_init = .FALSE.
938 48 : CALL auxbas_pw_pool%create_pw(rs)
939 48 : CALL auxbas_pw_pool%create_pw(gs)
940 :
941 : NULLIFY (stride)
942 48 : ALLOCATE (stride(3))
943 :
944 192 : DO dir = 1, 3
945 :
946 144 : CALL pw_zero(rs)
947 144 : CALL pw_zero(gs)
948 :
949 144 : CALL calculate_jrho_resp(zero, tmp, zero, zero, dir, dir, rs, gs, qs_env, current_env, retain_rsgrid=.TRUE.)
950 :
951 576 : stride = section_get_ivals(dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT%STRIDE")
952 :
953 144 : IF (dir == 1) THEN
954 48 : sdir = "-x"
955 96 : ELSE IF (dir == 2) THEN
956 48 : sdir = "-y"
957 : ELSE
958 48 : sdir = "-z"
959 : END IF
960 144 : WRITE (char_spin, "(I1)") spin
961 :
962 144 : ext = "-SPIN-"//char_spin//sdir//".cube"
963 144 : mpi_io = .TRUE.
964 : print_unit = cp_print_key_unit_nr(logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
965 : extension=ext, file_status="REPLACE", file_action="WRITE", &
966 144 : log_filename=.FALSE., mpi_io=mpi_io)
967 :
968 : CALL cp_pw_to_cube(rs, print_unit, "EMD current", particles=particles, stride=stride, &
969 144 : mpi_io=mpi_io)
970 :
971 : CALL cp_print_key_finished_output(print_unit, logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
972 192 : mpi_io=mpi_io)
973 :
974 : END DO
975 :
976 48 : CALL auxbas_pw_pool%give_back_pw(rs)
977 48 : CALL auxbas_pw_pool%give_back_pw(gs)
978 :
979 48 : CALL dbcsr_deallocate_matrix(zero)
980 48 : CALL dbcsr_deallocate_matrix(tmp)
981 :
982 48 : DEALLOCATE (stride)
983 :
984 48 : CALL timestop(handle)
985 :
986 3504 : END SUBROUTINE rt_current
987 :
988 : ! **************************************************************************************************
989 : !> \brief Interface routine to trigger writing of results available from normal
990 : !> SCF. Can write MO-dependent and MO free results (needed for call from
991 : !> the linear scaling code)
992 : !> Update: trigger also some of prints for time-dependent runs
993 : !> \param qs_env ...
994 : !> \param rtp ...
995 : !> \par History
996 : !> 2022-11 Update [Guillaume Le Breton]
997 : ! **************************************************************************************************
998 532 : SUBROUTINE write_available_results(qs_env, rtp)
999 : TYPE(qs_environment_type), POINTER :: qs_env
1000 : TYPE(rt_prop_type), POINTER :: rtp
1001 :
1002 : CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results'
1003 :
1004 : INTEGER :: handle
1005 : TYPE(qs_scf_env_type), POINTER :: scf_env
1006 :
1007 532 : CALL timeset(routineN, handle)
1008 :
1009 532 : CALL get_qs_env(qs_env, scf_env=scf_env)
1010 532 : IF (rtp%linear_scaling) THEN
1011 184 : CALL write_mo_free_results(qs_env)
1012 : ELSE
1013 348 : CALL write_mo_free_results(qs_env)
1014 348 : CALL write_mo_dependent_results(qs_env, scf_env)
1015 : ! Time-dependent MO print
1016 348 : CALL write_rtp_mos_to_output_unit(qs_env, rtp)
1017 348 : CALL write_rtp_mo_cubes(qs_env, rtp)
1018 : END IF
1019 :
1020 532 : CALL timestop(handle)
1021 :
1022 532 : END SUBROUTINE write_available_results
1023 :
1024 : ! **************************************************************************************************
1025 : !> \brief Print the field applied to the system. Either the electric
1026 : !> field or the vector potential depending on the gauge used
1027 : !> \param qs_env ...
1028 : !> \param dft_section ...
1029 : !> \par History
1030 : !> 2023-01 Created [Guillaume Le Breton]
1031 : ! **************************************************************************************************
1032 30 : SUBROUTINE print_field_applied(qs_env, dft_section)
1033 : TYPE(qs_environment_type), POINTER :: qs_env
1034 : TYPE(section_vals_type), POINTER :: dft_section
1035 :
1036 : CHARACTER(LEN=3), DIMENSION(3) :: rlab
1037 : CHARACTER(LEN=default_path_length) :: filename
1038 : INTEGER :: i, i_step, output_unit, unit_nr
1039 : LOGICAL :: new_file
1040 : REAL(kind=dp) :: field(3)
1041 : TYPE(cp_logger_type), POINTER :: logger
1042 : TYPE(dft_control_type), POINTER :: dft_control
1043 : TYPE(rt_prop_type), POINTER :: rtp
1044 :
1045 30 : NULLIFY (dft_control)
1046 :
1047 30 : logger => cp_get_default_logger()
1048 30 : output_unit = cp_logger_get_default_io_unit(logger)
1049 :
1050 30 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp)
1051 :
1052 30 : i_step = rtp%istep
1053 :
1054 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
1055 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat", is_new_file=new_file)
1056 :
1057 30 : IF (output_unit > 0) THEN
1058 60 : rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
1059 15 : IF (unit_nr /= output_unit) THEN
1060 15 : INQUIRE (UNIT=unit_nr, NAME=filename)
1061 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
1062 15 : "FIELD", "The field applied is written to the file:", &
1063 30 : TRIM(filename)
1064 : ELSE
1065 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED [a.u.]"
1066 : WRITE (UNIT=output_unit, FMT="(T5,3(A,A,E16.8,1X))") &
1067 0 : (TRIM(rlab(i)), "=", dft_control%rtp_control%field(i), i=1, 3)
1068 : END IF
1069 :
1070 15 : IF (new_file) THEN
1071 2 : IF (dft_control%apply_efield_field) THEN
1072 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,3(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z"
1073 1 : ELSE IF (dft_control%apply_vector_potential) THEN
1074 0 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,6(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z", &
1075 0 : " Vec. Pot. X", " Vec. Pot. Y", " Vec. Pot. Z"
1076 : END IF
1077 : END IF
1078 :
1079 15 : field = 0.0_dp
1080 15 : IF (dft_control%apply_efield_field) THEN
1081 4 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
1082 4 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1083 8 : field(1), field(2), field(3)
1084 : ! DO i=1,3
1085 : ! IF (ABS(field(i))< 10E-10) field(i) = 0.0_dp
1086 : ! END IF
1087 11 : ELSE IF (dft_control%apply_vector_potential) THEN
1088 9 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1089 9 : dft_control%rtp_control%field(1), dft_control%rtp_control%field(2), dft_control%rtp_control%field(3), &
1090 18 : dft_control%rtp_control%vec_pot(1), dft_control%rtp_control%vec_pot(2), dft_control%rtp_control%vec_pot(3)
1091 : END IF
1092 :
1093 : END IF
1094 :
1095 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
1096 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD")
1097 :
1098 30 : END SUBROUTINE print_field_applied
1099 :
1100 : ! **************************************************************************************************
1101 : !> \brief Print the components of the total energy used in an RTP calculation
1102 : !> \param qs_env ...
1103 : !> \param dft_section ...
1104 : !> \par History
1105 : !> 2024-02 Created [ANB]
1106 : ! **************************************************************************************************
1107 4 : SUBROUTINE print_rtp_energy_components(qs_env, dft_section)
1108 : TYPE(qs_environment_type), POINTER :: qs_env
1109 : TYPE(section_vals_type), POINTER :: dft_section
1110 :
1111 : CHARACTER(LEN=default_path_length) :: filename
1112 : INTEGER :: i_step, output_unit, unit_nr
1113 : LOGICAL :: new_file
1114 : TYPE(cp_logger_type), POINTER :: logger
1115 : TYPE(dft_control_type), POINTER :: dft_control
1116 : TYPE(qs_energy_type), POINTER :: energy
1117 : TYPE(rt_prop_type), POINTER :: rtp
1118 :
1119 4 : NULLIFY (dft_control, energy, rtp)
1120 :
1121 4 : logger => cp_get_default_logger()
1122 4 : output_unit = cp_logger_get_default_io_unit(logger)
1123 :
1124 4 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
1125 4 : i_step = rtp%istep
1126 :
1127 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
1128 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS", extension=".ener", &
1129 4 : file_action="WRITE", is_new_file=new_file)
1130 :
1131 4 : IF (output_unit > 0) THEN
1132 2 : IF (unit_nr /= output_unit) THEN
1133 2 : INQUIRE (UNIT=unit_nr, NAME=filename)
1134 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
1135 2 : "ENERGY_CONSTITUENTS", "Total Energy constituents written to file:", &
1136 4 : TRIM(filename)
1137 : ELSE
1138 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ENERGY_CONSTITUENTS"
1139 : END IF
1140 :
1141 2 : IF (new_file) THEN
1142 : ! NOTE that these are not all terms contributing to the total energy for RTP, only a selection of those
1143 : ! most significant / impactful. Therefore the printed components likely will not add up to the total energy.
1144 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,10(6X,A))') "Step Nr.", "Time[fs]", &
1145 1 : "Total ener.[a.u.]", "core[a.u.] ", " overlap [a.u.]", "hartree[a.u.]", " exc. [a.u.] ", &
1146 2 : " hartree 1c[a.u.]", "exc 1c[a.u.] ", "exc admm[a.u.]", "exc 1c admm[a.u.]", "efield LG"
1147 :
1148 : END IF
1149 : WRITE (UNIT=unit_nr, FMT="(I10,F20.6,10(F20.9))") &
1150 2 : qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1151 2 : energy%total, energy%core, energy%core_overlap, energy%hartree, energy%exc, &
1152 4 : energy%hartree_1c, energy%exc1, energy%exc_aux_fit, energy%exc1_aux_fit, energy%efield_core
1153 :
1154 : END IF
1155 :
1156 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
1157 4 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS")
1158 :
1159 4 : END SUBROUTINE print_rtp_energy_components
1160 :
1161 : ! **************************************************************************************************
1162 : !> \brief Print the dipole moments into a file
1163 : !> \param moments_section Section of the input defining the file/stream to print the moments to
1164 : !> \param info_unit Unit where standard output from the program is written - for add. identifiers
1165 : !> \param moments Actual moment values (for specific time step)
1166 : !> \param time Current simulation time
1167 : !> \param imag_opt Whether to calculate the imaginary part
1168 : !> \param append_opt ...
1169 : !> \par History
1170 : !> 10.2025 Created [Marek]
1171 : ! **************************************************************************************************
1172 1616 : SUBROUTINE print_moments(moments_section, info_unit, moments, time, imag_opt, append_opt)
1173 : TYPE(section_vals_type), POINTER :: moments_section
1174 : INTEGER :: info_unit
1175 : COMPLEX(kind=dp), DIMENSION(:, :) :: moments
1176 : REAL(kind=dp), OPTIONAL :: time
1177 : LOGICAL, OPTIONAL :: imag_opt, append_opt
1178 :
1179 : CHARACTER(len=14), DIMENSION(4) :: file_extensions
1180 : CHARACTER(len=21) :: prefix
1181 : COMPLEX(kind=dp), DIMENSION(3, 1) :: moment_t
1182 : INTEGER :: i, ndir, nspin, print_unit
1183 : LOGICAL :: append, imaginary
1184 : TYPE(cp_logger_type), POINTER :: logger
1185 :
1186 : ! Index 1 : spin, Index 2 : direction
1187 :
1188 1616 : nspin = SIZE(moments, 1)
1189 1616 : ndir = SIZE(moments, 2)
1190 :
1191 1616 : IF (nspin < 1) CPABORT("Zero spin index size in print moments!")
1192 1616 : IF (ndir < 1) CPABORT("Zero direction index size in print moments!")
1193 :
1194 1616 : imaginary = .TRUE.
1195 1616 : IF (PRESENT(imag_opt)) imaginary = imag_opt
1196 :
1197 1616 : append = .TRUE.
1198 1616 : IF (PRESENT(append_opt)) append = append_opt
1199 :
1200 : ! Get the program run info unit and target unit
1201 : ! If these are the same (most likely the case of __STD_OUT__), add
1202 : ! extra identifier to the printed output
1203 1616 : file_extensions(1) = "_SPIN_A_RE.dat"
1204 1616 : file_extensions(2) = "_SPIN_A_IM.dat"
1205 1616 : file_extensions(3) = "_SPIN_B_RE.dat"
1206 1616 : file_extensions(4) = "_SPIN_B_IM.dat"
1207 1616 : logger => cp_get_default_logger()
1208 3394 : DO i = 1, nspin
1209 7112 : moment_t(:, 1) = moments(i, :)
1210 : ! Real part
1211 1778 : print_unit = cp_print_key_unit_nr(logger, moments_section, extension=file_extensions(2*i - 1))
1212 1778 : IF (print_unit == info_unit) THEN
1213 : ! print with prefix
1214 1151 : prefix = " MOMENTS_TRACE_RE|"
1215 1151 : IF (append) THEN
1216 : ! Print without headers
1217 : CALL print_rt_file(print_unit, xvals=[time], yvals=moment_t, &
1218 : prefix=prefix, prefix_format="(A18)", &
1219 2236 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
1220 : ELSE
1221 : ! Print with headers
1222 : CALL print_rt_file(print_unit, &
1223 : headers=["# Time [fs]", " re(mom_t) x [at.u.]", &
1224 : " re(mom_t) y [at.u.]", " re(mom_t) z [at.u.]"], &
1225 : xvals=[time], yvals=moment_t, &
1226 : prefix=prefix, prefix_format="(A18)", &
1227 198 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
1228 : END IF
1229 : ELSE
1230 : ! Print without prefix
1231 627 : IF (append) THEN
1232 : ! Print without headers
1233 : CALL print_rt_file(print_unit, xvals=[time], yvals=moment_t, &
1234 1200 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
1235 : ELSE
1236 : ! Print with headers
1237 : CALL print_rt_file(print_unit, &
1238 : headers=["# Time [fs]", " re(mom_t) x [at.u.]", &
1239 : " re(mom_t) y [at.u.]", " re(mom_t) z [at.u.]"], &
1240 : xvals=[time], yvals=moment_t, &
1241 162 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
1242 : END IF
1243 : END IF
1244 1778 : CALL cp_print_key_finished_output(print_unit, logger, moments_section)
1245 : ! Same for imaginary part
1246 3394 : IF (imaginary) THEN
1247 1688 : print_unit = cp_print_key_unit_nr(logger, moments_section, extension=file_extensions(2*i))
1248 1688 : IF (print_unit == info_unit) THEN
1249 : ! print with prefix
1250 1106 : prefix = " MOMENTS_TRACE_IM|"
1251 1106 : IF (append) THEN
1252 : ! Print without headers
1253 : CALL print_rt_file(print_unit, xvals=[time], yvals=moment_t, &
1254 : prefix=prefix, prefix_format="(A18)", &
1255 2146 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_imag)
1256 : ELSE
1257 : ! Print with headers
1258 : CALL print_rt_file(print_unit, &
1259 : headers=["# Time [fs]", " im(mom_t) x [at.u.]", &
1260 : " im(mom_t) y [at.u.]", " im(mom_t) z [at.u.]"], &
1261 : xvals=[time], yvals=moment_t, &
1262 : prefix=prefix, prefix_format="(A18)", &
1263 198 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_imag)
1264 : END IF
1265 : ELSE
1266 : ! Print without prefix
1267 582 : IF (append) THEN
1268 : ! Print without headers
1269 : CALL print_rt_file(print_unit, xvals=[time], yvals=moment_t, &
1270 1110 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_imag)
1271 : ELSE
1272 : ! Print with headers
1273 : CALL print_rt_file(print_unit, &
1274 : headers=["# Time [fs]", " im(mom_t) x [at.u.]", &
1275 : " im(mom_t) y [at.u.]", " im(mom_t) z [at.u.]"], &
1276 : xvals=[time], yvals=moment_t, &
1277 162 : xscale_opt=femtoseconds, comp_opt=rt_file_comp_imag)
1278 : END IF
1279 : END IF
1280 1688 : CALL cp_print_key_finished_output(print_unit, logger, moments_section)
1281 : END IF
1282 : END DO
1283 :
1284 1616 : END SUBROUTINE print_moments
1285 :
1286 : ! **************************************************************************************************
1287 : !> \brief Calculate the values of real/imaginary parts of moments in all directions
1288 : !> \param moment_matrices Local matrix representations of dipole (position) operator
1289 : !> \param density_matrices Density matrices (spin and real+complex parts)
1290 : !> \param work Extra dbcsr matrix for work
1291 : !> \param moment Resulting moments (spin and direction)
1292 : !> \param imag_opt Whether to calculate the imaginary part of the moment
1293 : !> \par History
1294 : !> 10.2025 Created [Marek]
1295 : ! **************************************************************************************************
1296 54 : SUBROUTINE calc_local_moment(moment_matrices, density_matrices, work, moment, imag_opt)
1297 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moment_matrices, density_matrices
1298 : TYPE(dbcsr_type) :: work
1299 : COMPLEX(kind=dp), DIMENSION(:, :) :: moment
1300 : LOGICAL, OPTIONAL :: imag_opt
1301 :
1302 : INTEGER :: i, k, nspin
1303 : LOGICAL :: imag
1304 : REAL(kind=dp) :: real_moment
1305 :
1306 54 : imag = .FALSE.
1307 54 : IF (PRESENT(imag_opt)) imag = imag_opt
1308 54 : nspin = SIZE(density_matrices)/2
1309 :
1310 144 : DO i = 1, nspin
1311 414 : DO k = 1, 3
1312 : CALL dbcsr_multiply("N", "N", -1.0_dp, &
1313 : density_matrices(2*i - 1)%matrix, moment_matrices(k)%matrix, &
1314 270 : 0.0_dp, work)
1315 270 : CALL dbcsr_trace(work, real_moment)
1316 270 : moment(i, k) = CMPLX(real_moment, 0.0, kind=dp)
1317 360 : IF (imag) THEN
1318 : CALL dbcsr_multiply("N", "N", -1.0_dp, &
1319 : density_matrices(2*i)%matrix, moment_matrices(k)%matrix, &
1320 0 : 0.0_dp, work)
1321 0 : CALL dbcsr_trace(work, real_moment)
1322 0 : moment(i, k) = moment(i, k) + CMPLX(0.0, real_moment, kind=dp)
1323 : END IF
1324 : END DO
1325 : END DO
1326 :
1327 54 : END SUBROUTINE calc_local_moment
1328 :
1329 : ! **************************************************************************************************
1330 : !> \brief Calculate and print the Fourier transforms + polarizabilites from moment trace
1331 : !> \param rtp_section The RTP input section (needed to access PRINT configurations)
1332 : !> \param moments Moment trace
1333 : !> \param times Corresponding times
1334 : !> \param fields Corresponding fields
1335 : !> \param rtc rt_control_type that includes metadata
1336 : !> \param info_opt ...
1337 : !> \param cell If present, used to change the delta peak representation to be in units of reciprocal lattice
1338 : !> \par History
1339 : !> 10.2025 Created [Marek]
1340 : ! **************************************************************************************************
1341 90 : SUBROUTINE print_ft(rtp_section, moments, times, fields, rtc, info_opt, cell)
1342 : TYPE(section_vals_type), POINTER :: rtp_section
1343 : COMPLEX(kind=dp), DIMENSION(:, :, :), POINTER :: moments
1344 : REAL(kind=dp), DIMENSION(:), POINTER :: times
1345 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: fields
1346 : TYPE(rtp_control_type), POINTER :: rtc
1347 : INTEGER, OPTIONAL :: info_opt
1348 : TYPE(cell_type), OPTIONAL, POINTER :: cell
1349 :
1350 : CHARACTER(len=*), PARAMETER :: routineN = 'print_ft'
1351 :
1352 : CHARACTER(len=11), DIMENSION(2) :: file_extensions
1353 90 : CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: headers
1354 : CHARACTER(len=21) :: prefix
1355 : CHARACTER(len=5) :: prefix_format
1356 90 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: omegas_complex, omegas_pade
1357 90 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: field_results, field_results_pade, &
1358 90 : pol_results, pol_results_pade, pol_results_pade_spin_total, pol_results_spin_total, &
1359 90 : results, results_pade, results_pade_spin_total, results_spin_total, value_series
1360 : INTEGER :: ft_unit, handle, i, idx_omega_zero, &
1361 : info_unit, k, k_static, n, n_elems, &
1362 : n_pade, nspin
1363 : LOGICAL :: do_moments_ft, do_polarizability
1364 : REAL(kind=dp) :: damping, t0
1365 90 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: omegas, omegas_pade_real
1366 : REAL(kind=dp), DIMENSION(3) :: delta_vec
1367 : TYPE(cp_logger_type), POINTER :: logger
1368 : TYPE(section_vals_type), POINTER :: moment_ft_section, pol_section
1369 :
1370 90 : CALL timeset(routineN, handle)
1371 : ! For results, using spin * direction for first index, e.g. for nspin = 2
1372 : ! results(1,:) = (spin=1 and direction=1,:),
1373 : ! results(5,:) = (spin=2 and direction=2,:)
1374 :
1375 90 : logger => cp_get_default_logger()
1376 :
1377 90 : moment_ft_section => section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS_FT")
1378 90 : pol_section => section_vals_get_subs_vals(rtp_section, "PRINT%POLARIZABILITY")
1379 :
1380 90 : nspin = SIZE(moments, 1)
1381 90 : n = SIZE(times)
1382 90 : n_elems = SIZE(rtc%print_pol_elements, 1)
1383 :
1384 90 : info_unit = -1
1385 90 : IF (PRESENT(info_opt)) info_unit = info_opt
1386 :
1387 : ! NOTE : Allows for at most 2 spin species
1388 90 : file_extensions(1) = "_SPIN_A.dat"
1389 90 : file_extensions(2) = "_SPIN_B.dat"
1390 :
1391 : ! Determine whether MOMENTS_FT and/or polarizability needs to be calculated
1392 90 : do_moments_ft = cp_printkey_is_on(logger%iter_info, moment_ft_section)
1393 90 : do_polarizability = cp_printkey_is_on(logger%iter_info, pol_section)
1394 90 : do_polarizability = do_polarizability .AND. (n_elems > 0)
1395 :
1396 90 : damping = rtc%ft_damping
1397 90 : t0 = rtc%ft_t0
1398 :
1399 : ! Determine field ft if polarizability required
1400 90 : IF (do_polarizability) THEN
1401 204 : ALLOCATE (field_results(3, n))
1402 68 : IF (rtc%apply_delta_pulse) THEN
1403 : ! Constant real FT
1404 64 : IF (PRESENT(cell)) THEN
1405 : delta_vec(:) = (REAL(rtc%delta_pulse_direction(1), kind=dp)*cell%h_inv(1, :) + &
1406 : REAL(rtc%delta_pulse_direction(2), kind=dp)*cell%h_inv(2, :) + &
1407 : REAL(rtc%delta_pulse_direction(3), kind=dp)*cell%h_inv(3, :)) &
1408 0 : *twopi*rtc%delta_pulse_scale
1409 : ELSE
1410 256 : delta_vec(:) = REAL(rtc%delta_pulse_direction(:), kind=dp)*rtc%delta_pulse_scale
1411 : END IF
1412 256 : DO k = 1, 3
1413 4180 : field_results(k, :) = CMPLX(delta_vec(k), 0.0, kind=dp)
1414 : END DO
1415 : ELSE
1416 : ! Do explicit FT
1417 : CALL multi_fft(times, fields, field_results, &
1418 4 : damping_opt=damping, t0_opt=t0, subtract_initial_opt=.TRUE.)
1419 : END IF
1420 : END IF
1421 :
1422 90 : IF (do_moments_ft .OR. do_polarizability) THEN
1423 : ! We need to transform at least the moments
1424 : ! NOTE : Might be able to save some memory by only doing FT of actually
1425 : ! required moments, but for now, doing FT of all moment directions
1426 272 : ALLOCATE (results(3*nspin, n))
1427 204 : ALLOCATE (omegas(n))
1428 204 : ALLOCATE (value_series(3*nspin, n))
1429 144 : DO i = 1, nspin
1430 372 : DO k = 1, 3
1431 5284 : value_series(3*(i - 1) + k, :) = moments(i, k, :)
1432 : END DO
1433 : END DO
1434 : ! TODO : Choose whether the initial subtraction is applied in &FT section?
1435 : CALL multi_fft(times, value_series, results, omegas, &
1436 68 : damping_opt=damping, t0_opt=t0, subtract_initial_opt=.TRUE.)
1437 68 : DEALLOCATE (value_series)
1438 144 : DO i = 1, nspin
1439 : ! Output to FT file, if needed
1440 : ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension=file_extensions(i), &
1441 76 : file_form="FORMATTED", file_position="REWIND")
1442 : ! Print header
1443 76 : IF (ft_unit > 0) THEN
1444 38 : ALLOCATE (headers(7))
1445 38 : headers(2) = " x,real [at.u.]"
1446 38 : headers(3) = " x,imag [at.u.]"
1447 38 : headers(4) = " y,real [at.u.]"
1448 38 : headers(5) = " y,imag [at.u.]"
1449 38 : headers(6) = " z,real [at.u.]"
1450 38 : headers(7) = " z,imag [at.u.]"
1451 38 : IF (info_unit == ft_unit) THEN
1452 0 : headers(1) = "# Energy [eV]"
1453 0 : prefix = " MOMENTS_FT|"
1454 0 : prefix_format = "(A12)"
1455 : CALL print_rt_file(ft_unit, headers, omegas, results(3*(i - 1) + 1:3*(i - 1) + 3, :), &
1456 0 : prefix, prefix_format, evolt)
1457 : ELSE
1458 38 : headers(1) = "# omega [at.u.]"
1459 38 : CALL print_rt_file(ft_unit, headers, omegas, results(3*(i - 1) + 1:3*(i - 1) + 3, :))
1460 : END IF
1461 38 : DEALLOCATE (headers)
1462 : END IF
1463 144 : CALL cp_print_key_finished_output(ft_unit, logger, moment_ft_section)
1464 : END DO
1465 : ! Spin-summed total moments FT (open shell only; inert for nspin=1)
1466 68 : IF (nspin > 1) THEN
1467 24 : ALLOCATE (results_spin_total(3, n))
1468 8 : results_spin_total(:, :) = (0.0_dp, 0.0_dp)
1469 24 : DO i = 1, nspin
1470 72 : DO k = 1, 3
1471 952 : results_spin_total(k, :) = results_spin_total(k, :) + results(3*(i - 1) + k, :)
1472 : END DO
1473 : END DO
1474 : ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension="_SPIN_TOTAL.dat", &
1475 8 : file_form="FORMATTED", file_position="REWIND")
1476 8 : IF (ft_unit > 0) THEN
1477 4 : ALLOCATE (headers(7))
1478 4 : headers(2) = " x,real [at.u.]"
1479 4 : headers(3) = " x,imag [at.u.]"
1480 4 : headers(4) = " y,real [at.u.]"
1481 4 : headers(5) = " y,imag [at.u.]"
1482 4 : headers(6) = " z,real [at.u.]"
1483 4 : headers(7) = " z,imag [at.u.]"
1484 4 : IF (info_unit == ft_unit) THEN
1485 0 : headers(1) = "# Energy [eV]"
1486 0 : prefix = " MOMENTS_FT|"
1487 0 : prefix_format = "(A12)"
1488 : CALL print_rt_file(ft_unit, headers, omegas, results_spin_total, &
1489 0 : prefix, prefix_format, evolt)
1490 : ELSE
1491 4 : headers(1) = "# omega [at.u.]"
1492 4 : CALL print_rt_file(ft_unit, headers, omegas, results_spin_total)
1493 : END IF
1494 4 : DEALLOCATE (headers)
1495 : END IF
1496 8 : CALL cp_print_key_finished_output(ft_unit, logger, moment_ft_section)
1497 8 : DEALLOCATE (results_spin_total)
1498 : END IF
1499 : END IF
1500 :
1501 90 : IF (rtc%pade_requested .AND. (do_moments_ft .OR. do_polarizability)) THEN
1502 6 : ALLOCATE (omegas_complex(SIZE(omegas)))
1503 64 : omegas_complex(:) = CMPLX(omegas(:), 0.0, kind=dp)
1504 2 : n_pade = INT((rtc%pade_e_max - rtc%pade_e_min)/rtc%pade_e_step)
1505 6 : ALLOCATE (omegas_pade(n_pade))
1506 6 : ALLOCATE (omegas_pade_real(n_pade))
1507 : ! Construct omegas_pade and omegas_complex
1508 2000 : DO i = 1, n_pade
1509 1998 : omegas_pade_real(i) = (i - 1)*rtc%pade_e_step + rtc%pade_e_min
1510 2000 : omegas_pade(i) = CMPLX(omegas_pade_real(i), 0.0, kind=dp)
1511 : END DO
1512 8 : ALLOCATE (results_pade(nspin*3, n_pade), source=CMPLX(0.0, 0.0, kind=dp))
1513 4 : DO i = 1, nspin
1514 8 : DO k = 1, 3
1515 : CALL greenx_refine_ft(rtc%pade_fit_e_min, rtc%pade_fit_e_max, omegas_complex, results(3*(i - 1) + k, :), &
1516 8 : omegas_pade, results_pade(3*(i - 1) + k, :))
1517 : END DO
1518 : ! Print to a file
1519 : ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension="_PADE"//file_extensions(i), &
1520 2 : file_form="FORMATTED", file_position="REWIND")
1521 4 : IF (ft_unit > 0) THEN
1522 1 : ALLOCATE (headers(7))
1523 1 : headers(2) = " x,real,pade [at.u.]"
1524 1 : headers(3) = " x,imag,pade [at.u.]"
1525 1 : headers(4) = " y,real,pade [at.u.]"
1526 1 : headers(5) = " y,imag,pade [at.u.]"
1527 1 : headers(6) = " z,real,pade [at.u.]"
1528 1 : headers(7) = " z,imag,pade [at.u.]"
1529 1 : IF (info_unit == ft_unit) THEN
1530 0 : headers(1) = "# Energy [eV]"
1531 0 : prefix = " MOMENTS_FT_PADE|"
1532 0 : prefix_format = "(A17)"
1533 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade(3*(i - 1) + 1:3*(i - 1) + 3, :), &
1534 0 : prefix, prefix_format, evolt)
1535 : ELSE
1536 1 : headers(1) = "# omega [at.u.]"
1537 1 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade(3*(i - 1) + 1:3*(i - 1) + 3, :))
1538 : END IF
1539 1 : DEALLOCATE (headers)
1540 : END IF
1541 : END DO
1542 : ! Spin-summed total moments-FT Padé (open shell only; inert for nspin=1)
1543 2 : IF (nspin > 1) THEN
1544 0 : ALLOCATE (results_pade_spin_total(3, n_pade))
1545 0 : results_pade_spin_total(:, :) = (0.0_dp, 0.0_dp)
1546 0 : DO i = 1, nspin
1547 0 : DO k = 1, 3
1548 0 : results_pade_spin_total(k, :) = results_pade_spin_total(k, :) + results_pade(3*(i - 1) + k, :)
1549 : END DO
1550 : END DO
1551 : ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension="_PADE_SPIN_TOTAL.dat", &
1552 0 : file_form="FORMATTED", file_position="REWIND")
1553 0 : IF (ft_unit > 0) THEN
1554 0 : ALLOCATE (headers(7))
1555 0 : headers(2) = " x,real,pade [at.u.]"
1556 0 : headers(3) = " x,imag,pade [at.u.]"
1557 0 : headers(4) = " y,real,pade [at.u.]"
1558 0 : headers(5) = " y,imag,pade [at.u.]"
1559 0 : headers(6) = " z,real,pade [at.u.]"
1560 0 : headers(7) = " z,imag,pade [at.u.]"
1561 0 : IF (info_unit == ft_unit) THEN
1562 0 : headers(1) = "# Energy [eV]"
1563 0 : prefix = " MOMENTS_FT_PADE|"
1564 0 : prefix_format = "(A17)"
1565 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade_spin_total, &
1566 0 : prefix, prefix_format, evolt)
1567 : ELSE
1568 0 : headers(1) = "# omega [at.u.]"
1569 0 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade_spin_total)
1570 : END IF
1571 0 : DEALLOCATE (headers)
1572 : END IF
1573 0 : DEALLOCATE (results_pade_spin_total)
1574 : END IF
1575 : END IF
1576 :
1577 90 : IF (do_polarizability) THEN
1578 : ! get the polarizability elements, as required
1579 272 : ALLOCATE (pol_results(n_elems, n))
1580 144 : DO i = 1, nspin
1581 304 : DO k = 1, n_elems
1582 : ! NOTE - field is regularized to small value
1583 : pol_results(k, :) = results(3*(i - 1) + &
1584 : rtc%print_pol_elements(k, 1), :)/ &
1585 : (field_results(rtc%print_pol_elements(k, 2), :) + &
1586 5284 : 1.0e-10*field_results(rtc%print_pol_elements(k, 2), 2))
1587 : END DO
1588 : ! Print to the file
1589 : ft_unit = cp_print_key_unit_nr(logger, pol_section, extension=file_extensions(i), &
1590 76 : file_form="FORMATTED", file_position="REWIND")
1591 76 : IF (ft_unit > 0) THEN
1592 114 : ALLOCATE (headers(2*n_elems + 1))
1593 152 : DO k = 1, n_elems
1594 114 : WRITE (headers(2*k), "(A16,I2,I2)") "real pol. elem.", &
1595 114 : rtc%print_pol_elements(k, 1), &
1596 228 : rtc%print_pol_elements(k, 2)
1597 114 : WRITE (headers(2*k + 1), "(A16,I2,I2)") "imag pol. elem.", &
1598 114 : rtc%print_pol_elements(k, 1), &
1599 266 : rtc%print_pol_elements(k, 2)
1600 : END DO
1601 : ! Write header
1602 38 : IF (info_unit == ft_unit) THEN
1603 2 : headers(1) = "# Energy [eV]"
1604 2 : prefix = " POLARIZABILITY|"
1605 2 : prefix_format = "(A16)"
1606 : CALL print_rt_file(ft_unit, headers, omegas, pol_results, &
1607 2 : prefix, prefix_format, evolt)
1608 : ELSE
1609 36 : headers(1) = "# omega [at.u.]"
1610 36 : CALL print_rt_file(ft_unit, headers, omegas, pol_results)
1611 : END IF
1612 38 : DEALLOCATE (headers)
1613 : END IF
1614 76 : CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1615 : ! Static polarizability alpha(0): pol_results at the FFT-grid omega
1616 : ! closest to zero. Re is alpha(0); Im should be machine-zero (sanity).
1617 144 : IF (info_unit > 0) THEN
1618 868 : idx_omega_zero = MINLOC(ABS(omegas), DIM=1)
1619 38 : IF (i == 1) THEN
1620 : WRITE (info_unit, '(A,T22,A,T28,A,T36,A,T59,A)') &
1621 34 : " STATIC_POL|", "spin", "element", "Re [a.u.]", "Im [a.u.]"
1622 : END IF
1623 152 : DO k_static = 1, n_elems
1624 : WRITE (info_unit, '(A,T22,I4,T28,I3,",",I3,T36,ES22.10E3,T59,ES22.10E3)') &
1625 114 : " STATIC_POL|", i, &
1626 114 : rtc%print_pol_elements(k_static, 1), &
1627 114 : rtc%print_pol_elements(k_static, 2), &
1628 114 : REAL(pol_results(k_static, idx_omega_zero), kind=dp), &
1629 266 : AIMAG(pol_results(k_static, idx_omega_zero))
1630 : END DO
1631 : END IF
1632 : END DO
1633 : ! Spin-summed total polarizability (open shell only; inert for nspin=1).
1634 : ! Field is spin-independent, so (sum_s moments_s)/field == sum_s (moments_s/field).
1635 68 : IF (nspin > 1) THEN
1636 24 : ALLOCATE (pol_results_spin_total(n_elems, n))
1637 8 : pol_results_spin_total(:, :) = (0.0_dp, 0.0_dp)
1638 32 : DO k = 1, n_elems
1639 72 : DO i = 1, nspin
1640 : pol_results_spin_total(k, :) = pol_results_spin_total(k, :) + &
1641 960 : results(3*(i - 1) + rtc%print_pol_elements(k, 1), :)
1642 : END DO
1643 : pol_results_spin_total(k, :) = pol_results_spin_total(k, :)/ &
1644 : (field_results(rtc%print_pol_elements(k, 2), :) + &
1645 476 : 1.0e-10*field_results(rtc%print_pol_elements(k, 2), 2))
1646 : END DO
1647 : ft_unit = cp_print_key_unit_nr(logger, pol_section, extension="_SPIN_TOTAL.dat", &
1648 8 : file_form="FORMATTED", file_position="REWIND")
1649 8 : IF (ft_unit > 0) THEN
1650 12 : ALLOCATE (headers(2*n_elems + 1))
1651 16 : DO k = 1, n_elems
1652 12 : WRITE (headers(2*k), "(A16,I2,I2)") "real pol. elem.", &
1653 12 : rtc%print_pol_elements(k, 1), &
1654 24 : rtc%print_pol_elements(k, 2)
1655 12 : WRITE (headers(2*k + 1), "(A16,I2,I2)") "imag pol. elem.", &
1656 12 : rtc%print_pol_elements(k, 1), &
1657 28 : rtc%print_pol_elements(k, 2)
1658 : END DO
1659 4 : IF (info_unit == ft_unit) THEN
1660 0 : headers(1) = "# Energy [eV]"
1661 0 : prefix = " POLARIZABILITY|"
1662 0 : prefix_format = "(A16)"
1663 : CALL print_rt_file(ft_unit, headers, omegas, pol_results_spin_total, &
1664 0 : prefix, prefix_format, evolt)
1665 : ELSE
1666 4 : headers(1) = "# omega [at.u.]"
1667 4 : CALL print_rt_file(ft_unit, headers, omegas, pol_results_spin_total)
1668 : END IF
1669 4 : DEALLOCATE (headers)
1670 : END IF
1671 8 : CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1672 : ! Static polarizability total (header row already emitted by the per-spin block)
1673 8 : IF (info_unit > 0) THEN
1674 78 : idx_omega_zero = MINLOC(ABS(omegas), DIM=1)
1675 16 : DO k_static = 1, n_elems
1676 : WRITE (info_unit, '(A,T22,A,T28,I3,",",I3,T36,ES22.10E3,T59,ES22.10E3)') &
1677 12 : " STATIC_POL|", "TOT", &
1678 12 : rtc%print_pol_elements(k_static, 1), &
1679 12 : rtc%print_pol_elements(k_static, 2), &
1680 12 : REAL(pol_results_spin_total(k_static, idx_omega_zero), kind=dp), &
1681 28 : AIMAG(pol_results_spin_total(k_static, idx_omega_zero))
1682 : END DO
1683 : END IF
1684 8 : DEALLOCATE (pol_results_spin_total)
1685 : END IF
1686 : END IF
1687 :
1688 : ! Padé polarizability
1689 90 : IF (rtc%pade_requested .AND. do_polarizability) THEN
1690 : ! Start with the field pade
1691 6 : ALLOCATE (field_results_pade(3, n_pade))
1692 2 : IF (rtc%apply_delta_pulse) THEN
1693 8 : DO k = 1, 3
1694 6002 : field_results_pade(k, :) = CMPLX(delta_vec(k), 0.0, kind=dp)
1695 : END DO
1696 : ELSE
1697 0 : DO k = 1, 3
1698 : CALL greenx_refine_ft(rtc%pade_fit_e_min, rtc%pade_fit_e_max, &
1699 : omegas_complex, field_results(k, :), &
1700 0 : omegas_pade, field_results_pade(k, :))
1701 : END DO
1702 : END IF
1703 : ! Allocate polarisation pade
1704 8 : ALLOCATE (pol_results_pade(n_elems, n_pade))
1705 : ! Refine
1706 4 : DO i = 1, nspin
1707 8 : DO k = 1, n_elems
1708 : ! NOTE : Regularization to small value
1709 : pol_results_pade(k, :) = results_pade(3*(i - 1) + rtc%print_pol_elements(k, 1), :)/( &
1710 : field_results_pade(rtc%print_pol_elements(k, 2), :) + &
1711 6002 : field_results_pade(rtc%print_pol_elements(k, 2), 2)*1.0e-10_dp)
1712 : END DO
1713 : ! Print to the file
1714 : ft_unit = cp_print_key_unit_nr(logger, pol_section, extension="_PADE"//file_extensions(i), &
1715 2 : file_form="FORMATTED", file_position="REWIND")
1716 2 : IF (ft_unit > 0) THEN
1717 3 : ALLOCATE (headers(2*n_elems + 1))
1718 4 : DO k = 1, n_elems
1719 3 : WRITE (headers(2*k), "(A16,I2,I2)") "re,pade,pol.", &
1720 3 : rtc%print_pol_elements(k, 1), &
1721 6 : rtc%print_pol_elements(k, 2)
1722 3 : WRITE (headers(2*k + 1), "(A16,I2,I2)") "im,pade,pol.", &
1723 3 : rtc%print_pol_elements(k, 1), &
1724 7 : rtc%print_pol_elements(k, 2)
1725 : END DO
1726 : ! Write header
1727 1 : IF (info_unit == ft_unit) THEN
1728 1 : headers(1) = "# Energy [eV]"
1729 1 : prefix = " POLARIZABILITY_PADE|"
1730 1 : prefix_format = "(A21)"
1731 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade, &
1732 1 : prefix, prefix_format, evolt)
1733 : ELSE
1734 0 : headers(1) = "# omega [at.u.]"
1735 0 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade)
1736 : END IF
1737 1 : DEALLOCATE (headers)
1738 : END IF
1739 4 : CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1740 : END DO
1741 : ! Spin-summed total Padé polarizability (open shell only; inert for nspin=1)
1742 2 : IF (nspin > 1) THEN
1743 0 : ALLOCATE (pol_results_pade_spin_total(n_elems, n_pade))
1744 0 : pol_results_pade_spin_total(:, :) = (0.0_dp, 0.0_dp)
1745 0 : DO k = 1, n_elems
1746 0 : DO i = 1, nspin
1747 : pol_results_pade_spin_total(k, :) = pol_results_pade_spin_total(k, :) + &
1748 0 : results_pade(3*(i - 1) + rtc%print_pol_elements(k, 1), :)
1749 : END DO
1750 : pol_results_pade_spin_total(k, :) = pol_results_pade_spin_total(k, :)/( &
1751 : field_results_pade(rtc%print_pol_elements(k, 2), :) + &
1752 0 : field_results_pade(rtc%print_pol_elements(k, 2), 2)*1.0e-10_dp)
1753 : END DO
1754 : ft_unit = cp_print_key_unit_nr(logger, pol_section, extension="_PADE_SPIN_TOTAL.dat", &
1755 0 : file_form="FORMATTED", file_position="REWIND")
1756 0 : IF (ft_unit > 0) THEN
1757 0 : ALLOCATE (headers(2*n_elems + 1))
1758 0 : DO k = 1, n_elems
1759 0 : WRITE (headers(2*k), "(A16,I2,I2)") "re,pade,pol.", &
1760 0 : rtc%print_pol_elements(k, 1), &
1761 0 : rtc%print_pol_elements(k, 2)
1762 0 : WRITE (headers(2*k + 1), "(A16,I2,I2)") "im,pade,pol.", &
1763 0 : rtc%print_pol_elements(k, 1), &
1764 0 : rtc%print_pol_elements(k, 2)
1765 : END DO
1766 0 : IF (info_unit == ft_unit) THEN
1767 0 : headers(1) = "# Energy [eV]"
1768 0 : prefix = " POLARIZABILITY_PADE|"
1769 0 : prefix_format = "(A21)"
1770 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade_spin_total, &
1771 0 : prefix, prefix_format, evolt)
1772 : ELSE
1773 0 : headers(1) = "# omega [at.u.]"
1774 0 : CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade_spin_total)
1775 : END IF
1776 0 : DEALLOCATE (headers)
1777 : END IF
1778 0 : CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1779 0 : DEALLOCATE (pol_results_pade_spin_total)
1780 : END IF
1781 2 : DEALLOCATE (field_results_pade)
1782 2 : DEALLOCATE (pol_results_pade)
1783 : END IF
1784 :
1785 90 : IF (rtc%pade_requested .AND. (do_moments_ft .OR. do_polarizability)) THEN
1786 2 : DEALLOCATE (omegas_complex)
1787 2 : DEALLOCATE (omegas_pade)
1788 2 : DEALLOCATE (omegas_pade_real)
1789 2 : DEALLOCATE (results_pade)
1790 : END IF
1791 :
1792 90 : IF (do_polarizability) THEN
1793 68 : DEALLOCATE (pol_results)
1794 68 : DEALLOCATE (field_results)
1795 : END IF
1796 :
1797 90 : IF (do_moments_ft .OR. do_polarizability) THEN
1798 68 : DEALLOCATE (results)
1799 68 : DEALLOCATE (omegas)
1800 : END IF
1801 :
1802 90 : CALL timestop(handle)
1803 :
1804 90 : END SUBROUTINE print_ft
1805 :
1806 : ! **************************************************************************************************
1807 : !> \brief ...
1808 : !> \param rt_unit ...
1809 : !> \param headers ...
1810 : !> \param xvals ...
1811 : !> \param yvals ...
1812 : !> \param prefix ...
1813 : !> \param prefix_format ...
1814 : !> \param xscale_opt ...
1815 : !> \param comp_opt ...
1816 : ! **************************************************************************************************
1817 5114 : SUBROUTINE print_rt_file(rt_unit, headers, xvals, yvals, prefix, prefix_format, xscale_opt, comp_opt)
1818 : INTEGER, INTENT(IN) :: rt_unit
1819 : CHARACTER(len=20), DIMENSION(:), INTENT(IN), &
1820 : OPTIONAL :: headers
1821 : REAL(kind=dp), DIMENSION(:), INTENT(IN) :: xvals
1822 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: yvals
1823 : CHARACTER(len=21), INTENT(IN), OPTIONAL :: prefix
1824 : CHARACTER(len=5), INTENT(IN), OPTIONAL :: prefix_format
1825 : REAL(kind=dp), INTENT(IN), OPTIONAL :: xscale_opt
1826 : INTEGER, OPTIONAL :: comp_opt
1827 :
1828 : INTEGER :: do_comp, i, j, ncols, nrows
1829 : LOGICAL :: do_headers, do_prefix
1830 : REAL(kind=dp) :: xscale
1831 :
1832 5114 : do_prefix = .FALSE.
1833 5114 : IF (PRESENT(prefix)) THEN
1834 2260 : IF (PRESENT(prefix_format)) THEN
1835 : do_prefix = .TRUE.
1836 : ELSE
1837 0 : CPABORT("Printing of prefix with missing format!")
1838 : END IF
1839 : END IF
1840 :
1841 5114 : xscale = 1.0_dp
1842 5114 : IF (PRESENT(xscale_opt)) xscale = xscale_opt
1843 :
1844 5114 : ncols = SIZE(yvals, 1)
1845 5114 : nrows = SIZE(yvals, 2)
1846 :
1847 : ! Check whether printing complex data (default) or just a component
1848 5114 : do_comp = rt_file_comp_both
1849 5114 : IF (PRESENT(comp_opt)) do_comp = comp_opt
1850 :
1851 5114 : do_headers = PRESENT(headers)
1852 : ! Check whether enough headers for yvals and xvals is present
1853 5114 : IF (do_headers) THEN
1854 : IF ((do_comp /= rt_file_comp_both .AND. SIZE(headers) < ncols + 1) &
1855 216 : .OR. (do_comp == rt_file_comp_both .AND. SIZE(headers) < 2*ncols + 1)) THEN
1856 0 : CPABORT("Not enought headers to print the file!")
1857 : END IF
1858 : END IF
1859 :
1860 5114 : IF (SIZE(xvals) < nrows) THEN
1861 0 : CPABORT("Not enough xvals to print all yvals!")
1862 : END IF
1863 :
1864 5114 : IF (rt_unit > 0) THEN
1865 : ! Print headers
1866 672 : IF (do_headers) THEN
1867 : ! If prefix is present, write prefix
1868 94 : IF (do_prefix) THEN
1869 9 : WRITE (rt_unit, prefix_format, advance="no") prefix
1870 : END IF
1871 94 : WRITE (rt_unit, "(A20)", advance="no") headers(1)
1872 : ! Print the rest of the headers
1873 86 : SELECT CASE (do_comp)
1874 : CASE (rt_file_comp_both)
1875 : ! Complex case
1876 516 : DO j = 1, 2*ncols - 1
1877 516 : WRITE (rt_unit, "(A20)", advance="no") headers(j + 1)
1878 : END DO
1879 86 : WRITE (rt_unit, "(A20)") headers(2*ncols + 1)
1880 : CASE DEFAULT
1881 : ! For other cases, just one component is printed
1882 24 : DO j = 1, ncols - 1
1883 24 : WRITE (rt_unit, "(A20)", advance="no") headers(j + 1)
1884 : END DO
1885 102 : WRITE (rt_unit, "(A20)") headers(ncols + 1)
1886 : END SELECT
1887 : END IF
1888 : ! Done with the headers, print actual data
1889 5064 : DO i = 1, nrows
1890 : ! If prefix is present, write prefix
1891 4392 : IF (do_prefix) THEN
1892 1605 : WRITE (rt_unit, prefix_format, advance="no") prefix
1893 : END IF
1894 4392 : WRITE (rt_unit, "(E20.8E3)", advance="no") xvals(i)*xscale
1895 13176 : DO j = 1, ncols - 1
1896 4392 : SELECT CASE (do_comp)
1897 : CASE (rt_file_comp_real)
1898 : WRITE (rt_unit, "(E20.8E3)", advance="no") &
1899 586 : REAL(yvals(j, i))
1900 : CASE (rt_file_comp_imag)
1901 : WRITE (rt_unit, "(E20.8E3)", advance="no") &
1902 586 : AIMAG(yvals(j, i))
1903 : CASE DEFAULT
1904 : ! Print both components
1905 : WRITE (rt_unit, "(E20.8E3,E20.8E3)", advance="no") &
1906 8784 : REAL(yvals(j, i)), AIMAG(yvals(j, i))
1907 : END SELECT
1908 : END DO
1909 : ! Print the final column(s)
1910 672 : SELECT CASE (do_comp)
1911 : CASE (rt_file_comp_real)
1912 293 : WRITE (rt_unit, "(E20.8E3)") REAL(yvals(j, i))
1913 : CASE (rt_file_comp_imag)
1914 293 : WRITE (rt_unit, "(E20.8E3)") AIMAG(yvals(j, i))
1915 : CASE DEFAULT
1916 : ! Print both components
1917 : WRITE (rt_unit, "(E20.8E3,E20.8E3)") &
1918 4392 : REAL(yvals(j, i)), AIMAG(yvals(j, i))
1919 : END SELECT
1920 : END DO
1921 : END IF
1922 5114 : END SUBROUTINE print_rt_file
1923 :
1924 : END MODULE rt_propagation_output
|
