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