Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 cp_control_types, ONLY: dft_control_type
16 : USE cp_dbcsr_api, ONLY: &
17 : dbcsr_add, dbcsr_binary_write, dbcsr_checksum, dbcsr_copy, dbcsr_create, &
18 : dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_filter, dbcsr_get_info, &
19 : dbcsr_get_occupation, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
20 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, &
21 : dbcsr_set, dbcsr_type
22 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
23 : dbcsr_allocate_matrix_set,&
24 : dbcsr_deallocate_matrix_set
25 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
26 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27 : cp_fm_struct_double,&
28 : cp_fm_struct_release,&
29 : cp_fm_struct_type
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_get_info,&
32 : cp_fm_release,&
33 : cp_fm_set_all,&
34 : cp_fm_type
35 : USE cp_log_handling, ONLY: cp_get_default_logger,&
36 : cp_logger_get_default_io_unit,&
37 : cp_logger_get_default_unit_nr,&
38 : cp_logger_type
39 : USE cp_output_handling, ONLY: cp_iter_string,&
40 : cp_p_file,&
41 : cp_print_key_finished_output,&
42 : cp_print_key_should_output,&
43 : cp_print_key_unit_nr
44 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
45 : USE efield_utils, ONLY: make_field
46 : USE input_constants, ONLY: ehrenfest,&
47 : real_time_propagation
48 : USE input_section_types, ONLY: section_get_ivals,&
49 : section_vals_get_subs_vals,&
50 : section_vals_type
51 : USE kahan_sum, ONLY: accurate_sum
52 : USE kinds, ONLY: default_path_length,&
53 : dp
54 : USE machine, ONLY: m_flush
55 : USE message_passing, ONLY: mp_comm_type
56 : USE parallel_gemm_api, ONLY: parallel_gemm
57 : USE particle_list_types, ONLY: particle_list_type
58 : USE particle_types, ONLY: particle_type
59 : USE physcon, ONLY: femtoseconds
60 : USE pw_env_types, ONLY: pw_env_get,&
61 : pw_env_type
62 : USE pw_methods, ONLY: pw_zero
63 : USE pw_pool_types, ONLY: pw_pool_type
64 : USE pw_types, ONLY: pw_c1d_gs_type,&
65 : pw_r3d_rs_type
66 : USE qs_energy_types, ONLY: qs_energy_type
67 : USE qs_environment_types, ONLY: get_qs_env,&
68 : qs_environment_type
69 : USE qs_kind_types, ONLY: get_qs_kind_set,&
70 : qs_kind_type
71 : USE qs_linres_current, ONLY: calculate_jrho_resp
72 : USE qs_linres_types, ONLY: current_env_type
73 : USE qs_mo_io, ONLY: write_rt_mos_to_restart
74 : USE qs_rho_types, ONLY: qs_rho_get,&
75 : qs_rho_type
76 : USE qs_scf_post_gpw, ONLY: qs_scf_post_moments,&
77 : write_mo_dependent_results,&
78 : write_mo_free_results
79 : USE qs_scf_post_tb, ONLY: scf_post_calculation_tb
80 : USE qs_scf_types, ONLY: qs_scf_env_type
81 : USE qs_subsys_types, ONLY: qs_subsys_get,&
82 : qs_subsys_type
83 : USE rt_projection_mo_utils, ONLY: compute_and_write_proj_mo
84 : USE rt_propagation_types, ONLY: get_rtp,&
85 : rt_prop_type
86 : USE rt_propagation_utils, ONLY: calculate_P_imaginary,&
87 : write_rtp_mo_cubes,&
88 : write_rtp_mos_to_output_unit
89 : #include "../base/base_uses.f90"
90 :
91 : IMPLICIT NONE
92 :
93 : PRIVATE
94 :
95 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_output'
96 :
97 : PUBLIC :: rt_prop_output, &
98 : rt_convergence, &
99 : rt_convergence_density, &
100 : report_density_occupation
101 :
102 : CONTAINS
103 :
104 : ! **************************************************************************************************
105 : !> \brief ...
106 : !> \param qs_env ...
107 : !> \param run_type ...
108 : !> \param delta_iter ...
109 : !> \param used_time ...
110 : ! **************************************************************************************************
111 4588 : SUBROUTINE rt_prop_output(qs_env, run_type, delta_iter, used_time)
112 : TYPE(qs_environment_type), POINTER :: qs_env
113 : INTEGER, INTENT(in) :: run_type
114 : REAL(dp), INTENT(in), OPTIONAL :: delta_iter, used_time
115 :
116 : INTEGER :: n_electrons, n_proj, nspin, output_unit, &
117 : spin
118 : REAL(dp) :: orthonormality, tot_rho_r
119 2294 : REAL(KIND=dp), DIMENSION(:), POINTER :: qs_tot_rho_r
120 2294 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
121 2294 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
122 : TYPE(cp_logger_type), POINTER :: logger
123 2294 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, P_im, rho_new
124 : TYPE(dft_control_type), POINTER :: dft_control
125 2294 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
126 2294 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
127 : TYPE(qs_rho_type), POINTER :: rho
128 : TYPE(rt_prop_type), POINTER :: rtp
129 : TYPE(section_vals_type), POINTER :: dft_section, input, rtp_section
130 :
131 2294 : NULLIFY (logger, dft_control)
132 :
133 4588 : logger => cp_get_default_logger()
134 : CALL get_qs_env(qs_env, &
135 : rtp=rtp, &
136 : matrix_s=matrix_s, &
137 : input=input, &
138 : rho=rho, &
139 : particle_set=particle_set, &
140 : atomic_kind_set=atomic_kind_set, &
141 : qs_kind_set=qs_kind_set, &
142 2294 : dft_control=dft_control)
143 :
144 2294 : rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
145 :
146 2294 : CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
147 2294 : n_electrons = n_electrons - dft_control%charge
148 :
149 2294 : CALL qs_rho_get(rho_struct=rho, tot_rho_r=qs_tot_rho_r)
150 :
151 2294 : tot_rho_r = accurate_sum(qs_tot_rho_r)
152 :
153 : output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
154 2294 : extension=".scfLog")
155 :
156 2294 : IF (output_unit > 0) THEN
157 : WRITE (output_unit, FMT="(/,(T3,A,T40,I5))") &
158 1147 : "Information at iteration step:", rtp%iter
159 : WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
160 1147 : "Total electronic density (r-space): ", &
161 1147 : tot_rho_r, &
162 : tot_rho_r + &
163 2294 : REAL(n_electrons, dp)
164 : WRITE (UNIT=output_unit, FMT="((T3,A,T59,F22.14))") &
165 1147 : "Total energy:", rtp%energy_new
166 1147 : IF (run_type == ehrenfest) &
167 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
168 575 : "Energy difference to previous iteration step:", rtp%energy_new - rtp%energy_old
169 1147 : IF (run_type == real_time_propagation) &
170 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
171 572 : "Energy difference to initial state:", rtp%energy_new - rtp%energy_old
172 1147 : IF (PRESENT(delta_iter)) &
173 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,E20.6))") &
174 1147 : "Convergence:", delta_iter
175 1147 : IF (rtp%converged) THEN
176 307 : IF (run_type == real_time_propagation) &
177 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F12.2))") &
178 172 : "Time needed for propagation:", used_time
179 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,3X,F16.14))") &
180 307 : "CONVERGENCE REACHED", rtp%energy_new - rtp%energy_old
181 : END IF
182 : END IF
183 :
184 2294 : IF (rtp%converged) THEN
185 614 : IF (.NOT. rtp%linear_scaling) THEN
186 432 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
187 : CALL rt_calculate_orthonormality(orthonormality, &
188 432 : mos_new, matrix_s(1)%matrix)
189 432 : IF (output_unit > 0) &
190 : WRITE (output_unit, FMT="(/,(T3,A,T60,F20.10))") &
191 216 : "Max deviation from orthonormalization:", orthonormality
192 : END IF
193 : END IF
194 :
195 2294 : IF (output_unit > 0) &
196 1147 : CALL m_flush(output_unit)
197 : CALL cp_print_key_finished_output(output_unit, logger, rtp_section, &
198 2294 : "PRINT%PROGRAM_RUN_INFO")
199 :
200 2294 : IF (rtp%converged) THEN
201 614 : dft_section => section_vals_get_subs_vals(input, "DFT")
202 614 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
203 : dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) &
204 30 : CALL print_field_applied(qs_env, dft_section)
205 614 : CALL make_moment(qs_env)
206 614 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
207 : dft_section, "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS"), cp_p_file)) THEN
208 4 : CALL print_rtp_energy_components(qs_env, dft_section)
209 : END IF
210 614 : IF (.NOT. dft_control%qs_control%dftb) THEN
211 494 : CALL write_available_results(qs_env=qs_env, rtp=rtp)
212 : END IF
213 614 : IF (rtp%linear_scaling) THEN
214 182 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
215 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
216 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART"), cp_p_file)) THEN
217 96 : CALL write_rt_p_to_restart(rho_new, .FALSE.)
218 : END IF
219 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
220 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY"), cp_p_file)) THEN
221 2 : CALL write_rt_p_to_restart(rho_new, .TRUE.)
222 : END IF
223 182 : IF (.NOT. dft_control%qs_control%dftb) THEN
224 : !Not sure if these things could also work with dftb or not
225 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
226 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
227 64 : DO spin = 1, SIZE(rho_new)/2
228 64 : CALL rt_current(qs_env, rho_new(2*spin)%matrix, dft_section, spin, SIZE(rho_new)/2)
229 : END DO
230 : END IF
231 : END IF
232 : ELSE
233 432 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
234 432 : IF (.NOT. dft_control%qs_control%dftb) THEN
235 312 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
236 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
237 12 : NULLIFY (P_im)
238 12 : nspin = SIZE(mos_new)/2
239 12 : CALL dbcsr_allocate_matrix_set(P_im, nspin)
240 24 : DO spin = 1, nspin
241 12 : CALL dbcsr_init_p(P_im(spin)%matrix)
242 24 : CALL dbcsr_create(P_im(spin)%matrix, template=matrix_s(1)%matrix, matrix_type="N")
243 : END DO
244 12 : CALL calculate_P_imaginary(qs_env, rtp, P_im)
245 24 : DO spin = 1, nspin
246 24 : CALL rt_current(qs_env, P_im(spin)%matrix, dft_section, spin, nspin)
247 : END DO
248 12 : CALL dbcsr_deallocate_matrix_set(P_im)
249 : END IF
250 312 : IF (dft_control%rtp_control%is_proj_mo) THEN
251 112 : DO n_proj = 1, SIZE(dft_control%rtp_control%proj_mo_list)
252 : CALL compute_and_write_proj_mo(qs_env, mos_new, &
253 112 : dft_control%rtp_control%proj_mo_list(n_proj)%proj_mo, n_proj)
254 : END DO
255 : END IF
256 : END IF
257 : CALL write_rt_mos_to_restart(qs_env%mos, mos_new, particle_set, &
258 432 : dft_section, qs_kind_set)
259 : END IF
260 : END IF
261 :
262 2294 : rtp%energy_old = rtp%energy_new
263 :
264 2294 : IF (.NOT. rtp%converged .AND. rtp%iter >= dft_control%rtp_control%max_iter) &
265 : CALL cp_abort(__LOCATION__, "EMD did not converge, either increase MAX_ITER "// &
266 0 : "or use a smaller TIMESTEP")
267 :
268 2294 : END SUBROUTINE rt_prop_output
269 :
270 : ! **************************************************************************************************
271 : !> \brief computes the effective orthonormality of a set of mos given an s-matrix
272 : !> orthonormality is the max deviation from unity of the C^T S C
273 : !> \param orthonormality ...
274 : !> \param mos_new ...
275 : !> \param matrix_s ...
276 : !> \author Florian Schiffmann (02.09)
277 : ! **************************************************************************************************
278 432 : SUBROUTINE rt_calculate_orthonormality(orthonormality, mos_new, matrix_s)
279 : REAL(KIND=dp), INTENT(out) :: orthonormality
280 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
281 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
282 :
283 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_calculate_orthonormality'
284 :
285 : INTEGER :: handle, i, im, ispin, j, k, n, &
286 : ncol_local, nrow_local, nspin, re
287 432 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
288 : REAL(KIND=dp) :: alpha, max_alpha, max_beta
289 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
290 : TYPE(cp_fm_type) :: overlap_re, svec_im, svec_re
291 :
292 432 : NULLIFY (tmp_fm_struct)
293 :
294 432 : CALL timeset(routineN, handle)
295 :
296 432 : nspin = SIZE(mos_new)/2
297 432 : max_alpha = 0.0_dp
298 432 : max_beta = 0.0_dp
299 966 : DO ispin = 1, nspin
300 534 : re = ispin*2 - 1
301 534 : im = ispin*2
302 : ! get S*C
303 534 : CALL cp_fm_create(svec_re, mos_new(im)%matrix_struct)
304 534 : CALL cp_fm_create(svec_im, mos_new(im)%matrix_struct)
305 : CALL cp_fm_get_info(mos_new(im), &
306 534 : nrow_global=n, ncol_global=k)
307 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(re), &
308 534 : svec_re, k)
309 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(im), &
310 534 : svec_im, k)
311 :
312 : ! get C^T (S*C)
313 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
314 : para_env=mos_new(re)%matrix_struct%para_env, &
315 534 : context=mos_new(re)%matrix_struct%context)
316 534 : CALL cp_fm_create(overlap_re, tmp_fm_struct)
317 :
318 534 : CALL cp_fm_struct_release(tmp_fm_struct)
319 :
320 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(re), &
321 534 : svec_re, 0.0_dp, overlap_re)
322 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(im), &
323 534 : svec_im, 1.0_dp, overlap_re)
324 :
325 534 : CALL cp_fm_release(svec_re)
326 534 : CALL cp_fm_release(svec_im)
327 :
328 : CALL cp_fm_get_info(overlap_re, nrow_local=nrow_local, ncol_local=ncol_local, &
329 534 : row_indices=row_indices, col_indices=col_indices)
330 1413 : DO i = 1, nrow_local
331 5020 : DO j = 1, ncol_local
332 3607 : alpha = overlap_re%local_data(i, j)
333 3607 : IF (row_indices(i) .EQ. col_indices(j)) alpha = alpha - 1.0_dp
334 4486 : max_alpha = MAX(max_alpha, ABS(alpha))
335 : END DO
336 : END DO
337 2568 : CALL cp_fm_release(overlap_re)
338 : END DO
339 432 : CALL mos_new(1)%matrix_struct%para_env%max(max_alpha)
340 432 : CALL mos_new(1)%matrix_struct%para_env%max(max_beta)
341 432 : orthonormality = max_alpha
342 :
343 432 : CALL timestop(handle)
344 :
345 432 : END SUBROUTINE rt_calculate_orthonormality
346 :
347 : ! **************************************************************************************************
348 : !> \brief computes the convergence criterion for RTP and EMD
349 : !> \param rtp ...
350 : !> \param matrix_s Overlap matrix without the derivatives
351 : !> \param delta_mos ...
352 : !> \param delta_eps ...
353 : !> \author Florian Schiffmann (02.09)
354 : ! **************************************************************************************************
355 :
356 1478 : SUBROUTINE rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
357 : TYPE(rt_prop_type), POINTER :: rtp
358 : TYPE(dbcsr_type), POINTER :: matrix_s
359 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: delta_mos
360 : REAL(dp), INTENT(out) :: delta_eps
361 :
362 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence'
363 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
364 :
365 : INTEGER :: handle, i, icol, im, ispin, j, lcol, &
366 : lrow, nao, newdim, nmo, nspin, re
367 : LOGICAL :: double_col, double_row
368 : REAL(KIND=dp) :: alpha, max_alpha
369 : TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1, tmp_fm_struct
370 : TYPE(cp_fm_type) :: work, work1, work2
371 1478 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
372 :
373 1478 : NULLIFY (tmp_fm_struct)
374 :
375 1478 : CALL timeset(routineN, handle)
376 :
377 1478 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
378 :
379 1478 : nspin = SIZE(delta_mos)/2
380 1478 : max_alpha = 0.0_dp
381 :
382 5258 : DO i = 1, SIZE(mos_new)
383 5258 : CALL cp_fm_scale_and_add(-one, delta_mos(i), one, mos_new(i))
384 : END DO
385 :
386 3368 : DO ispin = 1, nspin
387 1890 : re = ispin*2 - 1
388 1890 : im = ispin*2
389 :
390 1890 : double_col = .TRUE.
391 1890 : double_row = .FALSE.
392 : CALL cp_fm_struct_double(newstruct, &
393 : delta_mos(re)%matrix_struct, &
394 : delta_mos(re)%matrix_struct%context, &
395 : double_col, &
396 1890 : double_row)
397 :
398 1890 : CALL cp_fm_create(work, matrix_struct=newstruct)
399 1890 : CALL cp_fm_create(work1, matrix_struct=newstruct)
400 :
401 : CALL cp_fm_get_info(delta_mos(re), ncol_local=lcol, ncol_global=nmo, &
402 1890 : nrow_global=nao)
403 1890 : CALL cp_fm_get_info(work, ncol_global=newdim)
404 :
405 1890 : CALL cp_fm_set_all(work, zero, zero)
406 :
407 8132 : DO icol = 1, lcol
408 51268 : work%local_data(:, icol) = delta_mos(re)%local_data(:, icol)
409 53158 : work%local_data(:, icol + lcol) = delta_mos(im)%local_data(:, icol)
410 : END DO
411 :
412 1890 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, work, work1, ncol=newdim)
413 :
414 1890 : CALL cp_fm_release(work)
415 :
416 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, &
417 : para_env=delta_mos(re)%matrix_struct%para_env, &
418 1890 : context=delta_mos(re)%matrix_struct%context)
419 : CALL cp_fm_struct_double(newstruct1, &
420 : tmp_fm_struct, &
421 : delta_mos(re)%matrix_struct%context, &
422 : double_col, &
423 1890 : double_row)
424 :
425 1890 : CALL cp_fm_create(work, matrix_struct=newstruct1)
426 1890 : CALL cp_fm_create(work2, matrix_struct=newstruct1)
427 :
428 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(re), &
429 1890 : work1, zero, work)
430 :
431 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(im), &
432 1890 : work1, zero, work2)
433 :
434 1890 : CALL cp_fm_get_info(work, nrow_local=lrow)
435 5011 : DO i = 1, lrow
436 17796 : DO j = 1, lcol
437 : alpha = SQRT((work%local_data(i, j) + work2%local_data(i, j + lcol))**2 + &
438 12785 : (work%local_data(i, j + lcol) - work2%local_data(i, j))**2)
439 15906 : max_alpha = MAX(max_alpha, ABS(alpha))
440 : END DO
441 : END DO
442 :
443 1890 : CALL cp_fm_release(work)
444 1890 : CALL cp_fm_release(work1)
445 1890 : CALL cp_fm_release(work2)
446 1890 : CALL cp_fm_struct_release(tmp_fm_struct)
447 1890 : CALL cp_fm_struct_release(newstruct)
448 9038 : CALL cp_fm_struct_release(newstruct1)
449 :
450 : END DO
451 :
452 1478 : CALL delta_mos(1)%matrix_struct%para_env%max(max_alpha)
453 1478 : delta_eps = SQRT(max_alpha)
454 :
455 1478 : CALL timestop(handle)
456 :
457 1478 : END SUBROUTINE rt_convergence
458 :
459 : ! **************************************************************************************************
460 : !> \brief computes the convergence criterion for RTP and EMD based on the density matrix
461 : !> \param rtp ...
462 : !> \param delta_P ...
463 : !> \param delta_eps ...
464 : !> \author Samuel Andermatt (02.14)
465 : ! **************************************************************************************************
466 :
467 1632 : SUBROUTINE rt_convergence_density(rtp, delta_P, delta_eps)
468 :
469 : TYPE(rt_prop_type), POINTER :: rtp
470 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P
471 : REAL(dp), INTENT(out) :: delta_eps
472 :
473 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence_density'
474 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
475 :
476 : INTEGER :: col_atom, group_handle, handle, i, &
477 : ispin, row_atom
478 : REAL(dp) :: alpha, max_alpha
479 816 : REAL(dp), DIMENSION(:), POINTER :: block_values
480 : TYPE(dbcsr_iterator_type) :: iter
481 816 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
482 : TYPE(dbcsr_type), POINTER :: tmp
483 : TYPE(mp_comm_type) :: group
484 :
485 816 : CALL timeset(routineN, handle)
486 :
487 816 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
488 :
489 2940 : DO i = 1, SIZE(rho_new)
490 2940 : CALL dbcsr_add(delta_P(i)%matrix, rho_new(i)%matrix, one, -one)
491 : END DO
492 : !get the maximum value of delta_P
493 2940 : DO i = 1, SIZE(delta_P)
494 : !square all entries of both matrices
495 2124 : CALL dbcsr_iterator_start(iter, delta_P(i)%matrix)
496 12356 : DO WHILE (dbcsr_iterator_blocks_left(iter))
497 10232 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
498 712140 : block_values = block_values*block_values
499 : END DO
500 5064 : CALL dbcsr_iterator_stop(iter)
501 : END DO
502 : NULLIFY (tmp)
503 816 : ALLOCATE (tmp)
504 816 : CALL dbcsr_create(tmp, template=delta_P(1)%matrix, matrix_type="N")
505 1878 : DO ispin = 1, SIZE(delta_P)/2
506 1062 : CALL dbcsr_desymmetrize(delta_P(2*ispin - 1)%matrix, tmp)
507 1878 : CALL dbcsr_add(delta_P(2*ispin)%matrix, tmp, one, one)
508 : END DO
509 : !the absolute values are now in the even entries of delta_P
510 816 : max_alpha = zero
511 1878 : DO ispin = 1, SIZE(delta_P)/2
512 1062 : CALL dbcsr_iterator_start(iter, delta_P(2*ispin)%matrix)
513 6275 : DO WHILE (dbcsr_iterator_blocks_left(iter))
514 5213 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
515 365095 : alpha = MAXVAL(block_values)
516 6275 : IF (alpha > max_alpha) max_alpha = alpha
517 : END DO
518 2940 : CALL dbcsr_iterator_stop(iter)
519 : END DO
520 816 : CALL dbcsr_get_info(delta_P(1)%matrix, group=group_handle)
521 816 : CALL group%set_handle(group_handle)
522 816 : CALL group%max(max_alpha)
523 816 : delta_eps = SQRT(max_alpha)
524 816 : CALL dbcsr_deallocate_matrix(tmp)
525 816 : CALL timestop(handle)
526 :
527 816 : END SUBROUTINE rt_convergence_density
528 :
529 : ! **************************************************************************************************
530 : !> \brief interface to qs_moments. Does only work for nonperiodic dipole
531 : !> \param qs_env ...
532 : !> \author Florian Schiffmann (02.09)
533 : ! **************************************************************************************************
534 :
535 614 : SUBROUTINE make_moment(qs_env)
536 :
537 : TYPE(qs_environment_type), POINTER :: qs_env
538 :
539 : CHARACTER(len=*), PARAMETER :: routineN = 'make_moment'
540 :
541 : INTEGER :: handle, output_unit
542 : TYPE(cp_logger_type), POINTER :: logger
543 : TYPE(dft_control_type), POINTER :: dft_control
544 :
545 614 : CALL timeset(routineN, handle)
546 :
547 614 : NULLIFY (dft_control)
548 :
549 614 : logger => cp_get_default_logger()
550 614 : output_unit = cp_logger_get_default_io_unit(logger)
551 614 : CALL get_qs_env(qs_env, dft_control=dft_control)
552 614 : IF (dft_control%qs_control%dftb) THEN
553 120 : CALL scf_post_calculation_tb(qs_env, "DFTB", .FALSE.)
554 494 : ELSE IF (dft_control%qs_control%xtb) THEN
555 60 : CALL scf_post_calculation_tb(qs_env, "xTB", .FALSE.)
556 : ELSE
557 434 : CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
558 : END IF
559 614 : CALL timestop(handle)
560 :
561 614 : END SUBROUTINE make_moment
562 :
563 : ! **************************************************************************************************
564 : !> \brief Reports the sparsity pattern of the complex density matrix
565 : !> \param filter_eps ...
566 : !> \param rho ...
567 : !> \author Samuel Andermatt (09.14)
568 : ! **************************************************************************************************
569 :
570 182 : SUBROUTINE report_density_occupation(filter_eps, rho)
571 :
572 : REAL(KIND=dp) :: filter_eps
573 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho
574 :
575 : CHARACTER(len=*), PARAMETER :: routineN = 'report_density_occupation'
576 :
577 : INTEGER :: handle, i, im, ispin, re, unit_nr
578 : REAL(KIND=dp) :: eps, occ
579 : TYPE(cp_logger_type), POINTER :: logger
580 182 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmp
581 :
582 182 : CALL timeset(routineN, handle)
583 :
584 182 : logger => cp_get_default_logger()
585 182 : unit_nr = cp_logger_get_default_io_unit(logger)
586 182 : NULLIFY (tmp)
587 182 : CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
588 686 : DO i = 1, SIZE(rho)
589 504 : CALL dbcsr_init_p(tmp(i)%matrix)
590 504 : CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
591 686 : CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
592 : END DO
593 434 : DO ispin = 1, SIZE(rho)/2
594 252 : re = 2*ispin - 1
595 252 : im = 2*ispin
596 252 : eps = MAX(filter_eps, 1.0E-11_dp)
597 2338 : DO WHILE (eps < 1.1_dp)
598 2086 : CALL dbcsr_filter(tmp(re)%matrix, eps)
599 2086 : occ = dbcsr_get_occupation(tmp(re)%matrix)
600 3129 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
601 2086 : ispin, " eps ", eps, " real: ", occ
602 2086 : eps = eps*10
603 : END DO
604 252 : eps = MAX(filter_eps, 1.0E-11_dp)
605 2520 : DO WHILE (eps < 1.1_dp)
606 2086 : CALL dbcsr_filter(tmp(im)%matrix, eps)
607 2086 : occ = dbcsr_get_occupation(tmp(im)%matrix)
608 3129 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
609 2086 : ispin, " eps ", eps, " imag: ", occ
610 2086 : eps = eps*10.0_dp
611 : END DO
612 : END DO
613 182 : CALL dbcsr_deallocate_matrix_set(tmp)
614 182 : CALL timestop(handle)
615 :
616 182 : END SUBROUTINE report_density_occupation
617 :
618 : ! **************************************************************************************************
619 : !> \brief Writes the density matrix and the atomic positions to a restart file
620 : !> \param rho_new ...
621 : !> \param history ...
622 : !> \author Samuel Andermatt (09.14)
623 : ! **************************************************************************************************
624 :
625 98 : SUBROUTINE write_rt_p_to_restart(rho_new, history)
626 :
627 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
628 : LOGICAL :: history
629 :
630 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_p_to_restart'
631 :
632 : CHARACTER(LEN=default_path_length) :: file_name, project_name
633 : INTEGER :: handle, im, ispin, re, unit_nr
634 : REAL(KIND=dp) :: cs_pos
635 : TYPE(cp_logger_type), POINTER :: logger
636 :
637 98 : CALL timeset(routineN, handle)
638 98 : logger => cp_get_default_logger()
639 98 : IF (logger%para_env%is_source()) THEN
640 49 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
641 : ELSE
642 : unit_nr = -1
643 : END IF
644 :
645 98 : project_name = logger%iter_info%project_name
646 234 : DO ispin = 1, SIZE(rho_new)/2
647 136 : re = 2*ispin - 1
648 136 : im = 2*ispin
649 136 : IF (history) THEN
650 : WRITE (file_name, '(A,I0,A)') &
651 2 : TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
652 : ELSE
653 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
654 : END IF
655 136 : cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.TRUE.)
656 136 : IF (unit_nr > 0) THEN
657 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
658 : END IF
659 136 : CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
660 136 : IF (history) THEN
661 : WRITE (file_name, '(A,I0,A)') &
662 2 : TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
663 : ELSE
664 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
665 : END IF
666 136 : cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.TRUE.)
667 136 : IF (unit_nr > 0) THEN
668 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
669 : END IF
670 234 : CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
671 : END DO
672 :
673 98 : CALL timestop(handle)
674 :
675 98 : END SUBROUTINE write_rt_p_to_restart
676 :
677 : ! **************************************************************************************************
678 : !> \brief Collocation of the current and printing of it in a cube file
679 : !> \param qs_env ...
680 : !> \param P_im ...
681 : !> \param dft_section ...
682 : !> \param spin ...
683 : !> \param nspin ...
684 : !> \author Samuel Andermatt (06.15)
685 : ! **************************************************************************************************
686 48 : SUBROUTINE rt_current(qs_env, P_im, dft_section, spin, nspin)
687 : TYPE(qs_environment_type), POINTER :: qs_env
688 : TYPE(dbcsr_type), POINTER :: P_im
689 : TYPE(section_vals_type), POINTER :: dft_section
690 : INTEGER :: spin, nspin
691 :
692 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_current'
693 :
694 : CHARACTER(len=1) :: char_spin
695 : CHARACTER(len=14) :: ext
696 : CHARACTER(len=2) :: sdir
697 : INTEGER :: dir, handle, print_unit
698 48 : INTEGER, DIMENSION(:), POINTER :: stride(:)
699 : LOGICAL :: mpi_io
700 : TYPE(cp_logger_type), POINTER :: logger
701 : TYPE(current_env_type) :: current_env
702 : TYPE(dbcsr_type), POINTER :: tmp, zero
703 : TYPE(particle_list_type), POINTER :: particles
704 : TYPE(pw_c1d_gs_type) :: gs
705 : TYPE(pw_env_type), POINTER :: pw_env
706 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
707 : TYPE(pw_r3d_rs_type) :: rs
708 : TYPE(qs_subsys_type), POINTER :: subsys
709 :
710 48 : CALL timeset(routineN, handle)
711 :
712 48 : logger => cp_get_default_logger()
713 48 : CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
714 48 : CALL qs_subsys_get(subsys, particles=particles)
715 48 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
716 :
717 48 : NULLIFY (zero, tmp)
718 48 : ALLOCATE (zero, tmp)
719 48 : CALL dbcsr_create(zero, template=P_im)
720 48 : CALL dbcsr_copy(zero, P_im)
721 48 : CALL dbcsr_set(zero, 0.0_dp)
722 48 : CALL dbcsr_create(tmp, template=P_im)
723 48 : CALL dbcsr_copy(tmp, P_im)
724 48 : IF (nspin == 1) THEN
725 32 : CALL dbcsr_scale(tmp, 0.5_dp)
726 : END IF
727 48 : current_env%gauge = -1
728 48 : current_env%gauge_init = .FALSE.
729 48 : CALL auxbas_pw_pool%create_pw(rs)
730 48 : CALL auxbas_pw_pool%create_pw(gs)
731 :
732 : NULLIFY (stride)
733 48 : ALLOCATE (stride(3))
734 :
735 192 : DO dir = 1, 3
736 :
737 144 : CALL pw_zero(rs)
738 144 : CALL pw_zero(gs)
739 :
740 144 : CALL calculate_jrho_resp(zero, tmp, zero, zero, dir, dir, rs, gs, qs_env, current_env, retain_rsgrid=.TRUE.)
741 :
742 576 : stride = section_get_ivals(dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT%STRIDE")
743 :
744 144 : IF (dir == 1) THEN
745 48 : sdir = "-x"
746 96 : ELSEIF (dir == 2) THEN
747 48 : sdir = "-y"
748 : ELSE
749 48 : sdir = "-z"
750 : END IF
751 144 : WRITE (char_spin, "(I1)") spin
752 :
753 144 : ext = "-SPIN-"//char_spin//sdir//".cube"
754 144 : mpi_io = .TRUE.
755 : print_unit = cp_print_key_unit_nr(logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
756 : extension=ext, file_status="REPLACE", file_action="WRITE", &
757 144 : log_filename=.FALSE., mpi_io=mpi_io)
758 :
759 : CALL cp_pw_to_cube(rs, print_unit, "EMD current", particles=particles, stride=stride, &
760 144 : mpi_io=mpi_io)
761 :
762 : CALL cp_print_key_finished_output(print_unit, logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
763 192 : mpi_io=mpi_io)
764 :
765 : END DO
766 :
767 48 : CALL auxbas_pw_pool%give_back_pw(rs)
768 48 : CALL auxbas_pw_pool%give_back_pw(gs)
769 :
770 48 : CALL dbcsr_deallocate_matrix(zero)
771 48 : CALL dbcsr_deallocate_matrix(tmp)
772 :
773 48 : DEALLOCATE (stride)
774 :
775 48 : CALL timestop(handle)
776 :
777 3504 : END SUBROUTINE rt_current
778 :
779 : ! **************************************************************************************************
780 : !> \brief Interface routine to trigger writing of results available from normal
781 : !> SCF. Can write MO-dependent and MO free results (needed for call from
782 : !> the linear scaling code)
783 : !> Update: trigger also some of prints for time-dependent runs
784 : !> \param qs_env ...
785 : !> \param rtp ...
786 : !> \par History
787 : !> 2022-11 Update [Guillaume Le Breton]
788 : ! **************************************************************************************************
789 494 : SUBROUTINE write_available_results(qs_env, rtp)
790 : TYPE(qs_environment_type), POINTER :: qs_env
791 : TYPE(rt_prop_type), POINTER :: rtp
792 :
793 : CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results'
794 :
795 : INTEGER :: handle
796 : TYPE(qs_scf_env_type), POINTER :: scf_env
797 :
798 494 : CALL timeset(routineN, handle)
799 :
800 494 : CALL get_qs_env(qs_env, scf_env=scf_env)
801 494 : IF (rtp%linear_scaling) THEN
802 182 : CALL write_mo_free_results(qs_env)
803 : ELSE
804 312 : CALL write_mo_free_results(qs_env)
805 312 : CALL write_mo_dependent_results(qs_env, scf_env)
806 : ! Time-dependent MO print
807 312 : CALL write_rtp_mos_to_output_unit(qs_env, rtp)
808 312 : CALL write_rtp_mo_cubes(qs_env, rtp)
809 : END IF
810 :
811 494 : CALL timestop(handle)
812 :
813 494 : END SUBROUTINE write_available_results
814 :
815 : ! **************************************************************************************************
816 : !> \brief Print the field applied to the system. Either the electric
817 : !> field or the vector potential depending on the gauge used
818 : !> \param qs_env ...
819 : !> \param dft_section ...
820 : !> \par History
821 : !> 2023-01 Created [Guillaume Le Breton]
822 : ! **************************************************************************************************
823 30 : SUBROUTINE print_field_applied(qs_env, dft_section)
824 : TYPE(qs_environment_type), POINTER :: qs_env
825 : TYPE(section_vals_type), POINTER :: dft_section
826 :
827 : CHARACTER(LEN=3), DIMENSION(3) :: rlab
828 : CHARACTER(LEN=default_path_length) :: filename
829 : INTEGER :: i, i_step, output_unit, unit_nr
830 : LOGICAL :: new_file
831 : REAL(kind=dp) :: field(3)
832 : TYPE(cp_logger_type), POINTER :: logger
833 : TYPE(dft_control_type), POINTER :: dft_control
834 : TYPE(rt_prop_type), POINTER :: rtp
835 :
836 30 : NULLIFY (dft_control)
837 :
838 30 : logger => cp_get_default_logger()
839 30 : output_unit = cp_logger_get_default_io_unit(logger)
840 :
841 30 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp)
842 :
843 30 : i_step = rtp%istep
844 :
845 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
846 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat", is_new_file=new_file)
847 :
848 30 : IF (output_unit > 0) THEN
849 60 : rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
850 15 : IF (unit_nr /= output_unit) THEN
851 15 : INQUIRE (UNIT=unit_nr, NAME=filename)
852 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
853 15 : "FIELD", "The field applied is written to the file:", &
854 30 : TRIM(filename)
855 : ELSE
856 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED [a.u.]"
857 : WRITE (UNIT=output_unit, FMT="(T5,3(A,A,E16.8,1X))") &
858 0 : (TRIM(rlab(i)), "=", dft_control%rtp_control%field(i), i=1, 3)
859 : END IF
860 :
861 15 : IF (new_file) THEN
862 3 : IF (dft_control%apply_efield_field) THEN
863 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,3(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z"
864 2 : ELSE IF (dft_control%apply_vector_potential) THEN
865 0 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,6(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z", &
866 0 : " Vec. Pot. X", " Vec. Pot. Y", " Vec. Pot. Z"
867 : END IF
868 : END IF
869 :
870 15 : field = 0.0_dp
871 15 : IF (dft_control%apply_efield_field) THEN
872 4 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
873 4 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
874 8 : field(1), field(2), field(3)
875 : ! DO i=1,3
876 : ! IF (ABS(field(i))< 10E-10) field(i) = 0.0_dp
877 : ! END IF
878 11 : ELSE IF (dft_control%apply_vector_potential) THEN
879 9 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
880 9 : dft_control%rtp_control%field(1), dft_control%rtp_control%field(2), dft_control%rtp_control%field(3), &
881 18 : dft_control%rtp_control%vec_pot(1), dft_control%rtp_control%vec_pot(2), dft_control%rtp_control%vec_pot(3)
882 : END IF
883 :
884 : END IF
885 :
886 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
887 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD")
888 :
889 30 : END SUBROUTINE print_field_applied
890 :
891 : ! **************************************************************************************************
892 : !> \brief Print the components of the total energy used in an RTP calculation
893 : !> \param qs_env ...
894 : !> \param dft_section ...
895 : !> \par History
896 : !> 2024-02 Created [ANB]
897 : ! **************************************************************************************************
898 4 : SUBROUTINE print_rtp_energy_components(qs_env, dft_section)
899 : TYPE(qs_environment_type), POINTER :: qs_env
900 : TYPE(section_vals_type), POINTER :: dft_section
901 :
902 : CHARACTER(LEN=default_path_length) :: filename
903 : INTEGER :: i_step, output_unit, unit_nr
904 : LOGICAL :: new_file
905 : TYPE(cp_logger_type), POINTER :: logger
906 : TYPE(dft_control_type), POINTER :: dft_control
907 : TYPE(qs_energy_type), POINTER :: energy
908 : TYPE(rt_prop_type), POINTER :: rtp
909 :
910 4 : NULLIFY (dft_control, energy, rtp)
911 :
912 4 : logger => cp_get_default_logger()
913 4 : output_unit = cp_logger_get_default_io_unit(logger)
914 :
915 4 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
916 4 : i_step = rtp%istep
917 :
918 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
919 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS", extension=".ener", &
920 4 : file_action="WRITE", is_new_file=new_file)
921 :
922 4 : IF (output_unit > 0) THEN
923 2 : IF (unit_nr /= output_unit) THEN
924 2 : INQUIRE (UNIT=unit_nr, NAME=filename)
925 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
926 2 : "ENERGY_CONSTITUENTS", "Total Energy constituents written to file:", &
927 4 : TRIM(filename)
928 : ELSE
929 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ENERGY_CONSTITUENTS"
930 : END IF
931 :
932 2 : IF (new_file) THEN
933 : ! NOTE that these are not all terms contributing to the total energy for RTP, only a selection of those
934 : ! most significant / impactful. Therefore the printed components likely will not add up to the total energy.
935 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,10(6X,A))') "Step Nr.", "Time[fs]", &
936 1 : "Total ener.[a.u.]", "core[a.u.] ", " overlap [a.u.]", "hartree[a.u.]", " exc. [a.u.] ", &
937 2 : " hartree 1c[a.u.]", "exc 1c[a.u.] ", "exc admm[a.u.]", "exc 1c admm[a.u.]", "efield LG"
938 :
939 : END IF
940 : WRITE (UNIT=unit_nr, FMT="(I10,F20.6,10(F20.9))") &
941 2 : qs_env%sim_step, qs_env%sim_time*femtoseconds, &
942 2 : energy%total, energy%core, energy%core_overlap, energy%hartree, energy%exc, &
943 4 : energy%hartree_1c, energy%exc1, energy%exc_aux_fit, energy%exc1_aux_fit, energy%efield_core
944 :
945 : END IF
946 :
947 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
948 4 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS")
949 :
950 4 : END SUBROUTINE print_rtp_energy_components
951 :
952 : END MODULE rt_propagation_output
|
