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