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