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