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