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