Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_copy, dbcsr_create, dbcsr_deallocate_matrix, &
18 : dbcsr_desymmetrize, dbcsr_filter, dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, &
19 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
20 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_checksum
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, handle, i, ispin, row_atom
477 : REAL(dp) :: alpha, max_alpha
478 816 : REAL(dp), DIMENSION(:, :), POINTER :: block_values
479 : TYPE(dbcsr_iterator_type) :: iter
480 816 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
481 : TYPE(dbcsr_type), POINTER :: tmp
482 : TYPE(mp_comm_type) :: group
483 :
484 816 : CALL timeset(routineN, handle)
485 :
486 816 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
487 :
488 2940 : DO i = 1, SIZE(rho_new)
489 2940 : CALL dbcsr_add(delta_P(i)%matrix, rho_new(i)%matrix, one, -one)
490 : END DO
491 : !get the maximum value of delta_P
492 2940 : DO i = 1, SIZE(delta_P)
493 : !square all entries of both matrices
494 2124 : CALL dbcsr_iterator_start(iter, delta_P(i)%matrix)
495 12356 : DO WHILE (dbcsr_iterator_blocks_left(iter))
496 10232 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
497 792248 : block_values = block_values*block_values
498 : END DO
499 5064 : CALL dbcsr_iterator_stop(iter)
500 : END DO
501 : NULLIFY (tmp)
502 816 : ALLOCATE (tmp)
503 816 : CALL dbcsr_create(tmp, template=delta_P(1)%matrix, matrix_type="N")
504 1878 : DO ispin = 1, SIZE(delta_P)/2
505 1062 : CALL dbcsr_desymmetrize(delta_P(2*ispin - 1)%matrix, tmp)
506 1878 : CALL dbcsr_add(delta_P(2*ispin)%matrix, tmp, one, one)
507 : END DO
508 : !the absolute values are now in the even entries of delta_P
509 816 : max_alpha = zero
510 1878 : DO ispin = 1, SIZE(delta_P)/2
511 1062 : CALL dbcsr_iterator_start(iter, delta_P(2*ispin)%matrix)
512 6275 : DO WHILE (dbcsr_iterator_blocks_left(iter))
513 5213 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
514 400573 : alpha = MAXVAL(block_values)
515 6275 : IF (alpha > max_alpha) max_alpha = alpha
516 : END DO
517 2940 : CALL dbcsr_iterator_stop(iter)
518 : END DO
519 816 : CALL dbcsr_get_info(delta_P(1)%matrix, group=group)
520 816 : CALL group%max(max_alpha)
521 816 : delta_eps = SQRT(max_alpha)
522 816 : CALL dbcsr_deallocate_matrix(tmp)
523 816 : CALL timestop(handle)
524 :
525 816 : END SUBROUTINE rt_convergence_density
526 :
527 : ! **************************************************************************************************
528 : !> \brief interface to qs_moments. Does only work for nonperiodic dipole
529 : !> \param qs_env ...
530 : !> \author Florian Schiffmann (02.09)
531 : ! **************************************************************************************************
532 :
533 614 : SUBROUTINE make_moment(qs_env)
534 :
535 : TYPE(qs_environment_type), POINTER :: qs_env
536 :
537 : CHARACTER(len=*), PARAMETER :: routineN = 'make_moment'
538 :
539 : INTEGER :: handle, output_unit
540 : TYPE(cp_logger_type), POINTER :: logger
541 : TYPE(dft_control_type), POINTER :: dft_control
542 :
543 614 : CALL timeset(routineN, handle)
544 :
545 614 : NULLIFY (dft_control)
546 :
547 614 : logger => cp_get_default_logger()
548 614 : output_unit = cp_logger_get_default_io_unit(logger)
549 614 : CALL get_qs_env(qs_env, dft_control=dft_control)
550 614 : IF (dft_control%qs_control%dftb) THEN
551 120 : CALL scf_post_calculation_tb(qs_env, "DFTB", .FALSE.)
552 494 : ELSE IF (dft_control%qs_control%xtb) THEN
553 60 : CALL scf_post_calculation_tb(qs_env, "xTB", .FALSE.)
554 : ELSE
555 434 : CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
556 : END IF
557 614 : CALL timestop(handle)
558 :
559 614 : END SUBROUTINE make_moment
560 :
561 : ! **************************************************************************************************
562 : !> \brief Reports the sparsity pattern of the complex density matrix
563 : !> \param filter_eps ...
564 : !> \param rho ...
565 : !> \author Samuel Andermatt (09.14)
566 : ! **************************************************************************************************
567 :
568 182 : SUBROUTINE report_density_occupation(filter_eps, rho)
569 :
570 : REAL(KIND=dp) :: filter_eps
571 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho
572 :
573 : CHARACTER(len=*), PARAMETER :: routineN = 'report_density_occupation'
574 :
575 : INTEGER :: handle, i, im, ispin, re, unit_nr
576 : REAL(KIND=dp) :: eps, occ
577 : TYPE(cp_logger_type), POINTER :: logger
578 182 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmp
579 :
580 182 : CALL timeset(routineN, handle)
581 :
582 182 : logger => cp_get_default_logger()
583 182 : unit_nr = cp_logger_get_default_io_unit(logger)
584 182 : NULLIFY (tmp)
585 182 : CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
586 686 : DO i = 1, SIZE(rho)
587 504 : CALL dbcsr_init_p(tmp(i)%matrix)
588 504 : CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
589 686 : CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
590 : END DO
591 434 : DO ispin = 1, SIZE(rho)/2
592 252 : re = 2*ispin - 1
593 252 : im = 2*ispin
594 252 : eps = MAX(filter_eps, 1.0E-11_dp)
595 2338 : DO WHILE (eps < 1.1_dp)
596 2086 : CALL dbcsr_filter(tmp(re)%matrix, eps)
597 2086 : occ = dbcsr_get_occupation(tmp(re)%matrix)
598 3129 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
599 2086 : ispin, " eps ", eps, " real: ", occ
600 2086 : eps = eps*10
601 : END DO
602 252 : eps = MAX(filter_eps, 1.0E-11_dp)
603 2520 : DO WHILE (eps < 1.1_dp)
604 2086 : CALL dbcsr_filter(tmp(im)%matrix, eps)
605 2086 : occ = dbcsr_get_occupation(tmp(im)%matrix)
606 3129 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
607 2086 : ispin, " eps ", eps, " imag: ", occ
608 2086 : eps = eps*10.0_dp
609 : END DO
610 : END DO
611 182 : CALL dbcsr_deallocate_matrix_set(tmp)
612 182 : CALL timestop(handle)
613 :
614 182 : END SUBROUTINE report_density_occupation
615 :
616 : ! **************************************************************************************************
617 : !> \brief Writes the density matrix and the atomic positions to a restart file
618 : !> \param rho_new ...
619 : !> \param history ...
620 : !> \author Samuel Andermatt (09.14)
621 : ! **************************************************************************************************
622 :
623 98 : SUBROUTINE write_rt_p_to_restart(rho_new, history)
624 :
625 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
626 : LOGICAL :: history
627 :
628 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_p_to_restart'
629 :
630 : CHARACTER(LEN=default_path_length) :: file_name, project_name
631 : INTEGER :: handle, im, ispin, re, unit_nr
632 : REAL(KIND=dp) :: cs_pos
633 : TYPE(cp_logger_type), POINTER :: logger
634 :
635 98 : CALL timeset(routineN, handle)
636 98 : logger => cp_get_default_logger()
637 98 : IF (logger%para_env%is_source()) THEN
638 49 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
639 : ELSE
640 : unit_nr = -1
641 : END IF
642 :
643 98 : project_name = logger%iter_info%project_name
644 234 : DO ispin = 1, SIZE(rho_new)/2
645 136 : re = 2*ispin - 1
646 136 : im = 2*ispin
647 136 : IF (history) THEN
648 : WRITE (file_name, '(A,I0,A)') &
649 2 : TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
650 : ELSE
651 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
652 : END IF
653 136 : cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.TRUE.)
654 136 : IF (unit_nr > 0) THEN
655 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
656 : END IF
657 136 : CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
658 136 : IF (history) THEN
659 : WRITE (file_name, '(A,I0,A)') &
660 2 : TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
661 : ELSE
662 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
663 : END IF
664 136 : cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.TRUE.)
665 136 : IF (unit_nr > 0) THEN
666 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
667 : END IF
668 234 : CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
669 : END DO
670 :
671 98 : CALL timestop(handle)
672 :
673 98 : END SUBROUTINE write_rt_p_to_restart
674 :
675 : ! **************************************************************************************************
676 : !> \brief Collocation of the current and printing of it in a cube file
677 : !> \param qs_env ...
678 : !> \param P_im ...
679 : !> \param dft_section ...
680 : !> \param spin ...
681 : !> \param nspin ...
682 : !> \author Samuel Andermatt (06.15)
683 : ! **************************************************************************************************
684 48 : SUBROUTINE rt_current(qs_env, P_im, dft_section, spin, nspin)
685 : TYPE(qs_environment_type), POINTER :: qs_env
686 : TYPE(dbcsr_type), POINTER :: P_im
687 : TYPE(section_vals_type), POINTER :: dft_section
688 : INTEGER :: spin, nspin
689 :
690 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_current'
691 :
692 : CHARACTER(len=1) :: char_spin
693 : CHARACTER(len=14) :: ext
694 : CHARACTER(len=2) :: sdir
695 : INTEGER :: dir, handle, print_unit
696 48 : INTEGER, DIMENSION(:), POINTER :: stride(:)
697 : LOGICAL :: mpi_io
698 : TYPE(cp_logger_type), POINTER :: logger
699 : TYPE(current_env_type) :: current_env
700 : TYPE(dbcsr_type), POINTER :: tmp, zero
701 : TYPE(particle_list_type), POINTER :: particles
702 : TYPE(pw_c1d_gs_type) :: gs
703 : TYPE(pw_env_type), POINTER :: pw_env
704 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
705 : TYPE(pw_r3d_rs_type) :: rs
706 : TYPE(qs_subsys_type), POINTER :: subsys
707 :
708 48 : CALL timeset(routineN, handle)
709 :
710 48 : logger => cp_get_default_logger()
711 48 : CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
712 48 : CALL qs_subsys_get(subsys, particles=particles)
713 48 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
714 :
715 48 : NULLIFY (zero, tmp)
716 48 : ALLOCATE (zero, tmp)
717 48 : CALL dbcsr_create(zero, template=P_im)
718 48 : CALL dbcsr_copy(zero, P_im)
719 48 : CALL dbcsr_set(zero, 0.0_dp)
720 48 : CALL dbcsr_create(tmp, template=P_im)
721 48 : CALL dbcsr_copy(tmp, P_im)
722 48 : IF (nspin == 1) THEN
723 32 : CALL dbcsr_scale(tmp, 0.5_dp)
724 : END IF
725 48 : current_env%gauge = -1
726 48 : current_env%gauge_init = .FALSE.
727 48 : CALL auxbas_pw_pool%create_pw(rs)
728 48 : CALL auxbas_pw_pool%create_pw(gs)
729 :
730 : NULLIFY (stride)
731 48 : ALLOCATE (stride(3))
732 :
733 192 : DO dir = 1, 3
734 :
735 144 : CALL pw_zero(rs)
736 144 : CALL pw_zero(gs)
737 :
738 144 : CALL calculate_jrho_resp(zero, tmp, zero, zero, dir, dir, rs, gs, qs_env, current_env, retain_rsgrid=.TRUE.)
739 :
740 576 : stride = section_get_ivals(dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT%STRIDE")
741 :
742 144 : IF (dir == 1) THEN
743 48 : sdir = "-x"
744 96 : ELSEIF (dir == 2) THEN
745 48 : sdir = "-y"
746 : ELSE
747 48 : sdir = "-z"
748 : END IF
749 144 : WRITE (char_spin, "(I1)") spin
750 :
751 144 : ext = "-SPIN-"//char_spin//sdir//".cube"
752 144 : mpi_io = .TRUE.
753 : print_unit = cp_print_key_unit_nr(logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
754 : extension=ext, file_status="REPLACE", file_action="WRITE", &
755 144 : log_filename=.FALSE., mpi_io=mpi_io)
756 :
757 : CALL cp_pw_to_cube(rs, print_unit, "EMD current", particles=particles, stride=stride, &
758 144 : mpi_io=mpi_io)
759 :
760 : CALL cp_print_key_finished_output(print_unit, logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
761 192 : mpi_io=mpi_io)
762 :
763 : END DO
764 :
765 48 : CALL auxbas_pw_pool%give_back_pw(rs)
766 48 : CALL auxbas_pw_pool%give_back_pw(gs)
767 :
768 48 : CALL dbcsr_deallocate_matrix(zero)
769 48 : CALL dbcsr_deallocate_matrix(tmp)
770 :
771 48 : DEALLOCATE (stride)
772 :
773 48 : CALL timestop(handle)
774 :
775 3504 : END SUBROUTINE rt_current
776 :
777 : ! **************************************************************************************************
778 : !> \brief Interface routine to trigger writing of results available from normal
779 : !> SCF. Can write MO-dependent and MO free results (needed for call from
780 : !> the linear scaling code)
781 : !> Update: trigger also some of prints for time-dependent runs
782 : !> \param qs_env ...
783 : !> \param rtp ...
784 : !> \par History
785 : !> 2022-11 Update [Guillaume Le Breton]
786 : ! **************************************************************************************************
787 494 : SUBROUTINE write_available_results(qs_env, rtp)
788 : TYPE(qs_environment_type), POINTER :: qs_env
789 : TYPE(rt_prop_type), POINTER :: rtp
790 :
791 : CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results'
792 :
793 : INTEGER :: handle
794 : TYPE(qs_scf_env_type), POINTER :: scf_env
795 :
796 494 : CALL timeset(routineN, handle)
797 :
798 494 : CALL get_qs_env(qs_env, scf_env=scf_env)
799 494 : IF (rtp%linear_scaling) THEN
800 182 : CALL write_mo_free_results(qs_env)
801 : ELSE
802 312 : CALL write_mo_free_results(qs_env)
803 312 : CALL write_mo_dependent_results(qs_env, scf_env)
804 : ! Time-dependent MO print
805 312 : CALL write_rtp_mos_to_output_unit(qs_env, rtp)
806 312 : CALL write_rtp_mo_cubes(qs_env, rtp)
807 : END IF
808 :
809 494 : CALL timestop(handle)
810 :
811 494 : END SUBROUTINE write_available_results
812 :
813 : ! **************************************************************************************************
814 : !> \brief Print the field applied to the system. Either the electric
815 : !> field or the vector potential depending on the gauge used
816 : !> \param qs_env ...
817 : !> \param dft_section ...
818 : !> \par History
819 : !> 2023-01 Created [Guillaume Le Breton]
820 : ! **************************************************************************************************
821 30 : SUBROUTINE print_field_applied(qs_env, dft_section)
822 : TYPE(qs_environment_type), POINTER :: qs_env
823 : TYPE(section_vals_type), POINTER :: dft_section
824 :
825 : CHARACTER(LEN=3), DIMENSION(3) :: rlab
826 : CHARACTER(LEN=default_path_length) :: filename
827 : INTEGER :: i, i_step, output_unit, unit_nr
828 : LOGICAL :: new_file
829 : REAL(kind=dp) :: field(3)
830 : TYPE(cp_logger_type), POINTER :: logger
831 : TYPE(dft_control_type), POINTER :: dft_control
832 : TYPE(rt_prop_type), POINTER :: rtp
833 :
834 30 : NULLIFY (dft_control)
835 :
836 30 : logger => cp_get_default_logger()
837 30 : output_unit = cp_logger_get_default_io_unit(logger)
838 :
839 30 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp)
840 :
841 30 : i_step = rtp%istep
842 :
843 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
844 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat", is_new_file=new_file)
845 :
846 30 : IF (output_unit > 0) THEN
847 60 : rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
848 15 : IF (unit_nr /= output_unit) THEN
849 15 : INQUIRE (UNIT=unit_nr, NAME=filename)
850 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
851 15 : "FIELD", "The field applied is written to the file:", &
852 30 : TRIM(filename)
853 : ELSE
854 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED [a.u.]"
855 : WRITE (UNIT=output_unit, FMT="(T5,3(A,A,E16.8,1X))") &
856 0 : (TRIM(rlab(i)), "=", dft_control%rtp_control%field(i), i=1, 3)
857 : END IF
858 :
859 15 : IF (new_file) THEN
860 3 : IF (dft_control%apply_efield_field) THEN
861 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,3(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z"
862 2 : ELSE IF (dft_control%apply_vector_potential) THEN
863 0 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,6(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z", &
864 0 : " Vec. Pot. X", " Vec. Pot. Y", " Vec. Pot. Z"
865 : END IF
866 : END IF
867 :
868 15 : field = 0.0_dp
869 15 : IF (dft_control%apply_efield_field) THEN
870 4 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
871 4 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
872 8 : field(1), field(2), field(3)
873 : ! DO i=1,3
874 : ! IF (ABS(field(i))< 10E-10) field(i) = 0.0_dp
875 : ! END IF
876 11 : ELSE IF (dft_control%apply_vector_potential) THEN
877 9 : WRITE (UNIT=unit_nr, FMT="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
878 9 : dft_control%rtp_control%field(1), dft_control%rtp_control%field(2), dft_control%rtp_control%field(3), &
879 18 : dft_control%rtp_control%vec_pot(1), dft_control%rtp_control%vec_pot(2), dft_control%rtp_control%vec_pot(3)
880 : END IF
881 :
882 : END IF
883 :
884 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
885 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD")
886 :
887 30 : END SUBROUTINE print_field_applied
888 :
889 : ! **************************************************************************************************
890 : !> \brief Print the components of the total energy used in an RTP calculation
891 : !> \param qs_env ...
892 : !> \param dft_section ...
893 : !> \par History
894 : !> 2024-02 Created [ANB]
895 : ! **************************************************************************************************
896 4 : SUBROUTINE print_rtp_energy_components(qs_env, dft_section)
897 : TYPE(qs_environment_type), POINTER :: qs_env
898 : TYPE(section_vals_type), POINTER :: dft_section
899 :
900 : CHARACTER(LEN=default_path_length) :: filename
901 : INTEGER :: i_step, output_unit, unit_nr
902 : LOGICAL :: new_file
903 : TYPE(cp_logger_type), POINTER :: logger
904 : TYPE(dft_control_type), POINTER :: dft_control
905 : TYPE(qs_energy_type), POINTER :: energy
906 : TYPE(rt_prop_type), POINTER :: rtp
907 :
908 4 : NULLIFY (dft_control, energy, rtp)
909 :
910 4 : logger => cp_get_default_logger()
911 4 : output_unit = cp_logger_get_default_io_unit(logger)
912 :
913 4 : CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
914 4 : i_step = rtp%istep
915 :
916 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
917 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS", extension=".ener", &
918 4 : file_action="WRITE", is_new_file=new_file)
919 :
920 4 : IF (output_unit > 0) THEN
921 2 : IF (unit_nr /= output_unit) THEN
922 2 : INQUIRE (UNIT=unit_nr, NAME=filename)
923 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
924 2 : "ENERGY_CONSTITUENTS", "Total Energy constituents written to file:", &
925 4 : TRIM(filename)
926 : ELSE
927 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ENERGY_CONSTITUENTS"
928 : END IF
929 :
930 2 : IF (new_file) THEN
931 : ! NOTE that these are not all terms contributing to the total energy for RTP, only a selection of those
932 : ! most significant / impactful. Therefore the printed components likely will not add up to the total energy.
933 1 : WRITE (UNIT=unit_nr, FMT='("#",5X,A,8X,A,10(6X,A))') "Step Nr.", "Time[fs]", &
934 1 : "Total ener.[a.u.]", "core[a.u.] ", " overlap [a.u.]", "hartree[a.u.]", " exc. [a.u.] ", &
935 2 : " hartree 1c[a.u.]", "exc 1c[a.u.] ", "exc admm[a.u.]", "exc 1c admm[a.u.]", "efield LG"
936 :
937 : END IF
938 : WRITE (UNIT=unit_nr, FMT="(I10,F20.6,10(F20.9))") &
939 2 : qs_env%sim_step, qs_env%sim_time*femtoseconds, &
940 2 : energy%total, energy%core, energy%core_overlap, energy%hartree, energy%exc, &
941 4 : energy%hartree_1c, energy%exc1, energy%exc_aux_fit, energy%exc1_aux_fit, energy%efield_core
942 :
943 : END IF
944 :
945 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
946 4 : "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS")
947 :
948 4 : END SUBROUTINE print_rtp_energy_components
949 :
950 : END MODULE rt_propagation_output
|
