Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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
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 : #include "../base/base_uses.f90"
105 :
106 : IMPLICIT NONE
107 :
108 : PRIVATE
109 :
110 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_output'
111 :
112 : PUBLIC :: rt_prop_output, &
113 : rt_convergence, &
114 : rt_convergence_density, &
115 : report_density_occupation, &
116 : print_moments, &
117 : calc_local_moment, &
118 : print_ft
119 :
120 : CONTAINS
121 :
122 : ! **************************************************************************************************
123 : !> \brief ...
124 : !> \param qs_env ...
125 : !> \param run_type ...
126 : !> \param delta_iter ...
127 : !> \param used_time ...
128 : ! **************************************************************************************************
129 2284 : SUBROUTINE rt_prop_output(qs_env, run_type, delta_iter, used_time)
130 : TYPE(qs_environment_type), POINTER :: qs_env
131 : INTEGER, INTENT(in) :: run_type
132 : REAL(dp), INTENT(in), OPTIONAL :: delta_iter, used_time
133 :
134 : INTEGER :: i, n_electrons, n_proj, natom, nspin, &
135 : output_unit, spin, unit_nr
136 2284 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
137 2284 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
138 : LOGICAL :: new_file
139 : REAL(dp) :: orthonormality, strace, tot_rho_r, trace
140 : REAL(kind=dp), DIMENSION(3) :: field, reference_point
141 2284 : REAL(KIND=dp), DIMENSION(:), POINTER :: qs_tot_rho_r
142 2284 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: j_int
143 2284 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
144 2284 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
145 : TYPE(cp_logger_type), POINTER :: logger
146 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
147 2284 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, P_im, p_xyz, rho_new
148 : TYPE(dbcsr_type), POINTER :: tmp_ao
149 : TYPE(dft_control_type), POINTER :: dft_control
150 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
151 2284 : POINTER :: sab_all, sab_orb
152 2284 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
153 2284 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
154 : TYPE(qs_rho_type), POINTER :: rho
155 : TYPE(rt_prop_type), POINTER :: rtp
156 : TYPE(section_vals_type), POINTER :: dft_section, input, rtp_section
157 :
158 2284 : NULLIFY (logger, dft_control)
159 :
160 4568 : logger => cp_get_default_logger()
161 : CALL get_qs_env(qs_env, &
162 : rtp=rtp, &
163 : matrix_s=matrix_s, &
164 : input=input, &
165 : rho=rho, &
166 : particle_set=particle_set, &
167 : atomic_kind_set=atomic_kind_set, &
168 : qs_kind_set=qs_kind_set, &
169 : dft_control=dft_control, sab_all=sab_all, sab_orb=sab_orb, &
170 2284 : dbcsr_dist=dbcsr_dist)
171 :
172 2284 : rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
173 :
174 2284 : CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
175 2284 : n_electrons = n_electrons - dft_control%charge
176 :
177 2284 : CALL qs_rho_get(rho_struct=rho, tot_rho_r=qs_tot_rho_r)
178 :
179 2284 : tot_rho_r = accurate_sum(qs_tot_rho_r)
180 :
181 : output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
182 2284 : extension=".scfLog")
183 :
184 2284 : IF (output_unit > 0) THEN
185 : WRITE (output_unit, FMT="(/,(T3,A,T40,I5))") &
186 1142 : "Information at iteration step:", rtp%iter
187 : WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
188 1142 : "Total electronic density (r-space): ", &
189 1142 : tot_rho_r, &
190 : tot_rho_r + &
191 2284 : REAL(n_electrons, dp)
192 : WRITE (UNIT=output_unit, FMT="((T3,A,T59,F22.14))") &
193 1142 : "Total energy:", rtp%energy_new
194 1142 : IF (run_type == ehrenfest) &
195 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
196 572 : "Energy difference to previous iteration step:", rtp%energy_new - rtp%energy_old
197 1142 : IF (run_type == real_time_propagation) &
198 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
199 570 : "Energy difference to initial state:", rtp%energy_new - rtp%energy_old
200 1142 : IF (PRESENT(delta_iter)) &
201 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,E20.6))") &
202 1142 : "Convergence:", delta_iter
203 1142 : IF (rtp%converged) THEN
204 307 : IF (run_type == real_time_propagation) &
205 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F12.2))") &
206 172 : "Time needed for propagation:", used_time
207 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,3X,F16.14))") &
208 307 : "CONVERGENCE REACHED", rtp%energy_new - rtp%energy_old
209 : END IF
210 : END IF
211 :
212 2284 : IF (rtp%converged) THEN
213 614 : IF (.NOT. rtp%linear_scaling) THEN
214 432 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
215 : CALL rt_calculate_orthonormality(orthonormality, &
216 432 : mos_new, matrix_s(1)%matrix)
217 432 : IF (output_unit > 0) &
218 : WRITE (output_unit, FMT="(/,(T3,A,T60,F20.10))") &
219 216 : "Max deviation from orthonormalization:", orthonormality
220 : END IF
221 : END IF
222 :
223 2284 : IF (output_unit > 0) &
224 1142 : CALL m_flush(output_unit)
225 : CALL cp_print_key_finished_output(output_unit, logger, rtp_section, &
226 2284 : "PRINT%PROGRAM_RUN_INFO")
227 :
228 2284 : IF (rtp%converged) THEN
229 614 : dft_section => section_vals_get_subs_vals(input, "DFT")
230 614 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
231 : dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) &
232 30 : CALL print_field_applied(qs_env, dft_section)
233 614 : CALL make_moment(qs_env)
234 614 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
235 : dft_section, "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS"), cp_p_file)) THEN
236 4 : CALL print_rtp_energy_components(qs_env, dft_section)
237 : END IF
238 614 : IF (.NOT. dft_control%qs_control%dftb) THEN
239 494 : CALL write_available_results(qs_env=qs_env, rtp=rtp)
240 : END IF
241 :
242 614 : IF (rtp%linear_scaling) THEN
243 182 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
244 :
245 : ! Probably have to rebuild the moment matrix, since atoms can also move, in principle
246 182 : IF (dft_control%rtp_control%save_local_moments) THEN
247 : ! Save the field value
248 36 : IF (dft_control%apply_efield_field) THEN
249 0 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
250 0 : rtp%fields(:, rtp%istep + rtp%i_start + 1) = CMPLX(field(:), 0.0, kind=dp)
251 : END IF
252 36 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
253 0 : CALL build_local_moment_matrix(qs_env, rtp%local_moments, 1, reference_point)
254 : END IF
255 : ! TODO : Is symmetric rho possible?
256 : ! Spin + complex parts
257 : ! Extensions setup
258 : CALL calc_local_moment(rtp%local_moments, rho_new, &
259 36 : rtp%local_moments_work, rtp%moments(:, :, rtp%istep + rtp%i_start + 1))
260 : ! Time 1 is zero (start) time
261 36 : rtp%times(rtp%istep + rtp%i_start + 1) = qs_env%sim_time
262 36 : output_unit = cp_logger_get_default_io_unit(logger)
263 : CALL print_moments(section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS"), output_unit, &
264 36 : rtp%moments(:, :, rtp%istep + rtp%i_start + 1), qs_env%sim_time, rtp%track_imag_density)
265 : END IF
266 :
267 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
268 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART"), cp_p_file)) THEN
269 96 : CALL write_rt_p_to_restart(rho_new, .FALSE.)
270 : END IF
271 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
272 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY"), cp_p_file)) THEN
273 2 : CALL write_rt_p_to_restart(rho_new, .TRUE.)
274 : END IF
275 182 : IF (.NOT. dft_control%qs_control%dftb) THEN
276 : !Not sure if these things could also work with dftb or not
277 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
278 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
279 64 : DO spin = 1, SIZE(rho_new)/2
280 64 : CALL rt_current(qs_env, rho_new(2*spin)%matrix, dft_section, spin, SIZE(rho_new)/2)
281 : END DO
282 : END IF
283 : END IF
284 : ELSE
285 432 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
286 432 : IF (.NOT. dft_control%qs_control%dftb .AND. .NOT. dft_control%qs_control%xtb) THEN
287 252 : IF (rtp%track_imag_density) THEN
288 174 : NULLIFY (P_im, p_xyz)
289 174 : CALL dbcsr_allocate_matrix_set(p_xyz, 3)
290 :
291 : ! Linear momentum operator
292 : ! prepare for allocation
293 174 : natom = SIZE(particle_set, 1)
294 522 : ALLOCATE (first_sgf(natom))
295 348 : ALLOCATE (last_sgf(natom))
296 : CALL get_particle_set(particle_set, qs_kind_set, &
297 : first_sgf=first_sgf, &
298 174 : last_sgf=last_sgf)
299 348 : ALLOCATE (row_blk_sizes(natom))
300 174 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
301 174 : DEALLOCATE (first_sgf)
302 174 : DEALLOCATE (last_sgf)
303 :
304 174 : ALLOCATE (p_xyz(1)%matrix)
305 : CALL dbcsr_create(matrix=p_xyz(1)%matrix, &
306 : name="p_xyz", &
307 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
308 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
309 174 : mutable_work=.TRUE.)
310 174 : CALL cp_dbcsr_alloc_block_from_nbl(p_xyz(1)%matrix, sab_orb)
311 174 : CALL dbcsr_set(p_xyz(1)%matrix, 0.0_dp)
312 522 : DO i = 2, 3
313 348 : ALLOCATE (p_xyz(i)%matrix)
314 348 : CALL dbcsr_copy(p_xyz(i)%matrix, p_xyz(1)%matrix, "p_xyz-"//TRIM(ADJUSTL(cp_to_string(i))))
315 522 : CALL dbcsr_set(p_xyz(i)%matrix, 0.0_dp)
316 : END DO
317 174 : CALL build_lin_mom_matrix(qs_env, p_xyz)
318 174 : DEALLOCATE (row_blk_sizes)
319 :
320 174 : nspin = SIZE(mos_new)/2
321 174 : CALL qs_rho_get(rho, rho_ao_im=P_im)
322 522 : ALLOCATE (j_int(nspin, 3))
323 1290 : j_int = 0.0_dp
324 :
325 174 : NULLIFY (tmp_ao)
326 174 : CALL dbcsr_init_p(tmp_ao)
327 174 : CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
328 174 : CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
329 174 : CALL dbcsr_set(tmp_ao, 0.0_dp)
330 :
331 696 : DO i = 1, 3
332 522 : strace = 0.0_dp
333 1290 : DO spin = 1, nspin
334 594 : CALL dbcsr_set(tmp_ao, 0.0_dp)
335 : CALL dbcsr_multiply("T", "N", 1.0_dp, P_im(spin)%matrix, p_xyz(i)%matrix, &
336 594 : 0.0_dp, tmp_ao)
337 594 : CALL dbcsr_trace(tmp_ao, trace)
338 594 : strace = strace + trace
339 1116 : j_int(spin, i) = strace
340 : END DO
341 : END DO
342 : ! The term -(1/V)\int{(1/2)\sum_n f_n [\psi_n^*(t)(A(t))\psi_n(t) +cc]} missing. Is it just A(t)*number of electrons?
343 :
344 174 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
345 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT"), cp_p_file)) THEN
346 :
347 8 : output_unit = cp_logger_get_default_io_unit(logger)
348 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
349 8 : "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT", extension=".dat", is_new_file=new_file)
350 8 : IF (output_unit > 0) THEN
351 4 : IF (nspin == 2) THEN
352 2 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
353 4 : j_int(1, 1:3), j_int(2, 1:3)
354 : ELSE
355 2 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
356 4 : j_int(1, 1:3)
357 : END IF
358 : END IF
359 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
360 8 : "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT")
361 : END IF
362 174 : DEALLOCATE (j_int)
363 :
364 174 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
365 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
366 24 : DO spin = 1, nspin
367 24 : CALL rt_current(qs_env, P_im(spin)%matrix, dft_section, spin, nspin)
368 : END DO
369 : END IF
370 174 : CALL dbcsr_deallocate_matrix(tmp_ao)
371 174 : CALL dbcsr_deallocate_matrix_set(p_xyz)
372 : END IF
373 :
374 : ! projection of molecular orbitals
375 252 : IF (dft_control%rtp_control%is_proj_mo) THEN
376 112 : DO n_proj = 1, SIZE(dft_control%rtp_control%proj_mo_list)
377 : CALL compute_and_write_proj_mo(qs_env, mos_new, &
378 112 : dft_control%rtp_control%proj_mo_list(n_proj)%proj_mo, n_proj)
379 : END DO
380 : END IF
381 : END IF
382 : CALL write_rt_mos_to_restart(qs_env%mos, mos_new, particle_set, &
383 432 : dft_section, qs_kind_set)
384 : END IF
385 : END IF
386 :
387 2284 : rtp%energy_old = rtp%energy_new
388 :
389 2284 : IF (.NOT. rtp%converged .AND. rtp%iter >= dft_control%rtp_control%max_iter) &
390 : CALL cp_abort(__LOCATION__, "EMD did not converge, either increase MAX_ITER "// &
391 0 : "or use a smaller TIMESTEP")
392 :
393 4568 : END SUBROUTINE rt_prop_output
394 :
395 : ! **************************************************************************************************
396 : !> \brief computes the effective orthonormality of a set of mos given an s-matrix
397 : !> orthonormality is the max deviation from unity of the C^T S C
398 : !> \param orthonormality ...
399 : !> \param mos_new ...
400 : !> \param matrix_s ...
401 : !> \author Florian Schiffmann (02.09)
402 : ! **************************************************************************************************
403 432 : SUBROUTINE rt_calculate_orthonormality(orthonormality, mos_new, matrix_s)
404 : REAL(KIND=dp), INTENT(out) :: orthonormality
405 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
406 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
407 :
408 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_calculate_orthonormality'
409 :
410 : INTEGER :: handle, i, im, ispin, j, k, n, &
411 : ncol_local, nrow_local, nspin, re
412 432 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
413 : REAL(KIND=dp) :: alpha, max_alpha, max_beta
414 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
415 : TYPE(cp_fm_type) :: overlap_re, svec_im, svec_re
416 :
417 432 : NULLIFY (tmp_fm_struct)
418 :
419 432 : CALL timeset(routineN, handle)
420 :
421 432 : nspin = SIZE(mos_new)/2
422 432 : max_alpha = 0.0_dp
423 432 : max_beta = 0.0_dp
424 966 : DO ispin = 1, nspin
425 534 : re = ispin*2 - 1
426 534 : im = ispin*2
427 : ! get S*C
428 534 : CALL cp_fm_create(svec_re, mos_new(im)%matrix_struct)
429 534 : CALL cp_fm_create(svec_im, mos_new(im)%matrix_struct)
430 : CALL cp_fm_get_info(mos_new(im), &
431 534 : nrow_global=n, ncol_global=k)
432 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(re), &
433 534 : svec_re, k)
434 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(im), &
435 534 : svec_im, k)
436 :
437 : ! get C^T (S*C)
438 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
439 : para_env=mos_new(re)%matrix_struct%para_env, &
440 534 : context=mos_new(re)%matrix_struct%context)
441 534 : CALL cp_fm_create(overlap_re, tmp_fm_struct)
442 :
443 534 : CALL cp_fm_struct_release(tmp_fm_struct)
444 :
445 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(re), &
446 534 : svec_re, 0.0_dp, overlap_re)
447 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(im), &
448 534 : svec_im, 1.0_dp, overlap_re)
449 :
450 534 : CALL cp_fm_release(svec_re)
451 534 : CALL cp_fm_release(svec_im)
452 :
453 : CALL cp_fm_get_info(overlap_re, nrow_local=nrow_local, ncol_local=ncol_local, &
454 534 : row_indices=row_indices, col_indices=col_indices)
455 1413 : DO i = 1, nrow_local
456 5020 : DO j = 1, ncol_local
457 3607 : alpha = overlap_re%local_data(i, j)
458 3607 : IF (row_indices(i) == col_indices(j)) alpha = alpha - 1.0_dp
459 4486 : max_alpha = MAX(max_alpha, ABS(alpha))
460 : END DO
461 : END DO
462 2568 : CALL cp_fm_release(overlap_re)
463 : END DO
464 432 : CALL mos_new(1)%matrix_struct%para_env%max(max_alpha)
465 432 : CALL mos_new(1)%matrix_struct%para_env%max(max_beta)
466 432 : orthonormality = max_alpha
467 :
468 432 : CALL timestop(handle)
469 :
470 432 : END SUBROUTINE rt_calculate_orthonormality
471 :
472 : ! **************************************************************************************************
473 : !> \brief computes the convergence criterion for RTP and EMD
474 : !> \param rtp ...
475 : !> \param matrix_s Overlap matrix without the derivatives
476 : !> \param delta_mos ...
477 : !> \param delta_eps ...
478 : !> \author Florian Schiffmann (02.09)
479 : ! **************************************************************************************************
480 :
481 1478 : SUBROUTINE rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
482 : TYPE(rt_prop_type), POINTER :: rtp
483 : TYPE(dbcsr_type), POINTER :: matrix_s
484 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: delta_mos
485 : REAL(dp), INTENT(out) :: delta_eps
486 :
487 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence'
488 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
489 :
490 : INTEGER :: handle, i, icol, im, ispin, j, lcol, &
491 : lrow, nao, newdim, nmo, nspin, re
492 : LOGICAL :: double_col, double_row
493 : REAL(KIND=dp) :: alpha, max_alpha
494 : TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1, tmp_fm_struct
495 : TYPE(cp_fm_type) :: work, work1, work2
496 1478 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
497 :
498 1478 : NULLIFY (tmp_fm_struct)
499 :
500 1478 : CALL timeset(routineN, handle)
501 :
502 1478 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
503 :
504 1478 : nspin = SIZE(delta_mos)/2
505 1478 : max_alpha = 0.0_dp
506 :
507 5258 : DO i = 1, SIZE(mos_new)
508 5258 : CALL cp_fm_scale_and_add(-one, delta_mos(i), one, mos_new(i))
509 : END DO
510 :
511 3368 : DO ispin = 1, nspin
512 1890 : re = ispin*2 - 1
513 1890 : im = ispin*2
514 :
515 1890 : double_col = .TRUE.
516 1890 : double_row = .FALSE.
517 : CALL cp_fm_struct_double(newstruct, &
518 : delta_mos(re)%matrix_struct, &
519 : delta_mos(re)%matrix_struct%context, &
520 : double_col, &
521 1890 : double_row)
522 :
523 1890 : CALL cp_fm_create(work, matrix_struct=newstruct)
524 1890 : CALL cp_fm_create(work1, matrix_struct=newstruct)
525 :
526 : CALL cp_fm_get_info(delta_mos(re), ncol_local=lcol, ncol_global=nmo, &
527 1890 : nrow_global=nao)
528 1890 : CALL cp_fm_get_info(work, ncol_global=newdim)
529 :
530 1890 : CALL cp_fm_set_all(work, zero, zero)
531 :
532 8132 : DO icol = 1, lcol
533 51268 : work%local_data(:, icol) = delta_mos(re)%local_data(:, icol)
534 53158 : work%local_data(:, icol + lcol) = delta_mos(im)%local_data(:, icol)
535 : END DO
536 :
537 1890 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, work, work1, ncol=newdim)
538 :
539 1890 : CALL cp_fm_release(work)
540 :
541 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, &
542 : para_env=delta_mos(re)%matrix_struct%para_env, &
543 1890 : context=delta_mos(re)%matrix_struct%context)
544 : CALL cp_fm_struct_double(newstruct1, &
545 : tmp_fm_struct, &
546 : delta_mos(re)%matrix_struct%context, &
547 : double_col, &
548 1890 : double_row)
549 :
550 1890 : CALL cp_fm_create(work, matrix_struct=newstruct1)
551 1890 : CALL cp_fm_create(work2, matrix_struct=newstruct1)
552 :
553 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(re), &
554 1890 : work1, zero, work)
555 :
556 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(im), &
557 1890 : work1, zero, work2)
558 :
559 1890 : CALL cp_fm_get_info(work, nrow_local=lrow)
560 5011 : DO i = 1, lrow
561 17796 : DO j = 1, lcol
562 : alpha = SQRT((work%local_data(i, j) + work2%local_data(i, j + lcol))**2 + &
563 12785 : (work%local_data(i, j + lcol) - work2%local_data(i, j))**2)
564 15906 : max_alpha = MAX(max_alpha, ABS(alpha))
565 : END DO
566 : END DO
567 :
568 1890 : CALL cp_fm_release(work)
569 1890 : CALL cp_fm_release(work1)
570 1890 : CALL cp_fm_release(work2)
571 1890 : CALL cp_fm_struct_release(tmp_fm_struct)
572 1890 : CALL cp_fm_struct_release(newstruct)
573 9038 : CALL cp_fm_struct_release(newstruct1)
574 :
575 : END DO
576 :
577 1478 : CALL delta_mos(1)%matrix_struct%para_env%max(max_alpha)
578 1478 : delta_eps = SQRT(max_alpha)
579 :
580 1478 : CALL timestop(handle)
581 :
582 1478 : END SUBROUTINE rt_convergence
583 :
584 : ! **************************************************************************************************
585 : !> \brief computes the convergence criterion for RTP and EMD based on the density matrix
586 : !> \param rtp ...
587 : !> \param delta_P ...
588 : !> \param delta_eps ...
589 : !> \author Samuel Andermatt (02.14)
590 : ! **************************************************************************************************
591 :
592 1612 : SUBROUTINE rt_convergence_density(rtp, delta_P, delta_eps)
593 :
594 : TYPE(rt_prop_type), POINTER :: rtp
595 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P
596 : REAL(dp), INTENT(out) :: delta_eps
597 :
598 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence_density'
599 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
600 :
601 : INTEGER :: col_atom, handle, i, ispin, row_atom
602 : REAL(dp) :: alpha, max_alpha
603 806 : REAL(dp), DIMENSION(:, :), POINTER :: block_values
604 : TYPE(dbcsr_iterator_type) :: iter
605 806 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
606 : TYPE(dbcsr_type), POINTER :: tmp
607 : TYPE(mp_comm_type) :: group
608 :
609 806 : CALL timeset(routineN, handle)
610 :
611 806 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
612 :
613 2910 : DO i = 1, SIZE(rho_new)
614 2910 : CALL dbcsr_add(delta_P(i)%matrix, rho_new(i)%matrix, one, -one)
615 : END DO
616 : !get the maximum value of delta_P
617 2910 : DO i = 1, SIZE(delta_P)
618 : !square all entries of both matrices
619 2104 : CALL dbcsr_iterator_start(iter, delta_P(i)%matrix)
620 12246 : DO WHILE (dbcsr_iterator_blocks_left(iter))
621 10142 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
622 788738 : block_values = block_values*block_values
623 : END DO
624 5014 : CALL dbcsr_iterator_stop(iter)
625 : END DO
626 : NULLIFY (tmp)
627 806 : ALLOCATE (tmp)
628 806 : CALL dbcsr_create(tmp, template=delta_P(1)%matrix, matrix_type="N")
629 1858 : DO ispin = 1, SIZE(delta_P)/2
630 1052 : CALL dbcsr_desymmetrize(delta_P(2*ispin - 1)%matrix, tmp)
631 1858 : CALL dbcsr_add(delta_P(2*ispin)%matrix, tmp, one, one)
632 : END DO
633 : !the absolute values are now in the even entries of delta_P
634 806 : max_alpha = zero
635 1858 : DO ispin = 1, SIZE(delta_P)/2
636 1052 : CALL dbcsr_iterator_start(iter, delta_P(2*ispin)%matrix)
637 6220 : DO WHILE (dbcsr_iterator_blocks_left(iter))
638 5168 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
639 398828 : alpha = MAXVAL(block_values)
640 6220 : IF (alpha > max_alpha) max_alpha = alpha
641 : END DO
642 2910 : CALL dbcsr_iterator_stop(iter)
643 : END DO
644 806 : CALL dbcsr_get_info(delta_P(1)%matrix, group=group)
645 806 : CALL group%max(max_alpha)
646 806 : delta_eps = SQRT(max_alpha)
647 806 : CALL dbcsr_deallocate_matrix(tmp)
648 806 : CALL timestop(handle)
649 :
650 806 : END SUBROUTINE rt_convergence_density
651 :
652 : ! **************************************************************************************************
653 : !> \brief interface to qs_moments. Does only work for nonperiodic dipole
654 : !> \param qs_env ...
655 : !> \author Florian Schiffmann (02.09)
656 : ! **************************************************************************************************
657 :
658 614 : SUBROUTINE make_moment(qs_env)
659 :
660 : TYPE(qs_environment_type), POINTER :: qs_env
661 :
662 : CHARACTER(len=*), PARAMETER :: routineN = 'make_moment'
663 :
664 : INTEGER :: handle, output_unit
665 : TYPE(cp_logger_type), POINTER :: logger
666 : TYPE(dft_control_type), POINTER :: dft_control
667 :
668 614 : CALL timeset(routineN, handle)
669 :
670 614 : NULLIFY (dft_control)
671 :
672 614 : logger => cp_get_default_logger()
673 614 : output_unit = cp_logger_get_default_io_unit(logger)
674 614 : CALL get_qs_env(qs_env, dft_control=dft_control)
675 614 : IF (dft_control%qs_control%dftb) THEN
676 120 : CALL scf_post_calculation_tb(qs_env, "DFTB", .FALSE.)
677 494 : ELSE IF (dft_control%qs_control%xtb) THEN
678 60 : CALL scf_post_calculation_tb(qs_env, "xTB", .FALSE.)
679 : ELSE
680 434 : CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
681 : END IF
682 614 : CALL timestop(handle)
683 :
684 614 : END SUBROUTINE make_moment
685 :
686 : ! **************************************************************************************************
687 : !> \brief Reports the sparsity pattern of the complex density matrix
688 : !> \param filter_eps ...
689 : !> \param rho ...
690 : !> \author Samuel Andermatt (09.14)
691 : ! **************************************************************************************************
692 :
693 182 : SUBROUTINE report_density_occupation(filter_eps, rho)
694 :
695 : REAL(KIND=dp) :: filter_eps
696 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho
697 :
698 : CHARACTER(len=*), PARAMETER :: routineN = 'report_density_occupation'
699 :
700 : INTEGER :: handle, i, im, ispin, re, unit_nr
701 : REAL(KIND=dp) :: eps, occ
702 : TYPE(cp_logger_type), POINTER :: logger
703 182 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmp
704 :
705 182 : CALL timeset(routineN, handle)
706 :
707 182 : logger => cp_get_default_logger()
708 182 : unit_nr = cp_logger_get_default_io_unit(logger)
709 182 : NULLIFY (tmp)
710 182 : CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
711 686 : DO i = 1, SIZE(rho)
712 504 : CALL dbcsr_init_p(tmp(i)%matrix)
713 504 : CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
714 686 : CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
715 : END DO
716 434 : DO ispin = 1, SIZE(rho)/2
717 252 : re = 2*ispin - 1
718 252 : im = 2*ispin
719 252 : eps = MAX(filter_eps, 1.0E-11_dp)
720 2338 : DO WHILE (eps < 1.1_dp)
721 2086 : CALL dbcsr_filter(tmp(re)%matrix, eps)
722 2086 : occ = dbcsr_get_occupation(tmp(re)%matrix)
723 3129 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
724 2086 : ispin, " eps ", eps, " real: ", occ
725 2086 : eps = eps*10
726 : END DO
727 252 : eps = MAX(filter_eps, 1.0E-11_dp)
728 2520 : DO WHILE (eps < 1.1_dp)
729 2086 : CALL dbcsr_filter(tmp(im)%matrix, eps)
730 2086 : occ = dbcsr_get_occupation(tmp(im)%matrix)
731 3129 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
732 2086 : ispin, " eps ", eps, " imag: ", occ
733 2086 : eps = eps*10.0_dp
734 : END DO
735 : END DO
736 182 : CALL dbcsr_deallocate_matrix_set(tmp)
737 182 : CALL timestop(handle)
738 :
739 182 : END SUBROUTINE report_density_occupation
740 :
741 : ! **************************************************************************************************
742 : !> \brief Writes the density matrix and the atomic positions to a restart file
743 : !> \param rho_new ...
744 : !> \param history ...
745 : !> \author Samuel Andermatt (09.14)
746 : ! **************************************************************************************************
747 :
748 98 : SUBROUTINE write_rt_p_to_restart(rho_new, history)
749 :
750 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
751 : LOGICAL :: history
752 :
753 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_p_to_restart'
754 :
755 : CHARACTER(LEN=default_path_length) :: file_name, project_name
756 : INTEGER :: handle, im, ispin, re, unit_nr
757 : REAL(KIND=dp) :: cs_pos
758 : TYPE(cp_logger_type), POINTER :: logger
759 :
760 98 : CALL timeset(routineN, handle)
761 98 : logger => cp_get_default_logger()
762 98 : IF (logger%para_env%is_source()) THEN
763 49 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
764 : ELSE
765 : unit_nr = -1
766 : END IF
767 :
768 98 : project_name = logger%iter_info%project_name
769 234 : DO ispin = 1, SIZE(rho_new)/2
770 136 : re = 2*ispin - 1
771 136 : im = 2*ispin
772 136 : IF (history) THEN
773 : WRITE (file_name, '(A,I0,A)') &
774 2 : TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
775 : ELSE
776 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
777 : END IF
778 136 : cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.TRUE.)
779 136 : IF (unit_nr > 0) THEN
780 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
781 : END IF
782 136 : CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
783 136 : IF (history) THEN
784 : WRITE (file_name, '(A,I0,A)') &
785 2 : TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
786 : ELSE
787 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
788 : END IF
789 136 : cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.TRUE.)
790 136 : IF (unit_nr > 0) THEN
791 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
792 : END IF
793 234 : CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
794 : END DO
795 :
796 98 : CALL timestop(handle)
797 :
798 98 : END SUBROUTINE write_rt_p_to_restart
799 :
800 : ! **************************************************************************************************
801 : !> \brief Collocation of the current and printing of it in a cube file
802 : !> \param qs_env ...
803 : !> \param P_im ...
804 : !> \param dft_section ...
805 : !> \param spin ...
806 : !> \param nspin ...
807 : !> \author Samuel Andermatt (06.15)
808 : ! **************************************************************************************************
809 48 : SUBROUTINE rt_current(qs_env, P_im, dft_section, spin, nspin)
810 : TYPE(qs_environment_type), POINTER :: qs_env
811 : TYPE(dbcsr_type), POINTER :: P_im
812 : TYPE(section_vals_type), POINTER :: dft_section
813 : INTEGER :: spin, nspin
814 :
815 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_current'
816 :
817 : CHARACTER(len=1) :: char_spin
818 : CHARACTER(len=14) :: ext
819 : CHARACTER(len=2) :: sdir
820 : INTEGER :: dir, handle, print_unit
821 48 : INTEGER, DIMENSION(:), POINTER :: stride(:)
822 : LOGICAL :: mpi_io
823 : TYPE(cp_logger_type), POINTER :: logger
824 : TYPE(current_env_type) :: current_env
825 : TYPE(dbcsr_type), POINTER :: tmp, zero
826 : TYPE(particle_list_type), POINTER :: particles
827 : TYPE(pw_c1d_gs_type) :: gs
828 : TYPE(pw_env_type), POINTER :: pw_env
829 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
830 : TYPE(pw_r3d_rs_type) :: rs
831 : TYPE(qs_subsys_type), POINTER :: subsys
832 :
833 48 : CALL timeset(routineN, handle)
834 :
835 48 : logger => cp_get_default_logger()
836 48 : CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
837 48 : CALL qs_subsys_get(subsys, particles=particles)
838 48 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
839 :
840 48 : NULLIFY (zero, tmp)
841 48 : ALLOCATE (zero, tmp)
842 48 : CALL dbcsr_create(zero, template=P_im)
843 48 : CALL dbcsr_copy(zero, P_im)
844 48 : CALL dbcsr_set(zero, 0.0_dp)
845 48 : CALL dbcsr_create(tmp, template=P_im)
846 48 : CALL dbcsr_copy(tmp, P_im)
847 48 : IF (nspin == 1) THEN
848 32 : CALL dbcsr_scale(tmp, 0.5_dp)
849 : END IF
850 48 : current_env%gauge = -1
851 48 : current_env%gauge_init = .FALSE.
852 48 : CALL auxbas_pw_pool%create_pw(rs)
853 48 : CALL auxbas_pw_pool%create_pw(gs)
854 :
855 : NULLIFY (stride)
856 48 : ALLOCATE (stride(3))
857 :
858 192 : DO dir = 1, 3
859 :
860 144 : CALL pw_zero(rs)
861 144 : CALL pw_zero(gs)
862 :
863 144 : CALL calculate_jrho_resp(zero, tmp, zero, zero, dir, dir, rs, gs, qs_env, current_env, retain_rsgrid=.TRUE.)
864 :
865 576 : stride = section_get_ivals(dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT%STRIDE")
866 :
867 144 : IF (dir == 1) THEN
868 48 : sdir = "-x"
869 96 : ELSEIF (dir == 2) THEN
870 48 : sdir = "-y"
871 : ELSE
872 48 : sdir = "-z"
873 : END IF
874 144 : WRITE (char_spin, "(I1)") spin
875 :
876 144 : ext = "-SPIN-"//char_spin//sdir//".cube"
877 144 : mpi_io = .TRUE.
878 : print_unit = cp_print_key_unit_nr(logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
879 : extension=ext, file_status="REPLACE", file_action="WRITE", &
880 144 : log_filename=.FALSE., mpi_io=mpi_io)
881 :
882 : CALL cp_pw_to_cube(rs, print_unit, "EMD current", particles=particles, stride=stride, &
883 144 : mpi_io=mpi_io)
884 :
885 : CALL cp_print_key_finished_output(print_unit, logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
886 192 : mpi_io=mpi_io)
887 :
888 : END DO
889 :
890 48 : CALL auxbas_pw_pool%give_back_pw(rs)
891 48 : CALL auxbas_pw_pool%give_back_pw(gs)
892 :
893 48 : CALL dbcsr_deallocate_matrix(zero)
894 48 : CALL dbcsr_deallocate_matrix(tmp)
895 :
896 48 : DEALLOCATE (stride)
897 :
898 48 : CALL timestop(handle)
899 :
900 3504 : END SUBROUTINE rt_current
901 :
902 : ! **************************************************************************************************
903 : !> \brief Interface routine to trigger writing of results available from normal
904 : !> SCF. Can write MO-dependent and MO free results (needed for call from
905 : !> the linear scaling code)
906 : !> Update: trigger also some of prints for time-dependent runs
907 : !> \param qs_env ...
908 : !> \param rtp ...
909 : !> \par History
910 : !> 2022-11 Update [Guillaume Le Breton]
911 : ! **************************************************************************************************
912 494 : SUBROUTINE write_available_results(qs_env, rtp)
913 : TYPE(qs_environment_type), POINTER :: qs_env
914 : TYPE(rt_prop_type), POINTER :: rtp
915 :
916 : CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results'
917 :
918 : INTEGER :: handle
919 : TYPE(qs_scf_env_type), POINTER :: scf_env
920 :
921 494 : CALL timeset(routineN, handle)
922 :
923 494 : CALL get_qs_env(qs_env, scf_env=scf_env)
924 494 : IF (rtp%linear_scaling) THEN
925 182 : CALL write_mo_free_results(qs_env)
926 : ELSE
927 312 : CALL write_mo_free_results(qs_env)
928 312 : CALL write_mo_dependent_results(qs_env, scf_env)
929 : ! Time-dependent MO print
930 312 : CALL write_rtp_mos_to_output_unit(qs_env, rtp)
931 312 : CALL write_rtp_mo_cubes(qs_env, rtp)
932 : END IF
933 :
934 494 : CALL timestop(handle)
935 :
936 494 : END SUBROUTINE write_available_results
937 :
938 : ! **************************************************************************************************
939 : !> \brief Print the field applied to the system. Either the electric
940 : !> field or the vector potential depending on the gauge used
941 : !> \param qs_env ...
942 : !> \param dft_section ...
943 : !> \par History
944 : !> 2023-01 Created [Guillaume Le Breton]
945 : ! **************************************************************************************************
946 30 : SUBROUTINE print_field_applied(qs_env, dft_section)
947 : TYPE(qs_environment_type), POINTER :: qs_env
948 : TYPE(section_vals_type), POINTER :: dft_section
949 :
950 : CHARACTER(LEN=3), DIMENSION(3) :: rlab
951 : CHARACTER(LEN=default_path_length) :: filename
952 : INTEGER :: i, i_step, output_unit, unit_nr
953 : LOGICAL :: new_file
954 : REAL(kind=dp) :: field(3)
955 : TYPE(cp_logger_type), POINTER :: logger
956 : TYPE(dft_control_type), POINTER :: dft_control
957 : TYPE(rt_prop_type), POINTER :: rtp
958 :
959 30 : NULLIFY (dft_control)
960 :
961 30 : logger => cp_get_default_logger()
962 30 : output_unit = cp_logger_get_default_io_unit(logger)
963 :
964 30 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp)
965 :
966 30 : i_step = rtp%istep
967 :
968 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
969 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat", is_new_file=new_file)
970 :
971 30 : IF (output_unit > 0) THEN
972 60 : rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
973 15 : IF (unit_nr /= output_unit) THEN
974 15 : INQUIRE (UNIT=unit_nr, NAME=filename)
975 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
976 15 : "FIELD", "The field applied is written to the file:", &
977 30 : TRIM(filename)
978 : ELSE
979 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED [a.u.]"
980 : WRITE (UNIT=output_unit, FMT="(T5,3(A,A,E16.8,1X))") &
981 0 : (TRIM(rlab(i)), "=", dft_control%rtp_control%field(i), i=1, 3)
982 : END IF
983 :
984 15 : IF (new_file) THEN
985 3 : IF (dft_control%apply_efield_field) THEN
986 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,3(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z"
987 2 : ELSE IF (dft_control%apply_vector_potential) THEN
988 0 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,6(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z", &
989 0 : " Vec. Pot. X", " Vec. Pot. Y", " Vec. Pot. Z"
990 : END IF
991 : END IF
992 :
993 15 : field = 0.0_dp
994 15 : IF (dft_control%apply_efield_field) THEN
995 4 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
996 4 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
997 8 : field(1), field(2), field(3)
998 : ! DO i=1,3
999 : ! IF (ABS(field(i))< 10E-10) field(i) = 0.0_dp
1000 : ! END IF
1001 11 : ELSE IF (dft_control%apply_vector_potential) THEN
1002 9 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1003 9 : dft_control%rtp_control%field(1), dft_control%rtp_control%field(2), dft_control%rtp_control%field(3), &
1004 18 : dft_control%rtp_control%vec_pot(1), dft_control%rtp_control%vec_pot(2), dft_control%rtp_control%vec_pot(3)
1005 : END IF
1006 :
1007 : END IF
1008 :
1009 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
1010 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD")
1011 :
1012 30 : END SUBROUTINE print_field_applied
1013 :
1014 : ! **************************************************************************************************
1015 : !> \brief Print the components of the total energy used in an RTP calculation
1016 : !> \param qs_env ...
1017 : !> \param dft_section ...
1018 : !> \par History
1019 : !> 2024-02 Created [ANB]
1020 : ! **************************************************************************************************
1021 4 : SUBROUTINE print_rtp_energy_components(qs_env, dft_section)
1022 : TYPE(qs_environment_type), POINTER :: qs_env
1023 : TYPE(section_vals_type), POINTER :: dft_section
1024 :
1025 : CHARACTER(LEN=default_path_length) :: filename
1026 : INTEGER :: i_step, output_unit, unit_nr
1027 : LOGICAL :: new_file
1028 : TYPE(cp_logger_type), POINTER :: logger
1029 : TYPE(dft_control_type), POINTER :: dft_control
1030 : TYPE(qs_energy_type), POINTER :: energy
1031 : TYPE(rt_prop_type), POINTER :: rtp
1032 :
1033 4 : NULLIFY (dft_control, energy, rtp)
1034 :
1035 4 : logger => cp_get_default_logger()
1036 4 : output_unit = cp_logger_get_default_io_unit(logger)
1037 :
1038 4 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
1039 4 : i_step = rtp%istep
1040 :
1041 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
1042 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS", extension=".ener", &
1043 4 : file_action="WRITE", is_new_file=new_file)
1044 :
1045 4 : IF (output_unit > 0) THEN
1046 2 : IF (unit_nr /= output_unit) THEN
1047 2 : INQUIRE (UNIT=unit_nr, NAME=filename)
1048 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
1049 2 : "ENERGY_CONSTITUENTS", "Total Energy constituents written to file:", &
1050 4 : TRIM(filename)
1051 : ELSE
1052 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ENERGY_CONSTITUENTS"
1053 : END IF
1054 :
1055 2 : IF (new_file) THEN
1056 : ! NOTE that these are not all terms contributing to the total energy for RTP, only a selection of those
1057 : ! most significant / impactful. Therefore the printed components likely will not add up to the total energy.
1058 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,10(6X,A))') "Step Nr.", "Time[fs]", &
1059 1 : "Total ener.[a.u.]", "core[a.u.] ", " overlap [a.u.]", "hartree[a.u.]", " exc. [a.u.] ", &
1060 2 : " hartree 1c[a.u.]", "exc 1c[a.u.] ", "exc admm[a.u.]", "exc 1c admm[a.u.]", "efield LG"
1061 :
1062 : END IF
1063 : WRITE (UNIT=unit_nr, FMT="(I10,F20.6,10(F20.9))") &
1064 2 : qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1065 2 : energy%total, energy%core, energy%core_overlap, energy%hartree, energy%exc, &
1066 4 : energy%hartree_1c, energy%exc1, energy%exc_aux_fit, energy%exc1_aux_fit, energy%efield_core
1067 :
1068 : END IF
1069 :
1070 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
1071 4 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS")
1072 :
1073 4 : END SUBROUTINE print_rtp_energy_components
1074 :
1075 : ! **************************************************************************************************
1076 : !> \brief Print the dipole moments into a file
1077 : !> \param moments_section Section of the input defining the file/stream to print the moments to
1078 : !> \param info_unit Unit where standard output from the program is written - for add. identifiers
1079 : !> \param moments Actual moment values (for specific time step)
1080 : !> \param time Current simulation time
1081 : !> \param imag_opt Whether to calculate the imaginary part
1082 : !> \par History
1083 : !> 10.2025 Created [Marek]
1084 : ! **************************************************************************************************
1085 616 : SUBROUTINE print_moments(moments_section, info_unit, moments, time, imag_opt)
1086 : TYPE(section_vals_type), POINTER :: moments_section
1087 : INTEGER :: info_unit
1088 : COMPLEX(kind=dp), DIMENSION(:, :) :: moments
1089 : REAL(kind=dp), OPTIONAL :: time
1090 : LOGICAL, OPTIONAL :: imag_opt
1091 :
1092 : CHARACTER(len=14), DIMENSION(4) :: file_extensions
1093 : INTEGER :: i, j, ndir, nspin, print_unit
1094 : LOGICAL :: imaginary
1095 : TYPE(cp_logger_type), POINTER :: logger
1096 :
1097 : ! Index 1 : spin, Index 2 : direction
1098 :
1099 616 : nspin = SIZE(moments, 1)
1100 616 : ndir = SIZE(moments, 2)
1101 :
1102 616 : IF (nspin < 1) CPABORT("Zero spin index size in print moments!")
1103 616 : IF (ndir < 1) CPABORT("Zero direction index size in print moments!")
1104 :
1105 616 : imaginary = .TRUE.
1106 616 : IF (PRESENT(imag_opt)) imaginary = imag_opt
1107 :
1108 : ! Get the program run info unit and target unit
1109 : ! If these are the same (most likely the case of __STD_OUT__), add
1110 : ! extra identifier to the printed output
1111 616 : file_extensions(1) = "_SPIN_A_RE.dat"
1112 616 : file_extensions(2) = "_SPIN_A_IM.dat"
1113 616 : file_extensions(3) = "_SPIN_B_RE.dat"
1114 616 : file_extensions(4) = "_SPIN_B_IM.dat"
1115 616 : logger => cp_get_default_logger()
1116 1268 : DO i = 1, nspin
1117 : ! Real part
1118 652 : print_unit = cp_print_key_unit_nr(logger, moments_section, extension=file_extensions(2*i - 1))
1119 652 : IF (print_unit == info_unit .AND. print_unit > 0) THEN
1120 : ! Do the output with additional suffix
1121 211 : WRITE (print_unit, "(A18)", advance="no") " MOMENTS_TRACE_RE|"
1122 : END IF
1123 652 : IF (print_unit > 0) THEN
1124 230 : IF (PRESENT(time)) WRITE (print_unit, "(E20.8E3)", advance="no") time*femtoseconds
1125 690 : DO j = 1, ndir - 1
1126 690 : WRITE (print_unit, "(E20.8E3)", advance="no") REAL(moments(i, j))
1127 : END DO
1128 : ! Write the last direction
1129 230 : WRITE (print_unit, "(E20.8E3)") REAL(moments(i, ndir))
1130 : END IF
1131 652 : CALL cp_print_key_finished_output(print_unit, logger, moments_section)
1132 : ! Same for imaginary part
1133 1268 : IF (imaginary) THEN
1134 562 : print_unit = cp_print_key_unit_nr(logger, moments_section, extension=file_extensions(2*i))
1135 562 : IF (print_unit == info_unit .AND. print_unit > 0) THEN
1136 : ! Do the output with additional suffix
1137 211 : WRITE (print_unit, "(A18)", advance="no") " MOMENTS_TRACE_IM|"
1138 : END IF
1139 562 : IF (print_unit > 0) THEN
1140 230 : IF (PRESENT(time)) WRITE (print_unit, "(E20.8E3)", advance="no") time*femtoseconds
1141 690 : DO j = 1, ndir - 1
1142 690 : WRITE (print_unit, "(E20.8E3)", advance="no") AIMAG(moments(i, j))
1143 : END DO
1144 : ! Write the last direction
1145 230 : WRITE (print_unit, "(E20.8E3)") AIMAG(moments(i, ndir))
1146 : END IF
1147 562 : CALL cp_print_key_finished_output(print_unit, logger, moments_section)
1148 : END IF
1149 : END DO
1150 :
1151 616 : END SUBROUTINE print_moments
1152 :
1153 : ! **************************************************************************************************
1154 : !> \brief Calculate the values of real/imaginary parts of moments in all directions
1155 : !> \param moment_matrices Local matrix representations of dipole (position) operator
1156 : !> \param density_matrices Density matrices (spin and real+complex parts)
1157 : !> \param work Extra dbcsr matrix for work
1158 : !> \param moment Resulting moments (spin and direction)
1159 : !> \param imag_opt Whether to calculate the imaginary part of the moment
1160 : !> \par History
1161 : !> 10.2025 Created [Marek]
1162 : ! **************************************************************************************************
1163 54 : SUBROUTINE calc_local_moment(moment_matrices, density_matrices, work, moment, imag_opt)
1164 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moment_matrices, density_matrices
1165 : TYPE(dbcsr_type) :: work
1166 : COMPLEX(kind=dp), DIMENSION(:, :) :: moment
1167 : LOGICAL, OPTIONAL :: imag_opt
1168 :
1169 : INTEGER :: i, k, nspin
1170 : LOGICAL :: imag
1171 : REAL(kind=dp) :: real_moment
1172 :
1173 54 : imag = .FALSE.
1174 54 : IF (PRESENT(imag_opt)) imag = imag_opt
1175 54 : nspin = SIZE(density_matrices)/2
1176 :
1177 144 : DO i = 1, nspin
1178 414 : DO k = 1, 3
1179 : CALL dbcsr_multiply("N", "N", -1.0_dp, &
1180 : density_matrices(2*i - 1)%matrix, moment_matrices(k)%matrix, &
1181 270 : 0.0_dp, work)
1182 270 : CALL dbcsr_trace(work, real_moment)
1183 270 : moment(i, k) = CMPLX(real_moment, 0.0, kind=dp)
1184 360 : IF (imag) THEN
1185 : CALL dbcsr_multiply("N", "N", -1.0_dp, &
1186 : density_matrices(2*i)%matrix, moment_matrices(k)%matrix, &
1187 0 : 0.0_dp, work)
1188 0 : CALL dbcsr_trace(work, real_moment)
1189 0 : moment(i, k) = moment(i, k) + CMPLX(0.0, real_moment, kind=dp)
1190 : END IF
1191 : END DO
1192 : END DO
1193 :
1194 54 : END SUBROUTINE calc_local_moment
1195 :
1196 : ! **************************************************************************************************
1197 : !> \brief Calculate and print the Fourier transforms + polarizabilites from moment trace
1198 : !> \param rtp_section The RTP input section (needed to access PRINT configurations)
1199 : !> \param moments Moment trace
1200 : !> \param times Corresponding times
1201 : !> \param fields Corresponding fields
1202 : !> \param rtc rt_control_type that includes metadata
1203 : !> \param info_opt ...
1204 : !> \param cell If present, used to change the delta peak representation to be in units of reciprocal lattice
1205 : !> \par History
1206 : !> 10.2025 Created [Marek]
1207 : ! **************************************************************************************************
1208 30 : SUBROUTINE print_ft(rtp_section, moments, times, fields, rtc, info_opt, cell)
1209 : TYPE(section_vals_type), POINTER :: rtp_section
1210 : COMPLEX(kind=dp), DIMENSION(:, :, :), POINTER :: moments
1211 : REAL(kind=dp), DIMENSION(:), POINTER :: times
1212 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: fields
1213 : TYPE(rtp_control_type), POINTER :: rtc
1214 : INTEGER, OPTIONAL :: info_opt
1215 : TYPE(cell_type), OPTIONAL, POINTER :: cell
1216 :
1217 : CHARACTER(len=11), DIMENSION(2) :: file_extensions
1218 30 : CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: headers
1219 : CHARACTER(len=21) :: prefix
1220 : CHARACTER(len=5) :: prefix_format
1221 30 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: omegas_complex, omegas_pade
1222 30 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: field_results, field_results_pade, &
1223 30 : pol_results, pol_results_pade, &
1224 30 : results, results_pade, value_series
1225 : INTEGER :: ft_unit, i, info_unit, k, n, n_elems, &
1226 : n_pade, nspin
1227 : LOGICAL :: do_moments_ft, do_polarizability
1228 : REAL(kind=dp) :: damping, t0
1229 30 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: omegas, omegas_pade_real
1230 : REAL(kind=dp), DIMENSION(3) :: delta_vec
1231 : TYPE(cp_logger_type), POINTER :: logger
1232 : TYPE(section_vals_type), POINTER :: moment_ft_section, pol_section
1233 :
1234 : ! For results, using spin * direction for first index, e.g. for nspin = 2
1235 : ! results(1,:) = (spin=1 and direction=1,:),
1236 : ! results(5,:) = (spin=2 and direction=2,:)
1237 :
1238 30 : logger => cp_get_default_logger()
1239 :
1240 30 : moment_ft_section => section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS_FT")
1241 30 : pol_section => section_vals_get_subs_vals(rtp_section, "PRINT%POLARIZABILITY")
1242 :
1243 30 : nspin = SIZE(moments, 1)
1244 30 : n = SIZE(times)
1245 30 : n_elems = SIZE(rtc%print_pol_elements, 1)
1246 :
1247 30 : info_unit = -1
1248 30 : IF (PRESENT(info_opt)) info_unit = info_opt
1249 :
1250 : ! NOTE : Allows for at most 2 spin species
1251 30 : file_extensions(1) = "_SPIN_A.dat"
1252 30 : file_extensions(2) = "_SPIN_B.dat"
1253 :
1254 : ! Determine whether MOMENTS_FT and/or polarizability needs to be calculated
1255 30 : do_moments_ft = cp_printkey_is_on(logger%iter_info, moment_ft_section)
1256 30 : do_polarizability = cp_printkey_is_on(logger%iter_info, pol_section)
1257 30 : do_polarizability = do_polarizability .AND. (n_elems > 0)
1258 :
1259 30 : damping = rtc%ft_damping
1260 30 : t0 = rtc%ft_t0
1261 :
1262 : ! Determine field ft if polarizability required
1263 30 : IF (do_polarizability) THEN
1264 30 : ALLOCATE (field_results(3, n))
1265 10 : IF (rtc%apply_delta_pulse) THEN
1266 : ! Constant real FT
1267 8 : IF (PRESENT(cell)) THEN
1268 : delta_vec(:) = (REAL(rtc%delta_pulse_direction(1), kind=dp)*cell%h_inv(1, :) + &
1269 : REAL(rtc%delta_pulse_direction(2), kind=dp)*cell%h_inv(2, :) + &
1270 : REAL(rtc%delta_pulse_direction(3), kind=dp)*cell%h_inv(3, :)) &
1271 0 : *twopi*rtc%delta_pulse_scale
1272 : ELSE
1273 32 : delta_vec(:) = REAL(rtc%delta_pulse_direction(:), kind=dp)*rtc%delta_pulse_scale
1274 : END IF
1275 32 : DO k = 1, 3
1276 4256 : field_results(k, :) = CMPLX(delta_vec(k), 0.0, kind=dp)
1277 : END DO
1278 : ELSE
1279 : ! Do explicit FT
1280 : CALL multi_fft(times, fields, field_results, &
1281 2 : damping_opt=damping, t0_opt=t0, subtract_initial_opt=.TRUE.)
1282 : END IF
1283 : END IF
1284 :
1285 30 : IF (do_moments_ft .OR. do_polarizability) THEN
1286 : ! We need to transform at least the moments
1287 : ! NOTE : Might be able to save some memory by only doing FT of actually
1288 : ! required moments, but for now, doing FT of all moment directions
1289 40 : ALLOCATE (results(3*nspin, n))
1290 30 : ALLOCATE (omegas(n))
1291 30 : ALLOCATE (value_series(3*nspin, n))
1292 20 : DO i = 1, nspin
1293 50 : DO k = 1, 3
1294 4570 : value_series(3*(i - 1) + k, :) = moments(i, k, :)
1295 : END DO
1296 : END DO
1297 : ! TODO : Choose whether the initial subtraction is applied in &FT section?
1298 : CALL multi_fft(times, value_series, results, omegas, &
1299 10 : damping_opt=damping, t0_opt=t0, subtract_initial_opt=.TRUE.)
1300 10 : DEALLOCATE (value_series)
1301 20 : DO i = 1, nspin
1302 : ! Output to FT file, if needed
1303 : ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension=file_extensions(i), &
1304 10 : file_form="FORMATTED", file_position="REWIND")
1305 : ! Print header
1306 10 : IF (ft_unit > 0) THEN
1307 5 : ALLOCATE (headers(7))
1308 5 : headers(2) = " x,real [at.u.]"
1309 5 : headers(3) = " x,imag [at.u.]"
1310 5 : headers(4) = " y,real [at.u.]"
1311 5 : headers(5) = " y,imag [at.u.]"
1312 5 : headers(6) = " z,real [at.u.]"
1313 5 : headers(7) = " z,imag [at.u.]"
1314 5 : IF (info_unit == ft_unit) THEN
1315 0 : headers(1) = "# Energy [eV]"
1316 0 : prefix = " MOMENTS_FT|"
1317 0 : prefix_format = "(A12)"
1318 : CALL print_ft_file(ft_unit, headers, omegas, results(3*(i - 1) + 1:3*(i - 1) + 3, :), &
1319 0 : prefix, prefix_format, evolt)
1320 : ELSE
1321 5 : headers(1) = "# omega [at.u.]"
1322 5 : CALL print_ft_file(ft_unit, headers, omegas, results(3*(i - 1) + 1:3*(i - 1) + 3, :))
1323 : END IF
1324 5 : DEALLOCATE (headers)
1325 : END IF
1326 20 : CALL cp_print_key_finished_output(ft_unit, logger, moment_ft_section)
1327 : END DO
1328 : END IF
1329 :
1330 30 : IF (rtc%pade_requested .AND. (do_moments_ft .OR. do_polarizability)) THEN
1331 6 : ALLOCATE (omegas_complex(SIZE(omegas)))
1332 1004 : omegas_complex(:) = CMPLX(omegas(:), 0.0, kind=dp)
1333 2 : n_pade = INT((rtc%pade_e_max - rtc%pade_e_min)/rtc%pade_e_step)
1334 6 : ALLOCATE (omegas_pade(n_pade))
1335 6 : ALLOCATE (omegas_pade_real(n_pade))
1336 : ! Construct omegas_pade and omegas_complex
1337 2000 : DO i = 1, n_pade
1338 1998 : omegas_pade_real(i) = (i - 1)*rtc%pade_e_step + rtc%pade_e_min
1339 2000 : omegas_pade(i) = CMPLX(omegas_pade_real(i), 0.0, kind=dp)
1340 : END DO
1341 8000 : ALLOCATE (results_pade(nspin*3, n_pade), source=CMPLX(0.0, 0.0, kind=dp))
1342 4 : DO i = 1, nspin
1343 8 : DO k = 1, 3
1344 : CALL greenx_refine_ft(rtc%pade_fit_e_min, rtc%pade_fit_e_max, omegas_complex, results(3*(i - 1) + k, :), &
1345 8 : omegas_pade, results_pade(3*(i - 1) + k, :))
1346 : END DO
1347 : ! Print to a file
1348 : ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension="_PADE"//file_extensions(i), &
1349 2 : file_form="FORMATTED", file_position="REWIND")
1350 4 : IF (ft_unit > 0) THEN
1351 1 : ALLOCATE (headers(7))
1352 1 : headers(2) = " x,real,pade [at.u.]"
1353 1 : headers(3) = " x,imag,pade [at.u.]"
1354 1 : headers(4) = " y,real,pade [at.u.]"
1355 1 : headers(5) = " y,imag,pade [at.u.]"
1356 1 : headers(6) = " z,real,pade [at.u.]"
1357 1 : headers(7) = " z,imag,pade [at.u.]"
1358 1 : IF (info_unit == ft_unit) THEN
1359 0 : headers(1) = "# Energy [eV]"
1360 0 : prefix = " MOMENTS_FT_PADE|"
1361 0 : prefix_format = "(A17)"
1362 : CALL print_ft_file(ft_unit, headers, omegas_pade_real, results_pade(3*(i - 1) + 1:3*(i - 1) + 3, :), &
1363 0 : prefix, prefix_format, evolt)
1364 : ELSE
1365 1 : headers(1) = "# omega [at.u.]"
1366 1 : CALL print_ft_file(ft_unit, headers, omegas_pade_real, results_pade(3*(i - 1) + 1:3*(i - 1) + 3, :))
1367 : END IF
1368 1 : DEALLOCATE (headers)
1369 : END IF
1370 : END DO
1371 : END IF
1372 :
1373 30 : IF (do_polarizability) THEN
1374 : ! get the polarizability elements, as required
1375 40 : ALLOCATE (pol_results(n_elems, n))
1376 20 : DO i = 1, nspin
1377 40 : DO k = 1, n_elems
1378 : ! NOTE - field is regularized to small value
1379 : pol_results(k, :) = results(3*(i - 1) + &
1380 : rtc%print_pol_elements(k, 1), :)/ &
1381 : (field_results(rtc%print_pol_elements(k, 2), :) + &
1382 4570 : 1.0e-10*field_results(rtc%print_pol_elements(k, 2), 2))
1383 : END DO
1384 : ! Print to the file
1385 : ft_unit = cp_print_key_unit_nr(logger, pol_section, extension=file_extensions(i), &
1386 10 : file_form="FORMATTED", file_position="REWIND")
1387 10 : IF (ft_unit > 0) THEN
1388 15 : ALLOCATE (headers(2*n_elems + 1))
1389 20 : DO k = 1, n_elems
1390 15 : WRITE (headers(2*k), "(A16,I2,I2)") "real pol. elem.", &
1391 15 : rtc%print_pol_elements(k, 1), &
1392 30 : rtc%print_pol_elements(k, 2)
1393 15 : WRITE (headers(2*k + 1), "(A16,I2,I2)") "imag pol. elem.", &
1394 15 : rtc%print_pol_elements(k, 1), &
1395 35 : rtc%print_pol_elements(k, 2)
1396 : END DO
1397 : ! Write header
1398 5 : IF (info_unit == ft_unit) THEN
1399 2 : headers(1) = "# Energy [eV]"
1400 2 : prefix = " POLARIZABILITY|"
1401 2 : prefix_format = "(A16)"
1402 : CALL print_ft_file(ft_unit, headers, omegas, pol_results, &
1403 2 : prefix, prefix_format, evolt)
1404 : ELSE
1405 3 : headers(1) = "# omega [at.u.]"
1406 3 : CALL print_ft_file(ft_unit, headers, omegas, pol_results)
1407 : END IF
1408 5 : DEALLOCATE (headers)
1409 : END IF
1410 20 : CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1411 : END DO
1412 : END IF
1413 :
1414 : ! Padé polarizability
1415 30 : IF (rtc%pade_requested .AND. do_polarizability) THEN
1416 : ! Start with the field pade
1417 6 : ALLOCATE (field_results_pade(3, n_pade))
1418 2 : IF (rtc%apply_delta_pulse) THEN
1419 8 : DO k = 1, 3
1420 6002 : field_results_pade(k, :) = CMPLX(delta_vec(k), 0.0, kind=dp)
1421 : END DO
1422 : ELSE
1423 0 : DO k = 1, 3
1424 : CALL greenx_refine_ft(rtc%pade_fit_e_min, rtc%pade_fit_e_max, &
1425 : omegas_complex, field_results(k, :), &
1426 0 : omegas_pade, field_results_pade(k, :))
1427 : END DO
1428 : END IF
1429 : ! Allocate polarisation pade
1430 8 : ALLOCATE (pol_results_pade(n_elems, n_pade))
1431 : ! Refine
1432 4 : DO i = 1, nspin
1433 8 : DO k = 1, n_elems
1434 : ! NOTE : Regularization to small value
1435 : pol_results_pade(k, :) = results_pade(3*(i - 1) + rtc%print_pol_elements(k, 1), :)/( &
1436 : field_results_pade(rtc%print_pol_elements(k, 2), :) + &
1437 6002 : field_results_pade(rtc%print_pol_elements(k, 2), 2)*1.0e-10_dp)
1438 : END DO
1439 : ! Print to the file
1440 : ft_unit = cp_print_key_unit_nr(logger, pol_section, extension="_PADE"//file_extensions(i), &
1441 2 : file_form="FORMATTED", file_position="REWIND")
1442 2 : IF (ft_unit > 0) THEN
1443 3 : ALLOCATE (headers(2*n_elems + 1))
1444 4 : DO k = 1, n_elems
1445 3 : WRITE (headers(2*k), "(A16,I2,I2)") "re,pade,pol.", &
1446 3 : rtc%print_pol_elements(k, 1), &
1447 6 : rtc%print_pol_elements(k, 2)
1448 3 : WRITE (headers(2*k + 1), "(A16,I2,I2)") "im,pade,pol.", &
1449 3 : rtc%print_pol_elements(k, 1), &
1450 7 : rtc%print_pol_elements(k, 2)
1451 : END DO
1452 : ! Write header
1453 1 : IF (info_unit == ft_unit) THEN
1454 1 : headers(1) = "# Energy [eV]"
1455 1 : prefix = " POLARIZABILITY_PADE|"
1456 1 : prefix_format = "(A21)"
1457 : CALL print_ft_file(ft_unit, headers, omegas_pade_real, pol_results_pade, &
1458 1 : prefix, prefix_format, evolt)
1459 : ELSE
1460 0 : headers(1) = "# omega [at.u.]"
1461 0 : CALL print_ft_file(ft_unit, headers, omegas_pade_real, pol_results_pade)
1462 : END IF
1463 1 : DEALLOCATE (headers)
1464 : END IF
1465 4 : CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1466 : END DO
1467 2 : DEALLOCATE (field_results_pade)
1468 2 : DEALLOCATE (pol_results_pade)
1469 : END IF
1470 :
1471 30 : IF (rtc%pade_requested .AND. (do_moments_ft .OR. do_polarizability)) THEN
1472 2 : DEALLOCATE (omegas_complex)
1473 2 : DEALLOCATE (omegas_pade)
1474 2 : DEALLOCATE (omegas_pade_real)
1475 2 : DEALLOCATE (results_pade)
1476 : END IF
1477 :
1478 30 : IF (do_polarizability) THEN
1479 10 : DEALLOCATE (pol_results)
1480 10 : DEALLOCATE (field_results)
1481 : END IF
1482 :
1483 30 : IF (do_moments_ft .OR. do_polarizability) THEN
1484 10 : DEALLOCATE (results)
1485 10 : DEALLOCATE (omegas)
1486 : END IF
1487 30 : END SUBROUTINE print_ft
1488 :
1489 : ! **************************************************************************************************
1490 : !> \brief ...
1491 : !> \param ft_unit ...
1492 : !> \param headers ...
1493 : !> \param xvals ...
1494 : !> \param yvals ...
1495 : !> \param prefix ...
1496 : !> \param prefix_format ...
1497 : !> \param xscale_opt ...
1498 : ! **************************************************************************************************
1499 12 : SUBROUTINE print_ft_file(ft_unit, headers, xvals, yvals, prefix, prefix_format, xscale_opt)
1500 : INTEGER, INTENT(IN) :: ft_unit
1501 : CHARACTER(len=20), DIMENSION(:), INTENT(IN) :: headers
1502 : REAL(kind=dp), DIMENSION(:), INTENT(IN) :: xvals
1503 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: yvals
1504 : CHARACTER(len=21), INTENT(IN), OPTIONAL :: prefix
1505 : CHARACTER(len=5), INTENT(IN), OPTIONAL :: prefix_format
1506 : REAL(kind=dp), INTENT(IN), OPTIONAL :: xscale_opt
1507 :
1508 : INTEGER :: i, j, ncols, nrows
1509 : LOGICAL :: do_prefix
1510 : REAL(kind=dp) :: xscale
1511 :
1512 12 : do_prefix = .FALSE.
1513 12 : IF (PRESENT(prefix)) THEN
1514 3 : IF (PRESENT(prefix_format)) THEN
1515 : do_prefix = .TRUE.
1516 : ELSE
1517 0 : CPABORT("Printing of prefix with missing format!")
1518 : END IF
1519 : END IF
1520 :
1521 12 : xscale = 1.0_dp
1522 12 : IF (PRESENT(xscale_opt)) xscale = xscale_opt
1523 :
1524 12 : ncols = SIZE(yvals, 1)
1525 12 : nrows = SIZE(yvals, 2)
1526 :
1527 : ! Check whether enough headers for yvals and xvals is present
1528 12 : IF (SIZE(headers) < 2*ncols + 1) THEN
1529 0 : CPABORT("Not enough headers to print the file!")
1530 : END IF
1531 :
1532 12 : IF (SIZE(xvals) < nrows) THEN
1533 0 : CPABORT("Not enough xvals to print all yvals!")
1534 : END IF
1535 :
1536 12 : IF (ft_unit > 0) THEN
1537 : ! If prefix is present, write prefix
1538 : ! Also write energy in eV instead of at.u.
1539 12 : IF (do_prefix) THEN
1540 3 : WRITE (ft_unit, prefix_format, advance="no") prefix
1541 : END IF
1542 12 : WRITE (ft_unit, "(A20)", advance="no") headers(1)
1543 : ! Print the rest of the headers
1544 72 : DO j = 1, 2*ncols - 1
1545 72 : WRITE (ft_unit, "(A20)", advance="no") headers(j + 1)
1546 : END DO
1547 12 : WRITE (ft_unit, "(A20)") headers(2*ncols + 1)
1548 : ! Done with the headers, print actual data
1549 3520 : DO i = 1, nrows
1550 3508 : IF (do_prefix) THEN
1551 1551 : WRITE (ft_unit, prefix_format, advance="no") prefix
1552 : END IF
1553 3508 : WRITE (ft_unit, "(E20.8E3)", advance="no") xvals(i)*xscale
1554 10524 : DO j = 1, ncols - 1
1555 : WRITE (ft_unit, "(E20.8E3,E20.8E3)", advance="no") &
1556 10524 : REAL(yvals(j, i)), AIMAG(yvals(j, i))
1557 : END DO
1558 : WRITE (ft_unit, "(E20.8E3,E20.8E3)") &
1559 3520 : REAL(yvals(ncols, i)), AIMAG(yvals(ncols, i))
1560 : END DO
1561 : END IF
1562 12 : END SUBROUTINE print_ft_file
1563 :
1564 : END MODULE rt_propagation_output
|
