(git:d3d49ac)
Loading...
Searching...
No Matches
rt_propagation_output.F
Go to the documentation of this file.
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
13MODULE 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
125CONTAINS
126
127! **************************************************************************************************
128!> \brief ...
129!> \param qs_env ...
130!> \param run_type ...
131!> \param delta_iter ...
132!> \param used_time ...
133! **************************************************************************************************
134 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
142 INTEGER, DIMENSION(2) :: nelectron_spin
143 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 REAL(kind=dp), DIMENSION(:), POINTER :: qs_tot_rho_r
148 REAL(kind=dp), DIMENSION(:, :), POINTER :: j_int
149 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
150 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
151 TYPE(cp_logger_type), POINTER :: logger
152 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
153 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 POINTER :: sab_all, sab_orb
158 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
159 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 NULLIFY (logger, dft_control)
165
166 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 dbcsr_dist=dbcsr_dist, nelectron_spin=nelectron_spin)
177
178 rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
179
180 CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
181 n_electrons = n_electrons - dft_control%charge
182
183 CALL qs_rho_get(rho_struct=rho, tot_rho_r=qs_tot_rho_r)
184
185 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 extension=".scfLog")
189
190 IF (output_unit > 0) THEN
191 WRITE (output_unit, fmt="(/,(T3,A,T40,I5))") &
192 "Information at iteration step:", rtp%iter
193 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
194 "Total electronic density (r-space): ", &
195 tot_rho_r, &
196 tot_rho_r + &
197 REAL(n_electrons, dp)
198 WRITE (unit=output_unit, fmt="((T3,A,T59,F22.14))") &
199 "Total energy:", rtp%energy_new
200 IF (run_type == ehrenfest) THEN
201 WRITE (unit=output_unit, fmt="((T3,A,T61,F20.14))") &
202 "Energy difference to previous iteration step:", rtp%energy_new - rtp%energy_old
203 END IF
204 IF (run_type == real_time_propagation) THEN
205 WRITE (unit=output_unit, fmt="((T3,A,T61,F20.14))") &
206 "Energy difference to initial state:", rtp%energy_new - rtp%energy_old
207 END IF
208 IF (PRESENT(delta_iter)) THEN
209 WRITE (unit=output_unit, fmt="((T3,A,T61,E20.6))") &
210 "Convergence:", delta_iter
211 END IF
212 IF (rtp%converged) THEN
213 IF (run_type == real_time_propagation) THEN
214 WRITE (unit=output_unit, fmt="((T3,A,T61,F12.2))") &
215 "Time needed for propagation:", used_time
216 END IF
217 WRITE (unit=output_unit, fmt="(/,(T3,A,3X,F16.14))") &
218 "CONVERGENCE REACHED", rtp%energy_new - rtp%energy_old
219 END IF
220 END IF
221
222 IF (rtp%converged) THEN
223 IF (.NOT. rtp%linear_scaling) THEN
224 CALL get_rtp(rtp=rtp, mos_new=mos_new)
225 CALL rt_calculate_orthonormality(orthonormality, &
226 mos_new, matrix_s(1)%matrix)
227 IF (output_unit > 0) THEN
228 WRITE (output_unit, fmt="(/,(T3,A,T60,F20.10))") &
229 "Max deviation from orthonormalization:", orthonormality
230 END IF
231 END IF
232 END IF
233
234 IF (output_unit > 0) THEN
235 CALL m_flush(output_unit)
236 END IF
237 CALL cp_print_key_finished_output(output_unit, logger, rtp_section, &
238 "PRINT%PROGRAM_RUN_INFO")
239
240 IF (rtp%converged) THEN
241 dft_section => section_vals_get_subs_vals(input, "DFT")
242 IF (btest(cp_print_key_should_output(logger%iter_info, &
243 dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) THEN
244 CALL print_field_applied(qs_env, dft_section)
245 END IF
246 CALL make_moment(qs_env)
247 IF (btest(cp_print_key_should_output(logger%iter_info, &
248 dft_section, "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS"), cp_p_file)) THEN
249 CALL print_rtp_energy_components(qs_env, dft_section)
250 END IF
251 IF (.NOT. dft_control%qs_control%dftb) THEN
252 CALL write_available_results(qs_env=qs_env, rtp=rtp)
253 END IF
254
255 IF (rtp%linear_scaling) THEN
256 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 IF (dft_control%rtp_control%save_local_moments) THEN
260 ! Save the field value
261 IF (dft_control%apply_efield_field) THEN
262 CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
263 rtp%fields(:, rtp%istep + rtp%i_start + 1) = cmplx(field(:), 0.0, kind=dp)
264 END IF
265 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
266 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 rtp%local_moments_work, rtp%moments(:, :, rtp%istep + rtp%i_start + 1))
273 ! Time 1 is zero (start) time
274 rtp%times(rtp%istep + rtp%i_start + 1) = qs_env%sim_time
275 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 rtp%moments(:, :, rtp%istep + rtp%i_start + 1), qs_env%sim_time, rtp%track_imag_density)
278 END IF
279
280 IF (btest(cp_print_key_should_output(logger%iter_info, &
281 dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART"), cp_p_file)) THEN
282 CALL write_rt_p_to_restart(rho_new, .false.)
283 END IF
284 IF (btest(cp_print_key_should_output(logger%iter_info, &
285 dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY"), cp_p_file)) THEN
286 CALL write_rt_p_to_restart(rho_new, .true.)
287 END IF
288 IF (.NOT. dft_control%qs_control%dftb) THEN
289 !Not sure if these things could also work with dftb or not
290 IF (btest(cp_print_key_should_output(logger%iter_info, &
291 dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
292 DO spin = 1, SIZE(rho_new)/2
293 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 CALL get_rtp(rtp=rtp, mos_new=mos_new)
299 IF (.NOT. dft_control%qs_control%dftb .AND. .NOT. dft_control%qs_control%xtb) THEN
300 IF (rtp%track_imag_density) THEN
301 NULLIFY (p_im, p_xyz)
302 CALL dbcsr_allocate_matrix_set(p_xyz, 3)
303
304! Linear momentum operator
305! prepare for allocation
306 natom = SIZE(particle_set, 1)
307 ALLOCATE (first_sgf(natom))
308 ALLOCATE (last_sgf(natom))
309 CALL get_particle_set(particle_set, qs_kind_set, &
310 first_sgf=first_sgf, &
311 last_sgf=last_sgf)
312 ALLOCATE (row_blk_sizes(natom))
313 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
314 DEALLOCATE (first_sgf)
315 DEALLOCATE (last_sgf)
316
317 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 mutable_work=.true.)
323 CALL cp_dbcsr_alloc_block_from_nbl(p_xyz(1)%matrix, sab_orb)
324 CALL dbcsr_set(p_xyz(1)%matrix, 0.0_dp)
325 DO i = 2, 3
326 ALLOCATE (p_xyz(i)%matrix)
327 CALL dbcsr_copy(p_xyz(i)%matrix, p_xyz(1)%matrix, "p_xyz-"//trim(adjustl(cp_to_string(i))))
328 CALL dbcsr_set(p_xyz(i)%matrix, 0.0_dp)
329 END DO
330 CALL build_lin_mom_matrix(qs_env, p_xyz)
331 DEALLOCATE (row_blk_sizes)
332
333 nspin = SIZE(mos_new)/2
334 CALL qs_rho_get(rho, rho_ao_im=p_im)
335 ALLOCATE (j_int(nspin, 3))
336 j_int = 0.0_dp
337
338 NULLIFY (tmp_ao)
339 CALL dbcsr_init_p(tmp_ao)
340 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
341 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
342 CALL dbcsr_set(tmp_ao, 0.0_dp)
343
344 DO i = 1, 3
345 strace = 0.0_dp
346 DO spin = 1, nspin
347 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 0.0_dp, tmp_ao)
350 CALL dbcsr_trace(tmp_ao, trace)
351 strace = strace + trace
352! dft_control%rtp_control%vec_pot(1)
353 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 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 output_unit = cp_logger_get_default_io_unit(logger)
363 unit_nr = cp_print_key_unit_nr(logger, dft_section, &
364 "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT", extension=".dat", is_new_file=new_file)
365
366 IF (output_unit > 0) THEN
367 IF (new_file) THEN
368 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 "Step Nr.", "Time[fs]", "ALPHA jint[X]", "jint[Y]", "jint[Z]", &
371 "BETA jint[X]", "jint[Y]", "jint[Z]"
372 ELSE
373 WRITE (unit=unit_nr, fmt='("#",5X,A,4X,A,8X,A,2(10X,A))') "Step Nr.", "Time[fs]", &
374 "jint[X]", "jint[Y]", "jint[Z]"
375 END IF
376 END IF
377
378 IF (nspin == 2) THEN
379 WRITE (unit=unit_nr, fmt="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
380 j_int(1, 1:3), j_int(2, 1:3)
381 ELSE
382 WRITE (unit=unit_nr, fmt="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
383 j_int(1, 1:3)
384 END IF
385 END IF
386 CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
387 "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT")
388 END IF
389 DEALLOCATE (j_int)
390
391 IF (btest(cp_print_key_should_output(logger%iter_info, &
392 dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
393 DO spin = 1, nspin
394 CALL rt_current(qs_env, p_im(spin)%matrix, dft_section, spin, nspin)
395 END DO
396 END IF
397 CALL dbcsr_deallocate_matrix(tmp_ao)
398 CALL dbcsr_deallocate_matrix_set(p_xyz)
399 END IF
400
401! projection of molecular orbitals
402 IF (dft_control%rtp_control%is_proj_mo) THEN
403 DO n_proj = 1, SIZE(dft_control%rtp_control%proj_mo_list)
404 CALL compute_and_write_proj_mo(qs_env, mos_new, &
405 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 dft_section, qs_kind_set)
411 END IF
412 END IF
413
414 rtp%energy_old = rtp%energy_new
415
416 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 "or use a smaller TIMESTEP")
419 END IF
420
421 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 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 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 NULLIFY (tmp_fm_struct)
446
447 CALL timeset(routinen, handle)
448
449 nspin = SIZE(mos_new)/2
450 max_alpha = 0.0_dp
451 max_beta = 0.0_dp
452 DO ispin = 1, nspin
453 re = ispin*2 - 1
454 im = ispin*2
455 ! get S*C
456 CALL cp_fm_create(svec_re, mos_new(im)%matrix_struct)
457 CALL cp_fm_create(svec_im, mos_new(im)%matrix_struct)
458 CALL cp_fm_get_info(mos_new(im), &
459 nrow_global=n, ncol_global=k)
460 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(re), &
461 svec_re, k)
462 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(im), &
463 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 context=mos_new(re)%matrix_struct%context)
469 CALL cp_fm_create(overlap_re, tmp_fm_struct)
470
471 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 svec_re, 0.0_dp, overlap_re)
475 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(im), &
476 svec_im, 1.0_dp, overlap_re)
477
478 CALL cp_fm_release(svec_re)
479 CALL cp_fm_release(svec_im)
480
481 CALL cp_fm_get_info(overlap_re, nrow_local=nrow_local, ncol_local=ncol_local, &
482 row_indices=row_indices, col_indices=col_indices)
483 DO i = 1, nrow_local
484 DO j = 1, ncol_local
485 alpha = overlap_re%local_data(i, j)
486 IF (row_indices(i) == col_indices(j)) alpha = alpha - 1.0_dp
487 max_alpha = max(max_alpha, abs(alpha))
488 END DO
489 END DO
490 CALL cp_fm_release(overlap_re)
491 END DO
492 CALL mos_new(1)%matrix_struct%para_env%max(max_alpha)
493 CALL mos_new(1)%matrix_struct%para_env%max(max_beta)
494 orthonormality = max_alpha
495
496 CALL timestop(handle)
497
498 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 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 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
525
526 NULLIFY (tmp_fm_struct)
527
528 CALL timeset(routinen, handle)
529
530 CALL get_rtp(rtp=rtp, mos_new=mos_new)
531
532 nspin = SIZE(delta_mos)/2
533 max_alpha = 0.0_dp
534
535 DO i = 1, SIZE(mos_new)
536 CALL cp_fm_scale_and_add(-one, delta_mos(i), one, mos_new(i))
537 END DO
538
539 DO ispin = 1, nspin
540 re = ispin*2 - 1
541 im = ispin*2
542
543 double_col = .true.
544 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 double_row)
550
551 CALL cp_fm_create(work, matrix_struct=newstruct)
552 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 nrow_global=nao)
556 CALL cp_fm_get_info(work, ncol_global=newdim)
557
558 CALL cp_fm_set_all(work, zero, zero)
559
560 DO icol = 1, lcol
561 work%local_data(:, icol) = delta_mos(re)%local_data(:, icol)
562 work%local_data(:, icol + lcol) = delta_mos(im)%local_data(:, icol)
563 END DO
564
565 CALL cp_dbcsr_sm_fm_multiply(matrix_s, work, work1, ncol=newdim)
566
567 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 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 double_row)
577
578 CALL cp_fm_create(work, matrix_struct=newstruct1)
579 CALL cp_fm_create(work2, matrix_struct=newstruct1)
580
581 CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(re), &
582 work1, zero, work)
583
584 CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(im), &
585 work1, zero, work2)
586
587 CALL cp_fm_get_info(work, nrow_local=lrow)
588 DO i = 1, lrow
589 DO j = 1, lcol
590 alpha = sqrt((work%local_data(i, j) + work2%local_data(i, j + lcol))**2 + &
591 (work%local_data(i, j + lcol) - work2%local_data(i, j))**2)
592 max_alpha = max(max_alpha, abs(alpha))
593 END DO
594 END DO
595
596 CALL cp_fm_release(work)
597 CALL cp_fm_release(work1)
598 CALL cp_fm_release(work2)
599 CALL cp_fm_struct_release(tmp_fm_struct)
600 CALL cp_fm_struct_release(newstruct)
601 CALL cp_fm_struct_release(newstruct1)
602
603 END DO
604
605 CALL delta_mos(1)%matrix_struct%para_env%max(max_alpha)
606 delta_eps = sqrt(max_alpha)
607
608 CALL timestop(handle)
609
610 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 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 REAL(dp), DIMENSION(:, :), POINTER :: block_values
632 TYPE(dbcsr_iterator_type) :: iter
633 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
634 TYPE(dbcsr_type), POINTER :: tmp
635 TYPE(mp_comm_type) :: group
636
637 CALL timeset(routinen, handle)
638
639 CALL get_rtp(rtp=rtp, rho_new=rho_new)
640
641 DO i = 1, SIZE(rho_new)
642 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 DO i = 1, SIZE(delta_p)
646 !square all entries of both matrices
647 CALL dbcsr_iterator_start(iter, delta_p(i)%matrix)
648 DO WHILE (dbcsr_iterator_blocks_left(iter))
649 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
650 block_values = block_values*block_values
651 END DO
652 CALL dbcsr_iterator_stop(iter)
653 END DO
654 NULLIFY (tmp)
655 ALLOCATE (tmp)
656 CALL dbcsr_create(tmp, template=delta_p(1)%matrix, matrix_type="N")
657 DO ispin = 1, SIZE(delta_p)/2
658 CALL dbcsr_desymmetrize(delta_p(2*ispin - 1)%matrix, tmp)
659 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 max_alpha = zero
663 DO ispin = 1, SIZE(delta_p)/2
664 CALL dbcsr_iterator_start(iter, delta_p(2*ispin)%matrix)
665 DO WHILE (dbcsr_iterator_blocks_left(iter))
666 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
667 alpha = maxval(block_values)
668 IF (alpha > max_alpha) max_alpha = alpha
669 END DO
670 CALL dbcsr_iterator_stop(iter)
671 END DO
672 CALL dbcsr_get_info(delta_p(1)%matrix, group=group)
673 CALL group%max(max_alpha)
674 delta_eps = sqrt(max_alpha)
675 CALL dbcsr_deallocate_matrix(tmp)
676 CALL timestop(handle)
677
678 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 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 CALL timeset(routinen, handle)
697
698 NULLIFY (dft_control)
699
700 logger => cp_get_default_logger()
701 output_unit = cp_logger_get_default_io_unit(logger)
702 CALL get_qs_env(qs_env, dft_control=dft_control)
703 IF (dft_control%qs_control%dftb) THEN
704 CALL scf_post_calculation_tb(qs_env, "DFTB", .false.)
705 ELSE IF (dft_control%qs_control%xtb) THEN
706 CALL scf_post_calculation_tb(qs_env, "xTB", .false.)
707 ELSE
708 CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
709 END IF
710 CALL timestop(handle)
711
712 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmp
732
733 CALL timeset(routinen, handle)
734
735 logger => cp_get_default_logger()
736 unit_nr = cp_logger_get_default_io_unit(logger)
737 NULLIFY (tmp)
738 CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
739 DO i = 1, SIZE(rho)
740 CALL dbcsr_init_p(tmp(i)%matrix)
741 CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
742 CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
743 END DO
744 DO ispin = 1, SIZE(rho)/2
745 re = 2*ispin - 1
746 im = 2*ispin
747 eps = max(filter_eps, 1.0e-11_dp)
748 DO WHILE (eps < 1.1_dp)
749 CALL dbcsr_filter(tmp(re)%matrix, eps)
750 occ = dbcsr_get_occupation(tmp(re)%matrix)
751 IF (unit_nr > 0) WRITE (unit_nr, fmt="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
752 ispin, " eps ", eps, " real: ", occ
753 eps = eps*10
754 END DO
755 eps = max(filter_eps, 1.0e-11_dp)
756 DO WHILE (eps < 1.1_dp)
757 CALL dbcsr_filter(tmp(im)%matrix, eps)
758 occ = dbcsr_get_occupation(tmp(im)%matrix)
759 IF (unit_nr > 0) WRITE (unit_nr, fmt="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
760 ispin, " eps ", eps, " imag: ", occ
761 eps = eps*10.0_dp
762 END DO
763 END DO
764 CALL dbcsr_deallocate_matrix_set(tmp)
765 CALL timestop(handle)
766
767 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 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 CALL timeset(routinen, handle)
789 logger => cp_get_default_logger()
790 IF (logger%para_env%is_source()) THEN
791 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
792 ELSE
793 unit_nr = -1
794 END IF
795
796 project_name = logger%iter_info%project_name
797 DO ispin = 1, SIZE(rho_new)/2
798 re = 2*ispin - 1
799 im = 2*ispin
800 IF (history) THEN
801 WRITE (file_name, '(A,I0,A)') &
802 trim(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//trim(cp_iter_string(logger%iter_info))//"_RESTART.dm"
803 ELSE
804 WRITE (file_name, '(A,I0,A)') trim(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
805 END IF
806 cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.true.)
807 IF (unit_nr > 0) THEN
808 WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//trim(file_name)//" with checksum: ", cs_pos
809 END IF
810 CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
811 IF (history) THEN
812 WRITE (file_name, '(A,I0,A)') &
813 trim(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//trim(cp_iter_string(logger%iter_info))//"_RESTART.dm"
814 ELSE
815 WRITE (file_name, '(A,I0,A)') trim(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
816 END IF
817 cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.true.)
818 IF (unit_nr > 0) THEN
819 WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//trim(file_name)//" with checksum: ", cs_pos
820 END IF
821 CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
822 END DO
823
824 CALL timestop(handle)
825
826 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 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 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 CALL timeset(routinen, handle)
862
863 logger => cp_get_default_logger()
864 CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
865 CALL qs_subsys_get(subsys, particles=particles)
866 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
867
868 NULLIFY (zero, tmp)
869 ALLOCATE (zero, tmp)
870 CALL dbcsr_create(zero, template=p_im)
871 CALL dbcsr_copy(zero, p_im)
872 CALL dbcsr_set(zero, 0.0_dp)
873 CALL dbcsr_create(tmp, template=p_im)
874 CALL dbcsr_copy(tmp, p_im)
875 IF (nspin == 1) THEN
876 CALL dbcsr_scale(tmp, 0.5_dp)
877 END IF
878 current_env%gauge = -1
879 current_env%gauge_init = .false.
880 CALL auxbas_pw_pool%create_pw(rs)
881 CALL auxbas_pw_pool%create_pw(gs)
882
883 NULLIFY (stride)
884 ALLOCATE (stride(3))
885
886 DO dir = 1, 3
887
888 CALL pw_zero(rs)
889 CALL pw_zero(gs)
890
891 CALL calculate_jrho_resp(zero, tmp, zero, zero, dir, dir, rs, gs, qs_env, current_env, retain_rsgrid=.true.)
892
893 stride = section_get_ivals(dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT%STRIDE")
894
895 IF (dir == 1) THEN
896 sdir = "-x"
897 ELSE IF (dir == 2) THEN
898 sdir = "-y"
899 ELSE
900 sdir = "-z"
901 END IF
902 WRITE (char_spin, "(I1)") spin
903
904 ext = "-SPIN-"//char_spin//sdir//".cube"
905 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 log_filename=.false., mpi_io=mpi_io)
909
910 CALL cp_pw_to_cube(rs, print_unit, "EMD current", particles=particles, stride=stride, &
911 mpi_io=mpi_io)
912
913 CALL cp_print_key_finished_output(print_unit, logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
914 mpi_io=mpi_io)
915
916 END DO
917
918 CALL auxbas_pw_pool%give_back_pw(rs)
919 CALL auxbas_pw_pool%give_back_pw(gs)
920
921 CALL dbcsr_deallocate_matrix(zero)
922 CALL dbcsr_deallocate_matrix(tmp)
923
924 DEALLOCATE (stride)
925
926 CALL timestop(handle)
927
928 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 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 CALL timeset(routinen, handle)
950
951 CALL get_qs_env(qs_env, scf_env=scf_env)
952 IF (rtp%linear_scaling) THEN
953 CALL write_mo_free_results(qs_env)
954 ELSE
955 CALL write_mo_free_results(qs_env)
956 CALL write_mo_dependent_results(qs_env, scf_env)
957 ! Time-dependent MO print
958 CALL write_rtp_mos_to_output_unit(qs_env, rtp)
959 CALL write_rtp_mo_cubes(qs_env, rtp)
960 END IF
961
962 CALL timestop(handle)
963
964 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 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 NULLIFY (dft_control)
988
989 logger => cp_get_default_logger()
990 output_unit = cp_logger_get_default_io_unit(logger)
991
992 CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp)
993
994 i_step = rtp%istep
995
996 unit_nr = cp_print_key_unit_nr(logger, dft_section, &
997 "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat", is_new_file=new_file)
998
999 IF (output_unit > 0) THEN
1000 rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
1001 IF (unit_nr /= output_unit) THEN
1002 INQUIRE (unit=unit_nr, name=filename)
1003 WRITE (unit=output_unit, fmt="(/,T2,A,2(/,T3,A),/)") &
1004 "FIELD", "The field applied is written to the file:", &
1005 trim(filename)
1006 ELSE
1007 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 (trim(rlab(i)), "=", dft_control%rtp_control%field(i), i=1, 3)
1010 END IF
1011
1012 IF (new_file) THEN
1013 IF (dft_control%apply_efield_field) THEN
1014 WRITE (unit=unit_nr, fmt='("#",5X,A,8X,A,3(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z"
1015 ELSE IF (dft_control%apply_vector_potential) THEN
1016 WRITE (unit=unit_nr, fmt='("#",5X,A,8X,A,6(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z", &
1017 " Vec. Pot. X", " Vec. Pot. Y", " Vec. Pot. Z"
1018 END IF
1019 END IF
1020
1021 field = 0.0_dp
1022 IF (dft_control%apply_efield_field) THEN
1023 CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
1024 WRITE (unit=unit_nr, fmt="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1025 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 ELSE IF (dft_control%apply_vector_potential) THEN
1030 WRITE (unit=unit_nr, fmt="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1031 dft_control%rtp_control%field(1), dft_control%rtp_control%field(2), dft_control%rtp_control%field(3), &
1032 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 "REAL_TIME_PROPAGATION%PRINT%FIELD")
1039
1040 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 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 NULLIFY (dft_control, energy, rtp)
1062
1063 logger => cp_get_default_logger()
1064 output_unit = cp_logger_get_default_io_unit(logger)
1065
1066 CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
1067 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 file_action="WRITE", is_new_file=new_file)
1072
1073 IF (output_unit > 0) THEN
1074 IF (unit_nr /= output_unit) THEN
1075 INQUIRE (unit=unit_nr, name=filename)
1076 WRITE (unit=output_unit, fmt="(/,T2,A,2(/,T3,A),/)") &
1077 "ENERGY_CONSTITUENTS", "Total Energy constituents written to file:", &
1078 trim(filename)
1079 ELSE
1080 WRITE (unit=output_unit, fmt="(/,T2,A)") "ENERGY_CONSTITUENTS"
1081 END IF
1082
1083 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 WRITE (unit=unit_nr, fmt='("#",5X,A,8X,A,10(6X,A))') "Step Nr.", "Time[fs]", &
1087 "Total ener.[a.u.]", "core[a.u.] ", " overlap [a.u.]", "hartree[a.u.]", " exc. [a.u.] ", &
1088 " 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 qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1093 energy%total, energy%core, energy%core_overlap, energy%hartree, energy%exc, &
1094 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 "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS")
1100
1101 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 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 nspin = SIZE(moments, 1)
1131 ndir = SIZE(moments, 2)
1132
1133 IF (nspin < 1) cpabort("Zero spin index size in print moments!")
1134 IF (ndir < 1) cpabort("Zero direction index size in print moments!")
1135
1136 imaginary = .true.
1137 IF (PRESENT(imag_opt)) imaginary = imag_opt
1138
1139 append = .true.
1140 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 file_extensions(1) = "_SPIN_A_RE.dat"
1146 file_extensions(2) = "_SPIN_A_IM.dat"
1147 file_extensions(3) = "_SPIN_B_RE.dat"
1148 file_extensions(4) = "_SPIN_B_IM.dat"
1149 logger => cp_get_default_logger()
1150 DO i = 1, nspin
1151 moment_t(:, 1) = moments(i, :)
1152 ! Real part
1153 print_unit = cp_print_key_unit_nr(logger, moments_section, extension=file_extensions(2*i - 1))
1154 IF (print_unit == info_unit) THEN
1155 ! print with prefix
1156 prefix = " MOMENTS_TRACE_RE|"
1157 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 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 xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
1170 END IF
1171 ELSE
1172 ! Print without prefix
1173 IF (append) THEN
1174 ! Print without headers
1175 CALL print_rt_file(print_unit, xvals=[time], yvals=moment_t, &
1176 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 xscale_opt=femtoseconds, comp_opt=rt_file_comp_real)
1184 END IF
1185 END IF
1186 CALL cp_print_key_finished_output(print_unit, logger, moments_section)
1187 ! Same for imaginary part
1188 IF (imaginary) THEN
1189 print_unit = cp_print_key_unit_nr(logger, moments_section, extension=file_extensions(2*i))
1190 IF (print_unit == info_unit) THEN
1191 ! print with prefix
1192 prefix = " MOMENTS_TRACE_IM|"
1193 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 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 xscale_opt=femtoseconds, comp_opt=rt_file_comp_imag)
1206 END IF
1207 ELSE
1208 ! Print without prefix
1209 IF (append) THEN
1210 ! Print without headers
1211 CALL print_rt_file(print_unit, xvals=[time], yvals=moment_t, &
1212 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 xscale_opt=femtoseconds, comp_opt=rt_file_comp_imag)
1220 END IF
1221 END IF
1222 CALL cp_print_key_finished_output(print_unit, logger, moments_section)
1223 END IF
1224 END DO
1225
1226 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 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 imag = .false.
1249 IF (PRESENT(imag_opt)) imag = imag_opt
1250 nspin = SIZE(density_matrices)/2
1251
1252 DO i = 1, nspin
1253 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 0.0_dp, work)
1257 CALL dbcsr_trace(work, real_moment)
1258 moment(i, k) = cmplx(real_moment, 0.0, kind=dp)
1259 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_dp, work)
1263 CALL dbcsr_trace(work, real_moment)
1264 moment(i, k) = moment(i, k) + cmplx(0.0, real_moment, kind=dp)
1265 END IF
1266 END DO
1267 END DO
1268
1269 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 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=*), PARAMETER :: routinen = 'print_ft'
1293
1294 CHARACTER(len=11), DIMENSION(2) :: file_extensions
1295 CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: headers
1296 CHARACTER(len=21) :: prefix
1297 CHARACTER(len=5) :: prefix_format
1298 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: omegas_complex, omegas_pade
1299 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: field_results, field_results_pade, &
1300 pol_results, pol_results_pade, pol_results_pade_spin_total, pol_results_spin_total, &
1301 results, results_pade, results_pade_spin_total, results_spin_total, value_series
1302 INTEGER :: ft_unit, handle, i, idx_omega_zero, &
1303 info_unit, k, k_static, n, n_elems, &
1304 n_pade, nspin
1305 LOGICAL :: do_moments_ft, do_polarizability
1306 REAL(kind=dp) :: damping, t0
1307 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: omegas, omegas_pade_real
1308 REAL(kind=dp), DIMENSION(3) :: delta_vec
1309 TYPE(cp_logger_type), POINTER :: logger
1310 TYPE(section_vals_type), POINTER :: moment_ft_section, pol_section
1311
1312 CALL timeset(routinen, handle)
1313! For results, using spin * direction for first index, e.g. for nspin = 2
1314! results(1,:) = (spin=1 and direction=1,:),
1315! results(5,:) = (spin=2 and direction=2,:)
1316
1317 logger => cp_get_default_logger()
1318
1319 moment_ft_section => section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS_FT")
1320 pol_section => section_vals_get_subs_vals(rtp_section, "PRINT%POLARIZABILITY")
1321
1322 nspin = SIZE(moments, 1)
1323 n = SIZE(times)
1324 n_elems = SIZE(rtc%print_pol_elements, 1)
1325
1326 info_unit = -1
1327 IF (PRESENT(info_opt)) info_unit = info_opt
1328
1329 ! NOTE : Allows for at most 2 spin species
1330 file_extensions(1) = "_SPIN_A.dat"
1331 file_extensions(2) = "_SPIN_B.dat"
1332
1333 ! Determine whether MOMENTS_FT and/or polarizability needs to be calculated
1334 do_moments_ft = cp_printkey_is_on(logger%iter_info, moment_ft_section)
1335 do_polarizability = cp_printkey_is_on(logger%iter_info, pol_section)
1336 do_polarizability = do_polarizability .AND. (n_elems > 0)
1337
1338 damping = rtc%ft_damping
1339 t0 = rtc%ft_t0
1340
1341 ! Determine field ft if polarizability required
1342 IF (do_polarizability) THEN
1343 ALLOCATE (field_results(3, n))
1344 IF (rtc%apply_delta_pulse) THEN
1345 ! Constant real FT
1346 IF (PRESENT(cell)) THEN
1347 delta_vec(:) = (real(rtc%delta_pulse_direction(1), kind=dp)*cell%h_inv(1, :) + &
1348 REAL(rtc%delta_pulse_direction(2), kind=dp)*cell%h_inv(2, :) + &
1349 REAL(rtc%delta_pulse_direction(3), kind=dp)*cell%h_inv(3, :)) &
1350 *twopi*rtc%delta_pulse_scale
1351 ELSE
1352 delta_vec(:) = real(rtc%delta_pulse_direction(:), kind=dp)*rtc%delta_pulse_scale
1353 END IF
1354 DO k = 1, 3
1355 field_results(k, :) = cmplx(delta_vec(k), 0.0, kind=dp)
1356 END DO
1357 ELSE
1358 ! Do explicit FT
1359 CALL multi_fft(times, fields, field_results, &
1360 damping_opt=damping, t0_opt=t0, subtract_initial_opt=.true.)
1361 END IF
1362 END IF
1363
1364 IF (do_moments_ft .OR. do_polarizability) THEN
1365 ! We need to transform at least the moments
1366 ! NOTE : Might be able to save some memory by only doing FT of actually
1367 ! required moments, but for now, doing FT of all moment directions
1368 ALLOCATE (results(3*nspin, n))
1369 ALLOCATE (omegas(n))
1370 ALLOCATE (value_series(3*nspin, n))
1371 DO i = 1, nspin
1372 DO k = 1, 3
1373 value_series(3*(i - 1) + k, :) = moments(i, k, :)
1374 END DO
1375 END DO
1376 ! TODO : Choose whether the initial subtraction is applied in &FT section?
1377 CALL multi_fft(times, value_series, results, omegas, &
1378 damping_opt=damping, t0_opt=t0, subtract_initial_opt=.true.)
1379 DEALLOCATE (value_series)
1380 DO i = 1, nspin
1381 ! Output to FT file, if needed
1382 ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension=file_extensions(i), &
1383 file_form="FORMATTED", file_position="REWIND")
1384 ! Print header
1385 IF (ft_unit > 0) THEN
1386 ALLOCATE (headers(7))
1387 headers(2) = " x,real [at.u.]"
1388 headers(3) = " x,imag [at.u.]"
1389 headers(4) = " y,real [at.u.]"
1390 headers(5) = " y,imag [at.u.]"
1391 headers(6) = " z,real [at.u.]"
1392 headers(7) = " z,imag [at.u.]"
1393 IF (info_unit == ft_unit) THEN
1394 headers(1) = "# Energy [eV]"
1395 prefix = " MOMENTS_FT|"
1396 prefix_format = "(A12)"
1397 CALL print_rt_file(ft_unit, headers, omegas, results(3*(i - 1) + 1:3*(i - 1) + 3, :), &
1398 prefix, prefix_format, evolt)
1399 ELSE
1400 headers(1) = "# omega [at.u.]"
1401 CALL print_rt_file(ft_unit, headers, omegas, results(3*(i - 1) + 1:3*(i - 1) + 3, :))
1402 END IF
1403 DEALLOCATE (headers)
1404 END IF
1405 CALL cp_print_key_finished_output(ft_unit, logger, moment_ft_section)
1406 END DO
1407 ! Spin-summed total moments FT (open shell only; inert for nspin=1)
1408 IF (nspin > 1) THEN
1409 ALLOCATE (results_spin_total(3, n))
1410 results_spin_total(:, :) = (0.0_dp, 0.0_dp)
1411 DO i = 1, nspin
1412 DO k = 1, 3
1413 results_spin_total(k, :) = results_spin_total(k, :) + results(3*(i - 1) + k, :)
1414 END DO
1415 END DO
1416 ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension="_SPIN_TOTAL.dat", &
1417 file_form="FORMATTED", file_position="REWIND")
1418 IF (ft_unit > 0) THEN
1419 ALLOCATE (headers(7))
1420 headers(2) = " x,real [at.u.]"
1421 headers(3) = " x,imag [at.u.]"
1422 headers(4) = " y,real [at.u.]"
1423 headers(5) = " y,imag [at.u.]"
1424 headers(6) = " z,real [at.u.]"
1425 headers(7) = " z,imag [at.u.]"
1426 IF (info_unit == ft_unit) THEN
1427 headers(1) = "# Energy [eV]"
1428 prefix = " MOMENTS_FT|"
1429 prefix_format = "(A12)"
1430 CALL print_rt_file(ft_unit, headers, omegas, results_spin_total, &
1431 prefix, prefix_format, evolt)
1432 ELSE
1433 headers(1) = "# omega [at.u.]"
1434 CALL print_rt_file(ft_unit, headers, omegas, results_spin_total)
1435 END IF
1436 DEALLOCATE (headers)
1437 END IF
1438 CALL cp_print_key_finished_output(ft_unit, logger, moment_ft_section)
1439 DEALLOCATE (results_spin_total)
1440 END IF
1441 END IF
1442
1443 IF (rtc%pade_requested .AND. (do_moments_ft .OR. do_polarizability)) THEN
1444 ALLOCATE (omegas_complex(SIZE(omegas)))
1445 omegas_complex(:) = cmplx(omegas(:), 0.0, kind=dp)
1446 n_pade = int((rtc%pade_e_max - rtc%pade_e_min)/rtc%pade_e_step)
1447 ALLOCATE (omegas_pade(n_pade))
1448 ALLOCATE (omegas_pade_real(n_pade))
1449 ! Construct omegas_pade and omegas_complex
1450 DO i = 1, n_pade
1451 omegas_pade_real(i) = (i - 1)*rtc%pade_e_step + rtc%pade_e_min
1452 omegas_pade(i) = cmplx(omegas_pade_real(i), 0.0, kind=dp)
1453 END DO
1454 ALLOCATE (results_pade(nspin*3, n_pade), source=cmplx(0.0, 0.0, kind=dp))
1455 DO i = 1, nspin
1456 DO k = 1, 3
1457 CALL greenx_refine_ft(rtc%pade_fit_e_min, rtc%pade_fit_e_max, omegas_complex, results(3*(i - 1) + k, :), &
1458 omegas_pade, results_pade(3*(i - 1) + k, :))
1459 END DO
1460 ! Print to a file
1461 ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension="_PADE"//file_extensions(i), &
1462 file_form="FORMATTED", file_position="REWIND")
1463 IF (ft_unit > 0) THEN
1464 ALLOCATE (headers(7))
1465 headers(2) = " x,real,pade [at.u.]"
1466 headers(3) = " x,imag,pade [at.u.]"
1467 headers(4) = " y,real,pade [at.u.]"
1468 headers(5) = " y,imag,pade [at.u.]"
1469 headers(6) = " z,real,pade [at.u.]"
1470 headers(7) = " z,imag,pade [at.u.]"
1471 IF (info_unit == ft_unit) THEN
1472 headers(1) = "# Energy [eV]"
1473 prefix = " MOMENTS_FT_PADE|"
1474 prefix_format = "(A17)"
1475 CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade(3*(i - 1) + 1:3*(i - 1) + 3, :), &
1476 prefix, prefix_format, evolt)
1477 ELSE
1478 headers(1) = "# omega [at.u.]"
1479 CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade(3*(i - 1) + 1:3*(i - 1) + 3, :))
1480 END IF
1481 DEALLOCATE (headers)
1482 END IF
1483 END DO
1484 ! Spin-summed total moments-FT Padé (open shell only; inert for nspin=1)
1485 IF (nspin > 1) THEN
1486 ALLOCATE (results_pade_spin_total(3, n_pade))
1487 results_pade_spin_total(:, :) = (0.0_dp, 0.0_dp)
1488 DO i = 1, nspin
1489 DO k = 1, 3
1490 results_pade_spin_total(k, :) = results_pade_spin_total(k, :) + results_pade(3*(i - 1) + k, :)
1491 END DO
1492 END DO
1493 ft_unit = cp_print_key_unit_nr(logger, moment_ft_section, extension="_PADE_SPIN_TOTAL.dat", &
1494 file_form="FORMATTED", file_position="REWIND")
1495 IF (ft_unit > 0) THEN
1496 ALLOCATE (headers(7))
1497 headers(2) = " x,real,pade [at.u.]"
1498 headers(3) = " x,imag,pade [at.u.]"
1499 headers(4) = " y,real,pade [at.u.]"
1500 headers(5) = " y,imag,pade [at.u.]"
1501 headers(6) = " z,real,pade [at.u.]"
1502 headers(7) = " z,imag,pade [at.u.]"
1503 IF (info_unit == ft_unit) THEN
1504 headers(1) = "# Energy [eV]"
1505 prefix = " MOMENTS_FT_PADE|"
1506 prefix_format = "(A17)"
1507 CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade_spin_total, &
1508 prefix, prefix_format, evolt)
1509 ELSE
1510 headers(1) = "# omega [at.u.]"
1511 CALL print_rt_file(ft_unit, headers, omegas_pade_real, results_pade_spin_total)
1512 END IF
1513 DEALLOCATE (headers)
1514 END IF
1515 DEALLOCATE (results_pade_spin_total)
1516 END IF
1517 END IF
1518
1519 IF (do_polarizability) THEN
1520 ! get the polarizability elements, as required
1521 ALLOCATE (pol_results(n_elems, n))
1522 DO i = 1, nspin
1523 DO k = 1, n_elems
1524 ! NOTE - field is regularized to small value
1525 pol_results(k, :) = results(3*(i - 1) + &
1526 rtc%print_pol_elements(k, 1), :)/ &
1527 (field_results(rtc%print_pol_elements(k, 2), :) + &
1528 1.0e-10*field_results(rtc%print_pol_elements(k, 2), 2))
1529 END DO
1530 ! Print to the file
1531 ft_unit = cp_print_key_unit_nr(logger, pol_section, extension=file_extensions(i), &
1532 file_form="FORMATTED", file_position="REWIND")
1533 IF (ft_unit > 0) THEN
1534 ALLOCATE (headers(2*n_elems + 1))
1535 DO k = 1, n_elems
1536 WRITE (headers(2*k), "(A16,I2,I2)") "real pol. elem.", &
1537 rtc%print_pol_elements(k, 1), &
1538 rtc%print_pol_elements(k, 2)
1539 WRITE (headers(2*k + 1), "(A16,I2,I2)") "imag pol. elem.", &
1540 rtc%print_pol_elements(k, 1), &
1541 rtc%print_pol_elements(k, 2)
1542 END DO
1543 ! Write header
1544 IF (info_unit == ft_unit) THEN
1545 headers(1) = "# Energy [eV]"
1546 prefix = " POLARIZABILITY|"
1547 prefix_format = "(A16)"
1548 CALL print_rt_file(ft_unit, headers, omegas, pol_results, &
1549 prefix, prefix_format, evolt)
1550 ELSE
1551 headers(1) = "# omega [at.u.]"
1552 CALL print_rt_file(ft_unit, headers, omegas, pol_results)
1553 END IF
1554 DEALLOCATE (headers)
1555 END IF
1556 CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1557 ! Static polarizability alpha(0): pol_results at the FFT-grid omega
1558 ! closest to zero. Re is alpha(0); Im should be machine-zero (sanity).
1559 IF (info_unit > 0) THEN
1560 idx_omega_zero = minloc(abs(omegas), dim=1)
1561 IF (i == 1) THEN
1562 WRITE (info_unit, '(A,T22,A,T28,A,T36,A,T59,A)') &
1563 " STATIC_POL|", "spin", "element", "Re [a.u.]", "Im [a.u.]"
1564 END IF
1565 DO k_static = 1, n_elems
1566 WRITE (info_unit, '(A,T22,I4,T28,I3,",",I3,T36,ES22.10E3,T59,ES22.10E3)') &
1567 " STATIC_POL|", i, &
1568 rtc%print_pol_elements(k_static, 1), &
1569 rtc%print_pol_elements(k_static, 2), &
1570 REAL(pol_results(k_static, idx_omega_zero), kind=dp), &
1571 aimag(pol_results(k_static, idx_omega_zero))
1572 END DO
1573 END IF
1574 END DO
1575 ! Spin-summed total polarizability (open shell only; inert for nspin=1).
1576 ! Field is spin-independent, so (sum_s moments_s)/field == sum_s (moments_s/field).
1577 IF (nspin > 1) THEN
1578 ALLOCATE (pol_results_spin_total(n_elems, n))
1579 pol_results_spin_total(:, :) = (0.0_dp, 0.0_dp)
1580 DO k = 1, n_elems
1581 DO i = 1, nspin
1582 pol_results_spin_total(k, :) = pol_results_spin_total(k, :) + &
1583 results(3*(i - 1) + rtc%print_pol_elements(k, 1), :)
1584 END DO
1585 pol_results_spin_total(k, :) = pol_results_spin_total(k, :)/ &
1586 (field_results(rtc%print_pol_elements(k, 2), :) + &
1587 1.0e-10*field_results(rtc%print_pol_elements(k, 2), 2))
1588 END DO
1589 ft_unit = cp_print_key_unit_nr(logger, pol_section, extension="_SPIN_TOTAL.dat", &
1590 file_form="FORMATTED", file_position="REWIND")
1591 IF (ft_unit > 0) THEN
1592 ALLOCATE (headers(2*n_elems + 1))
1593 DO k = 1, n_elems
1594 WRITE (headers(2*k), "(A16,I2,I2)") "real pol. elem.", &
1595 rtc%print_pol_elements(k, 1), &
1596 rtc%print_pol_elements(k, 2)
1597 WRITE (headers(2*k + 1), "(A16,I2,I2)") "imag pol. elem.", &
1598 rtc%print_pol_elements(k, 1), &
1599 rtc%print_pol_elements(k, 2)
1600 END DO
1601 IF (info_unit == ft_unit) THEN
1602 headers(1) = "# Energy [eV]"
1603 prefix = " POLARIZABILITY|"
1604 prefix_format = "(A16)"
1605 CALL print_rt_file(ft_unit, headers, omegas, pol_results_spin_total, &
1606 prefix, prefix_format, evolt)
1607 ELSE
1608 headers(1) = "# omega [at.u.]"
1609 CALL print_rt_file(ft_unit, headers, omegas, pol_results_spin_total)
1610 END IF
1611 DEALLOCATE (headers)
1612 END IF
1613 CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1614 ! Static polarizability total (header row already emitted by the per-spin block)
1615 IF (info_unit > 0) THEN
1616 idx_omega_zero = minloc(abs(omegas), dim=1)
1617 DO k_static = 1, n_elems
1618 WRITE (info_unit, '(A,T22,A,T28,I3,",",I3,T36,ES22.10E3,T59,ES22.10E3)') &
1619 " STATIC_POL|", "TOT", &
1620 rtc%print_pol_elements(k_static, 1), &
1621 rtc%print_pol_elements(k_static, 2), &
1622 REAL(pol_results_spin_total(k_static, idx_omega_zero), kind=dp), &
1623 aimag(pol_results_spin_total(k_static, idx_omega_zero))
1624 END DO
1625 END IF
1626 DEALLOCATE (pol_results_spin_total)
1627 END IF
1628 END IF
1629
1630 ! Padé polarizability
1631 IF (rtc%pade_requested .AND. do_polarizability) THEN
1632 ! Start with the field pade
1633 ALLOCATE (field_results_pade(3, n_pade))
1634 IF (rtc%apply_delta_pulse) THEN
1635 DO k = 1, 3
1636 field_results_pade(k, :) = cmplx(delta_vec(k), 0.0, kind=dp)
1637 END DO
1638 ELSE
1639 DO k = 1, 3
1640 CALL greenx_refine_ft(rtc%pade_fit_e_min, rtc%pade_fit_e_max, &
1641 omegas_complex, field_results(k, :), &
1642 omegas_pade, field_results_pade(k, :))
1643 END DO
1644 END IF
1645 ! Allocate polarisation pade
1646 ALLOCATE (pol_results_pade(n_elems, n_pade))
1647 ! Refine
1648 DO i = 1, nspin
1649 DO k = 1, n_elems
1650 ! NOTE : Regularization to small value
1651 pol_results_pade(k, :) = results_pade(3*(i - 1) + rtc%print_pol_elements(k, 1), :)/( &
1652 field_results_pade(rtc%print_pol_elements(k, 2), :) + &
1653 field_results_pade(rtc%print_pol_elements(k, 2), 2)*1.0e-10_dp)
1654 END DO
1655 ! Print to the file
1656 ft_unit = cp_print_key_unit_nr(logger, pol_section, extension="_PADE"//file_extensions(i), &
1657 file_form="FORMATTED", file_position="REWIND")
1658 IF (ft_unit > 0) THEN
1659 ALLOCATE (headers(2*n_elems + 1))
1660 DO k = 1, n_elems
1661 WRITE (headers(2*k), "(A16,I2,I2)") "re,pade,pol.", &
1662 rtc%print_pol_elements(k, 1), &
1663 rtc%print_pol_elements(k, 2)
1664 WRITE (headers(2*k + 1), "(A16,I2,I2)") "im,pade,pol.", &
1665 rtc%print_pol_elements(k, 1), &
1666 rtc%print_pol_elements(k, 2)
1667 END DO
1668 ! Write header
1669 IF (info_unit == ft_unit) THEN
1670 headers(1) = "# Energy [eV]"
1671 prefix = " POLARIZABILITY_PADE|"
1672 prefix_format = "(A21)"
1673 CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade, &
1674 prefix, prefix_format, evolt)
1675 ELSE
1676 headers(1) = "# omega [at.u.]"
1677 CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade)
1678 END IF
1679 DEALLOCATE (headers)
1680 END IF
1681 CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1682 END DO
1683 ! Spin-summed total Padé polarizability (open shell only; inert for nspin=1)
1684 IF (nspin > 1) THEN
1685 ALLOCATE (pol_results_pade_spin_total(n_elems, n_pade))
1686 pol_results_pade_spin_total(:, :) = (0.0_dp, 0.0_dp)
1687 DO k = 1, n_elems
1688 DO i = 1, nspin
1689 pol_results_pade_spin_total(k, :) = pol_results_pade_spin_total(k, :) + &
1690 results_pade(3*(i - 1) + rtc%print_pol_elements(k, 1), :)
1691 END DO
1692 pol_results_pade_spin_total(k, :) = pol_results_pade_spin_total(k, :)/( &
1693 field_results_pade(rtc%print_pol_elements(k, 2), :) + &
1694 field_results_pade(rtc%print_pol_elements(k, 2), 2)*1.0e-10_dp)
1695 END DO
1696 ft_unit = cp_print_key_unit_nr(logger, pol_section, extension="_PADE_SPIN_TOTAL.dat", &
1697 file_form="FORMATTED", file_position="REWIND")
1698 IF (ft_unit > 0) THEN
1699 ALLOCATE (headers(2*n_elems + 1))
1700 DO k = 1, n_elems
1701 WRITE (headers(2*k), "(A16,I2,I2)") "re,pade,pol.", &
1702 rtc%print_pol_elements(k, 1), &
1703 rtc%print_pol_elements(k, 2)
1704 WRITE (headers(2*k + 1), "(A16,I2,I2)") "im,pade,pol.", &
1705 rtc%print_pol_elements(k, 1), &
1706 rtc%print_pol_elements(k, 2)
1707 END DO
1708 IF (info_unit == ft_unit) THEN
1709 headers(1) = "# Energy [eV]"
1710 prefix = " POLARIZABILITY_PADE|"
1711 prefix_format = "(A21)"
1712 CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade_spin_total, &
1713 prefix, prefix_format, evolt)
1714 ELSE
1715 headers(1) = "# omega [at.u.]"
1716 CALL print_rt_file(ft_unit, headers, omegas_pade_real, pol_results_pade_spin_total)
1717 END IF
1718 DEALLOCATE (headers)
1719 END IF
1720 CALL cp_print_key_finished_output(ft_unit, logger, pol_section)
1721 DEALLOCATE (pol_results_pade_spin_total)
1722 END IF
1723 DEALLOCATE (field_results_pade)
1724 DEALLOCATE (pol_results_pade)
1725 END IF
1726
1727 IF (rtc%pade_requested .AND. (do_moments_ft .OR. do_polarizability)) THEN
1728 DEALLOCATE (omegas_complex)
1729 DEALLOCATE (omegas_pade)
1730 DEALLOCATE (omegas_pade_real)
1731 DEALLOCATE (results_pade)
1732 END IF
1733
1734 IF (do_polarizability) THEN
1735 DEALLOCATE (pol_results)
1736 DEALLOCATE (field_results)
1737 END IF
1738
1739 IF (do_moments_ft .OR. do_polarizability) THEN
1740 DEALLOCATE (results)
1741 DEALLOCATE (omegas)
1742 END IF
1743
1744 CALL timestop(handle)
1745
1746 END SUBROUTINE print_ft
1747
1748! **************************************************************************************************
1749!> \brief ...
1750!> \param rt_unit ...
1751!> \param headers ...
1752!> \param xvals ...
1753!> \param yvals ...
1754!> \param prefix ...
1755!> \param prefix_format ...
1756!> \param xscale_opt ...
1757!> \param comp_opt ...
1758! **************************************************************************************************
1759 SUBROUTINE print_rt_file(rt_unit, headers, xvals, yvals, prefix, prefix_format, xscale_opt, comp_opt)
1760 INTEGER, INTENT(IN) :: rt_unit
1761 CHARACTER(len=20), DIMENSION(:), INTENT(IN), &
1762 OPTIONAL :: headers
1763 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: xvals
1764 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: yvals
1765 CHARACTER(len=21), INTENT(IN), OPTIONAL :: prefix
1766 CHARACTER(len=5), INTENT(IN), OPTIONAL :: prefix_format
1767 REAL(kind=dp), INTENT(IN), OPTIONAL :: xscale_opt
1768 INTEGER, OPTIONAL :: comp_opt
1769
1770 INTEGER :: do_comp, i, j, ncols, nrows
1771 LOGICAL :: do_headers, do_prefix
1772 REAL(kind=dp) :: xscale
1773
1774 do_prefix = .false.
1775 IF (PRESENT(prefix)) THEN
1776 IF (PRESENT(prefix_format)) THEN
1777 do_prefix = .true.
1778 ELSE
1779 cpabort("Printing of prefix with missing format!")
1780 END IF
1781 END IF
1782
1783 xscale = 1.0_dp
1784 IF (PRESENT(xscale_opt)) xscale = xscale_opt
1785
1786 ncols = SIZE(yvals, 1)
1787 nrows = SIZE(yvals, 2)
1788
1789 ! Check whether printing complex data (default) or just a component
1790 do_comp = rt_file_comp_both
1791 IF (PRESENT(comp_opt)) do_comp = comp_opt
1792
1793 do_headers = PRESENT(headers)
1794 ! Check whether enough headers for yvals and xvals is present
1795 IF (do_headers) THEN
1796 IF ((do_comp /= rt_file_comp_both .AND. SIZE(headers) < ncols + 1) &
1797 .OR. (do_comp == rt_file_comp_both .AND. SIZE(headers) < 2*ncols + 1)) THEN
1798 cpabort("Not enought headers to print the file!")
1799 END IF
1800 END IF
1801
1802 IF (SIZE(xvals) < nrows) THEN
1803 cpabort("Not enough xvals to print all yvals!")
1804 END IF
1805
1806 IF (rt_unit > 0) THEN
1807 ! Print headers
1808 IF (do_headers) THEN
1809 ! If prefix is present, write prefix
1810 IF (do_prefix) THEN
1811 WRITE (rt_unit, prefix_format, advance="no") prefix
1812 END IF
1813 WRITE (rt_unit, "(A20)", advance="no") headers(1)
1814 ! Print the rest of the headers
1815 SELECT CASE (do_comp)
1816 CASE (rt_file_comp_both)
1817 ! Complex case
1818 DO j = 1, 2*ncols - 1
1819 WRITE (rt_unit, "(A20)", advance="no") headers(j + 1)
1820 END DO
1821 WRITE (rt_unit, "(A20)") headers(2*ncols + 1)
1822 CASE DEFAULT
1823 ! For other cases, just one component is printed
1824 DO j = 1, ncols - 1
1825 WRITE (rt_unit, "(A20)", advance="no") headers(j + 1)
1826 END DO
1827 WRITE (rt_unit, "(A20)") headers(ncols + 1)
1828 END SELECT
1829 END IF
1830 ! Done with the headers, print actual data
1831 DO i = 1, nrows
1832 ! If prefix is present, write prefix
1833 IF (do_prefix) THEN
1834 WRITE (rt_unit, prefix_format, advance="no") prefix
1835 END IF
1836 WRITE (rt_unit, "(E20.8E3)", advance="no") xvals(i)*xscale
1837 DO j = 1, ncols - 1
1838 SELECT CASE (do_comp)
1839 CASE (rt_file_comp_real)
1840 WRITE (rt_unit, "(E20.8E3)", advance="no") &
1841 REAL(yvals(j, i))
1842 CASE (rt_file_comp_imag)
1843 WRITE (rt_unit, "(E20.8E3)", advance="no") &
1844 aimag(yvals(j, i))
1845 CASE DEFAULT
1846 ! Print both components
1847 WRITE (rt_unit, "(E20.8E3,E20.8E3)", advance="no") &
1848 REAL(yvals(j, i)), aimag(yvals(j, i))
1849 END SELECT
1850 END DO
1851 ! Print the final column(s)
1852 SELECT CASE (do_comp)
1853 CASE (rt_file_comp_real)
1854 WRITE (rt_unit, "(E20.8E3)") real(yvals(j, i))
1855 CASE (rt_file_comp_imag)
1856 WRITE (rt_unit, "(E20.8E3)") aimag(yvals(j, i))
1857 CASE DEFAULT
1858 ! Print both components
1859 WRITE (rt_unit, "(E20.8E3,E20.8E3)") &
1860 REAL(yvals(j, i)), aimag(yvals(j, i))
1861 END SELECT
1862 END DO
1863 END IF
1864 END SUBROUTINE print_rt_file
1865
1866END MODULE rt_propagation_output
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
kahan_sum::accurate_sum
Definition kahan_sum.F:41
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
pw_methods::pw_zero
Definition pw_methods.F:83
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:658
cp_dbcsr_api::dbcsr_binary_write
subroutine, public dbcsr_binary_write(matrix, filepath)
...
Definition cp_dbcsr_api.F:290
cp_dbcsr_api::dbcsr_get_occupation
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
Definition cp_dbcsr_api.F:893
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_checksum
real(kind=dp) function, public dbcsr_checksum(matrix, pos)
Calculates the checksum of a DBCSR matrix.
Definition cp_dbcsr_contrib.F:628
cp_dbcsr_contrib::dbcsr_trace
subroutine, public dbcsr_trace(matrix, trace)
Computes the trace of the given matrix, also known as the sum of its diagonal elements.
Definition cp_dbcsr_contrib.F:412
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_double
subroutine, public cp_fm_struct_double(fmstruct, struct, context, col, row)
creates a struct with twice the number of blocks on each core. If matrix A has to be multiplied with ...
Definition cp_fm_struct.F:527
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_iter_string
character(len=default_string_length) function, public cp_iter_string(iter_info, print_key, for_file)
returns the iteration string, a string that is useful to create unique filenames (once you trim it)
Definition cp_output_handling.F:517
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_printkey_is_on
logical function, public cp_printkey_is_on(iteration_info, print_key)
returns true if the printlevel activates this printkey does not look if this iteration it should be p...
Definition cp_output_handling.F:429
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:48
efield_utils
all routins needed for a nonperiodic electric field
Definition efield_utils.F:12
efield_utils::make_field
subroutine, public make_field(dft_control, field, sim_step, sim_time)
computes the amplitude of the efield within a given envelop
Definition efield_utils.F:120
greenx_interface
Interface to the Greenx library.
Definition greenx_interface.F:13
greenx_interface::greenx_refine_ft
subroutine, public greenx_refine_ft(fit_e_min, fit_e_max, x_fit, y_fit, x_eval, y_eval, n_pade_opt)
Refines the FT grid using Padé approximants.
Definition greenx_interface.F:205
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::ehrenfest
integer, parameter, public ehrenfest
Definition input_constants.F:131
input_constants::real_time_propagation
integer, parameter, public real_time_propagation
Definition input_constants.F:131
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition input_section_types.F:1001
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::one
real(kind=dp), parameter, public one
Definition mathconstants.F:123
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::zero
real(kind=dp), parameter, public zero
Definition mathconstants.F:123
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::femtoseconds
real(kind=dp), parameter, public femtoseconds
Definition physcon.F:153
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_linres_current
given the response wavefunctions obtained by the application of the (rxp), p, and ((dk-dl)xp) operato...
Definition qs_linres_current.F:17
qs_linres_current::calculate_jrho_resp
subroutine, public calculate_jrho_resp(mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, ib, idir, current_rs, current_gs, qs_env, current_env, soft_valid, retain_rsgrid)
Calculation of the idir component of the response current density in the presence of a constant magne...
Definition qs_linres_current.F:578
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_mo_io
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
qs_mo_io::write_rt_mos_to_restart
subroutine, public write_rt_mos_to_restart(mo_array, rt_mos, particle_set, dft_section, qs_kind_set)
...
Definition qs_mo_io.F:244
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:593
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_operators_ao
Definition qs_operators_ao.F:13
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition qs_operators_ao.F:74
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_post_gpw
Does all kind of post scf calculations for GPW/GAPW.
Definition qs_scf_post_gpw.F:15
qs_scf_post_gpw::write_mo_free_results
subroutine, public write_mo_free_results(qs_env)
Write QS results always available (if switched on through the print_keys) Can be called from ls_scf.
Definition qs_scf_post_gpw.F:2555
qs_scf_post_gpw::qs_scf_post_moments
subroutine, public qs_scf_post_moments(input, logger, qs_env, output_unit)
Computes and prints electric moments.
Definition qs_scf_post_gpw.F:1629
qs_scf_post_gpw::write_mo_dependent_results
subroutine, public write_mo_dependent_results(qs_env, scf_env)
Write QS results available if MO's are present (if switched on through the print_keys) Writes only MO...
Definition qs_scf_post_gpw.F:2220
qs_scf_post_tb
Does all kind of post scf calculations for DFTB.
Definition qs_scf_post_tb.F:15
qs_scf_post_tb::scf_post_calculation_tb
subroutine, public scf_post_calculation_tb(qs_env, tb_type, no_mos)
collects possible post - scf calculations and prints info / computes properties.
Definition qs_scf_post_tb.F:171
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:137
rt_projection_mo_utils
Function related to MO projection in RTP calculations.
Definition rt_projection_mo_utils.F:12
rt_projection_mo_utils::compute_and_write_proj_mo
subroutine, public compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
Compute the projection of the current MO coefficients on reference ones and write the results.
Definition rt_projection_mo_utils.F:352
rt_propagation_ft
Separation of Fourier transform utilities into separate file.
Definition rt_propagation_ft.F:12
rt_propagation_ft::multi_fft
subroutine, public multi_fft(time_series, value_series, result_series, omega_series, damping_opt, t0_opt, subtract_initial_opt)
Calculates the Fourier transform - couples to FFT libraries in CP2K, if available.
Definition rt_propagation_ft.F:100
rt_propagation_output
Routine for the real time propagation output.
Definition rt_propagation_output.F:13
rt_propagation_output::print_ft
subroutine, public print_ft(rtp_section, moments, times, fields, rtc, info_opt, cell)
Calculate and print the Fourier transforms + polarizabilites from moment trace.
Definition rt_propagation_output.F:1284
rt_propagation_output::report_density_occupation
subroutine, public report_density_occupation(filter_eps, rho)
Reports the sparsity pattern of the complex density matrix.
Definition rt_propagation_output.F:722
rt_propagation_output::print_rt_file
subroutine, public print_rt_file(rt_unit, headers, xvals, yvals, prefix, prefix_format, xscale_opt, comp_opt)
...
Definition rt_propagation_output.F:1760
rt_propagation_output::rt_prop_output
subroutine, public rt_prop_output(qs_env, run_type, delta_iter, used_time)
...
Definition rt_propagation_output.F:135
rt_propagation_output::print_moments
subroutine, public print_moments(moments_section, info_unit, moments, time, imag_opt, append_opt)
Print the dipole moments into a file.
Definition rt_propagation_output.F:1115
rt_propagation_output::rt_convergence
subroutine, public rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
computes the convergence criterion for RTP and EMD
Definition rt_propagation_output.F:510
rt_propagation_output::rt_file_comp_both
integer, parameter, public rt_file_comp_both
Definition rt_propagation_output.F:121
rt_propagation_output::calc_local_moment
subroutine, public calc_local_moment(moment_matrices, density_matrices, work, moment, imag_opt)
Calculate the values of real/imaginary parts of moments in all directions.
Definition rt_propagation_output.F:1239
rt_propagation_output::rt_file_comp_imag
integer, parameter, public rt_file_comp_imag
Definition rt_propagation_output.F:121
rt_propagation_output::rt_file_comp_real
integer, parameter, public rt_file_comp_real
Definition rt_propagation_output.F:121
rt_propagation_output::rt_convergence_density
subroutine, public rt_convergence_density(rtp, delta_p, delta_eps)
computes the convergence criterion for RTP and EMD based on the density matrix
Definition rt_propagation_output.F:621
rt_propagation_types
Types and set_get for real time propagation depending on runtype and diagonalization method different...
Definition rt_propagation_types.F:21
rt_propagation_types::get_rtp
subroutine, public get_rtp(rtp, exp_h_old, exp_h_new, h_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, s_inv, s_half, s_minus_half, b_mat, c_mat, propagator_matrix, mixing, mixing_factor, s_der, dt, nsteps, sinvh, sinvh_imag, sinvb, admm_mos)
...
Definition rt_propagation_types.F:416
rt_propagation_utils
Routines needed for EMD.
Definition rt_propagation_utils.F:13
rt_propagation_utils::write_rtp_mos_to_output_unit
subroutine, public write_rtp_mos_to_output_unit(qs_env, rtp)
...
Definition rt_propagation_utils.F:635
rt_propagation_utils::write_rtp_mo_cubes
subroutine, public write_rtp_mo_cubes(qs_env, rtp)
Write the time dependent amplitude of the MOs in real grid. Very close to qs_scf_post_gpw/qs_scf_post...
Definition rt_propagation_utils.F:735
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_control_types::rtp_control_type
Definition cp_control_types.F:94
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_comm_type
Definition message_passing.F:161
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_linres_types::current_env_type
Definition qs_linres_types.F:94
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76