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_projection_mo_utils, ONLY: compute_and_write_proj_mo
87 : USE rt_propagation_types, ONLY: get_rtp,&
88 : rt_prop_type
89 : USE rt_propagation_utils, ONLY: calculate_P_imaginary,&
90 : write_rtp_mo_cubes,&
91 : write_rtp_mos_to_output_unit
92 : #include "../base/base_uses.f90"
93 :
94 : IMPLICIT NONE
95 :
96 : PRIVATE
97 :
98 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_output'
99 :
100 : PUBLIC :: rt_prop_output, &
101 : rt_convergence, &
102 : rt_convergence_density, &
103 : report_density_occupation
104 :
105 : CONTAINS
106 :
107 : ! **************************************************************************************************
108 : !> \brief ...
109 : !> \param qs_env ...
110 : !> \param run_type ...
111 : !> \param delta_iter ...
112 : !> \param used_time ...
113 : ! **************************************************************************************************
114 4652 : SUBROUTINE rt_prop_output(qs_env, run_type, delta_iter, used_time)
115 : TYPE(qs_environment_type), POINTER :: qs_env
116 : INTEGER, INTENT(in) :: run_type
117 : REAL(dp), INTENT(in), OPTIONAL :: delta_iter, used_time
118 :
119 : INTEGER :: n_electrons, n_proj, nspin, output_unit, &
120 : spin
121 : REAL(dp) :: orthonormality, tot_rho_r
122 2326 : REAL(KIND=dp), DIMENSION(:), POINTER :: qs_tot_rho_r
123 2326 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
124 2326 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
125 : TYPE(cp_logger_type), POINTER :: logger
126 2326 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, P_im, rho_new
127 : TYPE(dft_control_type), POINTER :: dft_control
128 2326 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
129 2326 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
130 : TYPE(qs_rho_type), POINTER :: rho
131 : TYPE(rt_prop_type), POINTER :: rtp
132 : TYPE(section_vals_type), POINTER :: dft_section, input, rtp_section
133 :
134 2326 : NULLIFY (logger, dft_control)
135 :
136 4652 : logger => cp_get_default_logger()
137 : CALL get_qs_env(qs_env, &
138 : rtp=rtp, &
139 : matrix_s=matrix_s, &
140 : input=input, &
141 : rho=rho, &
142 : particle_set=particle_set, &
143 : atomic_kind_set=atomic_kind_set, &
144 : qs_kind_set=qs_kind_set, &
145 2326 : dft_control=dft_control)
146 :
147 2326 : rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
148 :
149 2326 : CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
150 2326 : n_electrons = n_electrons - dft_control%charge
151 :
152 2326 : CALL qs_rho_get(rho_struct=rho, tot_rho_r=qs_tot_rho_r)
153 :
154 2326 : tot_rho_r = accurate_sum(qs_tot_rho_r)
155 :
156 : output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
157 2326 : extension=".scfLog")
158 :
159 2326 : IF (output_unit > 0) THEN
160 : WRITE (output_unit, FMT="(/,(T3,A,T40,I5))") &
161 1163 : "Information at iteration step:", rtp%iter
162 : WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
163 1163 : "Total electronic density (r-space): ", &
164 1163 : tot_rho_r, &
165 : tot_rho_r + &
166 2326 : REAL(n_electrons, dp)
167 : WRITE (UNIT=output_unit, FMT="((T3,A,T59,F22.14))") &
168 1163 : "Total energy:", rtp%energy_new
169 1163 : IF (run_type == ehrenfest) &
170 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
171 628 : "Energy difference to previous iteration step:", rtp%energy_new - rtp%energy_old
172 1163 : IF (run_type == real_time_propagation) &
173 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F20.14))") &
174 535 : "Energy difference to initial state:", rtp%energy_new - rtp%energy_old
175 1163 : IF (PRESENT(delta_iter)) &
176 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,E20.6))") &
177 1163 : "Convergence:", delta_iter
178 1163 : IF (rtp%converged) THEN
179 294 : IF (run_type == real_time_propagation) &
180 : WRITE (UNIT=output_unit, FMT="((T3,A,T61,F12.2))") &
181 166 : "Time needed for propagation:", used_time
182 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,3X,F16.14))") &
183 294 : "CONVERGENCE REACHED", rtp%energy_new - rtp%energy_old
184 : END IF
185 : END IF
186 :
187 2326 : IF (rtp%converged) THEN
188 588 : IF (.NOT. rtp%linear_scaling) THEN
189 406 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
190 : CALL rt_calculate_orthonormality(orthonormality, &
191 406 : mos_new, matrix_s(1)%matrix)
192 406 : IF (output_unit > 0) &
193 : WRITE (output_unit, FMT="(/,(T3,A,T60,F20.10))") &
194 203 : "Max deviation from orthonormalization:", orthonormality
195 : END IF
196 : END IF
197 :
198 2326 : IF (output_unit > 0) &
199 1163 : CALL m_flush(output_unit)
200 : CALL cp_print_key_finished_output(output_unit, logger, rtp_section, &
201 2326 : "PRINT%PROGRAM_RUN_INFO")
202 :
203 2326 : IF (rtp%converged) THEN
204 588 : dft_section => section_vals_get_subs_vals(input, "DFT")
205 588 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
206 : dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) &
207 30 : CALL print_field_applied(qs_env, dft_section)
208 588 : CALL make_moment(qs_env)
209 588 : IF (.NOT. dft_control%qs_control%dftb) THEN
210 468 : CALL write_available_results(qs_env=qs_env, rtp=rtp)
211 : END IF
212 588 : IF (rtp%linear_scaling) THEN
213 182 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
214 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
215 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART"), cp_p_file)) THEN
216 96 : CALL write_rt_p_to_restart(rho_new, .FALSE.)
217 : END IF
218 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
219 : dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY"), cp_p_file)) THEN
220 2 : CALL write_rt_p_to_restart(rho_new, .TRUE.)
221 : END IF
222 182 : IF (.NOT. dft_control%qs_control%dftb) THEN
223 : !Not sure if these things could also work with dftb or not
224 182 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
225 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
226 64 : DO spin = 1, SIZE(rho_new)/2
227 64 : CALL rt_current(qs_env, rho_new(2*spin)%matrix, dft_section, spin, SIZE(rho_new)/2)
228 : END DO
229 : END IF
230 : END IF
231 : ELSE
232 406 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
233 406 : IF (.NOT. dft_control%qs_control%dftb) THEN
234 286 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
235 : dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
236 12 : NULLIFY (P_im)
237 12 : nspin = SIZE(mos_new)/2
238 12 : CALL dbcsr_allocate_matrix_set(P_im, nspin)
239 24 : DO spin = 1, nspin
240 12 : CALL dbcsr_init_p(P_im(spin)%matrix)
241 24 : CALL dbcsr_create(P_im(spin)%matrix, template=matrix_s(1)%matrix, matrix_type="N")
242 : END DO
243 12 : CALL calculate_P_imaginary(qs_env, rtp, P_im)
244 24 : DO spin = 1, nspin
245 24 : CALL rt_current(qs_env, P_im(spin)%matrix, dft_section, spin, nspin)
246 : END DO
247 12 : CALL dbcsr_deallocate_matrix_set(P_im)
248 : END IF
249 286 : IF (dft_control%rtp_control%is_proj_mo) THEN
250 44 : DO n_proj = 1, SIZE(dft_control%rtp_control%proj_mo_list)
251 : CALL compute_and_write_proj_mo(qs_env, mos_new, &
252 44 : dft_control%rtp_control%proj_mo_list(n_proj)%proj_mo, n_proj)
253 : END DO
254 : END IF
255 : END IF
256 : CALL write_rt_mos_to_restart(qs_env%mos, mos_new, particle_set, &
257 406 : dft_section, qs_kind_set)
258 : END IF
259 : END IF
260 :
261 2326 : rtp%energy_old = rtp%energy_new
262 :
263 2326 : IF (.NOT. rtp%converged .AND. rtp%iter >= dft_control%rtp_control%max_iter) &
264 : CALL cp_abort(__LOCATION__, "EMD did not converge, either increase MAX_ITER "// &
265 0 : "or use a smaller TIMESTEP")
266 :
267 2326 : END SUBROUTINE rt_prop_output
268 :
269 : ! **************************************************************************************************
270 : !> \brief computes the effective orthonormality of a set of mos given an s-matrix
271 : !> orthonormality is the max deviation from unity of the C^T S C
272 : !> \param orthonormality ...
273 : !> \param mos_new ...
274 : !> \param matrix_s ...
275 : !> \author Florian Schiffmann (02.09)
276 : ! **************************************************************************************************
277 406 : SUBROUTINE rt_calculate_orthonormality(orthonormality, mos_new, matrix_s)
278 : REAL(KIND=dp), INTENT(out) :: orthonormality
279 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
280 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
281 :
282 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_calculate_orthonormality'
283 :
284 : INTEGER :: handle, i, im, ispin, j, k, n, &
285 : ncol_local, nrow_local, nspin, re
286 406 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
287 : REAL(KIND=dp) :: alpha, max_alpha, max_beta
288 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
289 : TYPE(cp_fm_type) :: overlap_re, svec_im, svec_re
290 :
291 406 : NULLIFY (tmp_fm_struct)
292 :
293 406 : CALL timeset(routineN, handle)
294 :
295 406 : nspin = SIZE(mos_new)/2
296 406 : max_alpha = 0.0_dp
297 406 : max_beta = 0.0_dp
298 898 : DO ispin = 1, nspin
299 492 : re = ispin*2 - 1
300 492 : im = ispin*2
301 : ! get S*C
302 492 : CALL cp_fm_create(svec_re, mos_new(im)%matrix_struct)
303 492 : CALL cp_fm_create(svec_im, mos_new(im)%matrix_struct)
304 : CALL cp_fm_get_info(mos_new(im), &
305 492 : nrow_global=n, ncol_global=k)
306 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(re), &
307 492 : svec_re, k)
308 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(im), &
309 492 : svec_im, k)
310 :
311 : ! get C^T (S*C)
312 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
313 : para_env=mos_new(re)%matrix_struct%para_env, &
314 492 : context=mos_new(re)%matrix_struct%context)
315 492 : CALL cp_fm_create(overlap_re, tmp_fm_struct)
316 :
317 492 : CALL cp_fm_struct_release(tmp_fm_struct)
318 :
319 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(re), &
320 492 : svec_re, 0.0_dp, overlap_re)
321 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(im), &
322 492 : svec_im, 1.0_dp, overlap_re)
323 :
324 492 : CALL cp_fm_release(svec_re)
325 492 : CALL cp_fm_release(svec_im)
326 :
327 : CALL cp_fm_get_info(overlap_re, nrow_local=nrow_local, ncol_local=ncol_local, &
328 492 : row_indices=row_indices, col_indices=col_indices)
329 1279 : DO i = 1, nrow_local
330 4470 : DO j = 1, ncol_local
331 3191 : alpha = overlap_re%local_data(i, j)
332 3191 : IF (row_indices(i) .EQ. col_indices(j)) alpha = alpha - 1.0_dp
333 3978 : max_alpha = MAX(max_alpha, ABS(alpha))
334 : END DO
335 : END DO
336 1882 : CALL cp_fm_release(overlap_re)
337 : END DO
338 406 : CALL mos_new(1)%matrix_struct%para_env%max(max_alpha)
339 406 : CALL mos_new(1)%matrix_struct%para_env%max(max_beta)
340 406 : orthonormality = max_alpha
341 :
342 406 : CALL timestop(handle)
343 :
344 406 : END SUBROUTINE rt_calculate_orthonormality
345 :
346 : ! **************************************************************************************************
347 : !> \brief computes the convergence criterion for RTP and EMD
348 : !> \param rtp ...
349 : !> \param matrix_s Overlap matrix without the derivatives
350 : !> \param delta_mos ...
351 : !> \param delta_eps ...
352 : !> \author Florian Schiffmann (02.09)
353 : ! **************************************************************************************************
354 :
355 1510 : SUBROUTINE rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
356 : TYPE(rt_prop_type), POINTER :: rtp
357 : TYPE(dbcsr_type), POINTER :: matrix_s
358 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: delta_mos
359 : REAL(dp), INTENT(out) :: delta_eps
360 :
361 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence'
362 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
363 :
364 : INTEGER :: handle, i, icol, im, ispin, j, lcol, &
365 : lrow, nao, newdim, nmo, nspin, re
366 : LOGICAL :: double_col, double_row
367 : REAL(KIND=dp) :: alpha, max_alpha
368 : TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1, tmp_fm_struct
369 : TYPE(cp_fm_type) :: work, work1, work2
370 1510 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
371 :
372 1510 : NULLIFY (tmp_fm_struct)
373 :
374 1510 : CALL timeset(routineN, handle)
375 :
376 1510 : CALL get_rtp(rtp=rtp, mos_new=mos_new)
377 :
378 1510 : nspin = SIZE(delta_mos)/2
379 1510 : max_alpha = 0.0_dp
380 :
381 5190 : DO i = 1, SIZE(mos_new)
382 5190 : CALL cp_fm_scale_and_add(-one, delta_mos(i), one, mos_new(i))
383 : END DO
384 :
385 3350 : DO ispin = 1, nspin
386 1840 : re = ispin*2 - 1
387 1840 : im = ispin*2
388 :
389 1840 : double_col = .TRUE.
390 1840 : double_row = .FALSE.
391 : CALL cp_fm_struct_double(newstruct, &
392 : delta_mos(re)%matrix_struct, &
393 : delta_mos(re)%matrix_struct%context, &
394 : double_col, &
395 1840 : double_row)
396 :
397 1840 : CALL cp_fm_create(work, matrix_struct=newstruct)
398 1840 : CALL cp_fm_create(work1, matrix_struct=newstruct)
399 :
400 : CALL cp_fm_get_info(delta_mos(re), ncol_local=lcol, ncol_global=nmo, &
401 1840 : nrow_global=nao)
402 1840 : CALL cp_fm_get_info(work, ncol_global=newdim)
403 :
404 1840 : CALL cp_fm_set_all(work, zero, zero)
405 :
406 7440 : DO icol = 1, lcol
407 42312 : work%local_data(:, icol) = delta_mos(re)%local_data(:, icol)
408 44152 : work%local_data(:, icol + lcol) = delta_mos(im)%local_data(:, icol)
409 : END DO
410 :
411 1840 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, work, work1, ncol=newdim)
412 :
413 1840 : CALL cp_fm_release(work)
414 :
415 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, &
416 : para_env=delta_mos(re)%matrix_struct%para_env, &
417 1840 : context=delta_mos(re)%matrix_struct%context)
418 : CALL cp_fm_struct_double(newstruct1, &
419 : tmp_fm_struct, &
420 : delta_mos(re)%matrix_struct%context, &
421 : double_col, &
422 1840 : double_row)
423 :
424 1840 : CALL cp_fm_create(work, matrix_struct=newstruct1)
425 1840 : CALL cp_fm_create(work2, matrix_struct=newstruct1)
426 :
427 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(re), &
428 1840 : work1, zero, work)
429 :
430 : CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(im), &
431 1840 : work1, zero, work2)
432 :
433 1840 : CALL cp_fm_get_info(work, nrow_local=lrow)
434 4640 : DO i = 1, lrow
435 15820 : DO j = 1, lcol
436 : alpha = SQRT((work%local_data(i, j) + work2%local_data(i, j + lcol))**2 + &
437 11180 : (work%local_data(i, j + lcol) - work2%local_data(i, j))**2)
438 13980 : max_alpha = MAX(max_alpha, ABS(alpha))
439 : END DO
440 : END DO
441 :
442 1840 : CALL cp_fm_release(work)
443 1840 : CALL cp_fm_release(work1)
444 1840 : CALL cp_fm_release(work2)
445 1840 : CALL cp_fm_struct_release(tmp_fm_struct)
446 1840 : CALL cp_fm_struct_release(newstruct)
447 7030 : CALL cp_fm_struct_release(newstruct1)
448 :
449 : END DO
450 :
451 1510 : CALL delta_mos(1)%matrix_struct%para_env%max(max_alpha)
452 1510 : delta_eps = SQRT(max_alpha)
453 :
454 1510 : CALL timestop(handle)
455 :
456 1510 : END SUBROUTINE rt_convergence
457 :
458 : ! **************************************************************************************************
459 : !> \brief computes the convergence criterion for RTP and EMD based on the density matrix
460 : !> \param rtp ...
461 : !> \param delta_P ...
462 : !> \param delta_eps ...
463 : !> \author Samuel Andermatt (02.14)
464 : ! **************************************************************************************************
465 :
466 1632 : SUBROUTINE rt_convergence_density(rtp, delta_P, delta_eps)
467 :
468 : TYPE(rt_prop_type), POINTER :: rtp
469 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P
470 : REAL(dp), INTENT(out) :: delta_eps
471 :
472 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_convergence_density'
473 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
474 :
475 : INTEGER :: col_atom, group_handle, handle, i, &
476 : 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 712140 : block_values = block_values*block_values
498 : END DO
499 2940 : 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 365095 : alpha = MAXVAL(block_values)
515 5213 : IF (alpha > max_alpha) max_alpha = alpha
516 : END DO
517 1878 : CALL dbcsr_iterator_stop(iter)
518 : END DO
519 816 : CALL dbcsr_get_info(delta_P(1)%matrix, group=group_handle)
520 816 : CALL group%set_handle(group_handle)
521 816 : CALL group%max(max_alpha)
522 816 : delta_eps = SQRT(max_alpha)
523 816 : CALL dbcsr_deallocate_matrix(tmp)
524 816 : CALL timestop(handle)
525 :
526 816 : END SUBROUTINE rt_convergence_density
527 :
528 : ! **************************************************************************************************
529 : !> \brief interface to qs_moments. Does only work for nonperiodic dipole
530 : !> \param qs_env ...
531 : !> \author Florian Schiffmann (02.09)
532 : ! **************************************************************************************************
533 :
534 588 : SUBROUTINE make_moment(qs_env)
535 :
536 : TYPE(qs_environment_type), POINTER :: qs_env
537 :
538 : CHARACTER(len=*), PARAMETER :: routineN = 'make_moment'
539 :
540 : INTEGER :: handle, output_unit
541 : TYPE(cp_logger_type), POINTER :: logger
542 : TYPE(dft_control_type), POINTER :: dft_control
543 :
544 588 : CALL timeset(routineN, handle)
545 :
546 588 : NULLIFY (dft_control)
547 :
548 588 : logger => cp_get_default_logger()
549 588 : output_unit = cp_logger_get_default_io_unit(logger)
550 588 : CALL get_qs_env(qs_env, dft_control=dft_control)
551 588 : IF (dft_control%qs_control%dftb) THEN
552 120 : CALL scf_post_calculation_tb(qs_env, "DFTB", .FALSE.)
553 468 : ELSE IF (dft_control%qs_control%xtb) THEN
554 60 : CALL scf_post_calculation_tb(qs_env, "xTB", .FALSE.)
555 : ELSE
556 408 : CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
557 : END IF
558 588 : CALL timestop(handle)
559 :
560 588 : END SUBROUTINE make_moment
561 :
562 : ! **************************************************************************************************
563 : !> \brief Reports the sparsity pattern of the complex density matrix
564 : !> \param filter_eps ...
565 : !> \param rho ...
566 : !> \author Samuel Andermatt (09.14)
567 : ! **************************************************************************************************
568 :
569 182 : SUBROUTINE report_density_occupation(filter_eps, rho)
570 :
571 : REAL(KIND=dp) :: filter_eps
572 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho
573 :
574 : CHARACTER(len=*), PARAMETER :: routineN = 'report_density_occupation'
575 :
576 : INTEGER :: handle, i, im, ispin, re, unit_nr
577 : REAL(KIND=dp) :: eps, occ
578 : TYPE(cp_logger_type), POINTER :: logger
579 182 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmp
580 :
581 182 : CALL timeset(routineN, handle)
582 :
583 182 : logger => cp_get_default_logger()
584 182 : unit_nr = cp_logger_get_default_io_unit(logger)
585 182 : NULLIFY (tmp)
586 182 : CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
587 686 : DO i = 1, SIZE(rho)
588 504 : CALL dbcsr_init_p(tmp(i)%matrix)
589 504 : CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
590 686 : CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
591 : END DO
592 434 : DO ispin = 1, SIZE(rho)/2
593 252 : re = 2*ispin - 1
594 252 : im = 2*ispin
595 252 : eps = MAX(filter_eps, 1.0E-11_dp)
596 2338 : DO WHILE (eps < 1.1_dp)
597 2086 : CALL dbcsr_filter(tmp(re)%matrix, eps)
598 2086 : occ = dbcsr_get_occupation(tmp(re)%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, " real: ", occ
601 2086 : eps = eps*10
602 : END DO
603 252 : eps = MAX(filter_eps, 1.0E-11_dp)
604 2520 : DO WHILE (eps < 1.1_dp)
605 2086 : CALL dbcsr_filter(tmp(im)%matrix, eps)
606 2086 : occ = dbcsr_get_occupation(tmp(im)%matrix)
607 3129 : IF (unit_nr > 0) WRITE (unit_nr, FMT="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
608 2086 : ispin, " eps ", eps, " imag: ", occ
609 2086 : eps = eps*10.0_dp
610 : END DO
611 : END DO
612 182 : CALL dbcsr_deallocate_matrix_set(tmp)
613 182 : CALL timestop(handle)
614 :
615 182 : END SUBROUTINE report_density_occupation
616 :
617 : ! **************************************************************************************************
618 : !> \brief Writes the density matrix and the atomic positions to a restart file
619 : !> \param rho_new ...
620 : !> \param history ...
621 : !> \author Samuel Andermatt (09.14)
622 : ! **************************************************************************************************
623 :
624 98 : SUBROUTINE write_rt_p_to_restart(rho_new, history)
625 :
626 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
627 : LOGICAL :: history
628 :
629 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_p_to_restart'
630 :
631 : CHARACTER(LEN=default_path_length) :: file_name, project_name
632 : INTEGER :: handle, im, ispin, re, unit_nr
633 : REAL(KIND=dp) :: cs_pos
634 : TYPE(cp_logger_type), POINTER :: logger
635 :
636 98 : CALL timeset(routineN, handle)
637 98 : logger => cp_get_default_logger()
638 98 : IF (logger%para_env%is_source()) THEN
639 49 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
640 : ELSE
641 : unit_nr = -1
642 : END IF
643 :
644 98 : project_name = logger%iter_info%project_name
645 234 : DO ispin = 1, SIZE(rho_new)/2
646 136 : re = 2*ispin - 1
647 136 : im = 2*ispin
648 136 : IF (history) THEN
649 : WRITE (file_name, '(A,I0,A)') &
650 2 : TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
651 : ELSE
652 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
653 : END IF
654 136 : cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.TRUE.)
655 136 : IF (unit_nr > 0) THEN
656 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
657 : END IF
658 136 : CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
659 136 : IF (history) THEN
660 : WRITE (file_name, '(A,I0,A)') &
661 2 : TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//TRIM(cp_iter_string(logger%iter_info))//"_RESTART.dm"
662 : ELSE
663 134 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
664 : END IF
665 136 : cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.TRUE.)
666 136 : IF (unit_nr > 0) THEN
667 68 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
668 : END IF
669 234 : CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
670 : END DO
671 :
672 98 : CALL timestop(handle)
673 :
674 98 : END SUBROUTINE write_rt_p_to_restart
675 :
676 : ! **************************************************************************************************
677 : !> \brief Collocation of the current and printing of it in a cube file
678 : !> \param qs_env ...
679 : !> \param P_im ...
680 : !> \param dft_section ...
681 : !> \param spin ...
682 : !> \param nspin ...
683 : !> \author Samuel Andermatt (06.15)
684 : ! **************************************************************************************************
685 48 : SUBROUTINE rt_current(qs_env, P_im, dft_section, spin, nspin)
686 : TYPE(qs_environment_type), POINTER :: qs_env
687 : TYPE(dbcsr_type), POINTER :: P_im
688 : TYPE(section_vals_type), POINTER :: dft_section
689 : INTEGER :: spin, nspin
690 :
691 : CHARACTER(len=*), PARAMETER :: routineN = 'rt_current'
692 :
693 : CHARACTER(len=1) :: char_spin
694 : CHARACTER(len=14) :: ext
695 : CHARACTER(len=2) :: sdir
696 : INTEGER :: dir, handle, print_unit
697 48 : INTEGER, DIMENSION(:), POINTER :: stride(:)
698 : LOGICAL :: mpi_io
699 : TYPE(cp_logger_type), POINTER :: logger
700 : TYPE(current_env_type) :: current_env
701 : TYPE(dbcsr_type), POINTER :: tmp, zero
702 : TYPE(particle_list_type), POINTER :: particles
703 : TYPE(pw_env_type), POINTER :: pw_env
704 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
705 : TYPE(pw_type) :: gs, 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 pw_pool_create_pw(auxbas_pw_pool, rs, use_data=REALDATA3D, in_space=REALSPACE)
728 48 : CALL pw_pool_create_pw(auxbas_pw_pool, gs, use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
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 pw_pool_give_back_pw(auxbas_pw_pool, rs)
766 48 : CALL pw_pool_give_back_pw(auxbas_pw_pool, 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 468 : 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 468 : CALL timeset(routineN, handle)
797 :
798 468 : CALL get_qs_env(qs_env, scf_env=scf_env)
799 468 : IF (rtp%linear_scaling) THEN
800 182 : CALL write_mo_free_results(qs_env)
801 : ELSE
802 286 : CALL write_mo_free_results(qs_env)
803 286 : CALL write_mo_dependent_results(qs_env, scf_env)
804 : ! Time-dependent MO print
805 286 : CALL write_rtp_mos_to_output_unit(qs_env, rtp)
806 286 : CALL write_rtp_mo_cubes(qs_env, rtp)
807 : END IF
808 :
809 468 : CALL timestop(handle)
810 :
811 468 : 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, output_unit, unit_nr
828 : REAL(kind=dp) :: field(3), to_write(3)
829 : TYPE(cp_logger_type), POINTER :: logger
830 : TYPE(dft_control_type), POINTER :: dft_control
831 :
832 30 : NULLIFY (dft_control)
833 :
834 30 : logger => cp_get_default_logger()
835 30 : output_unit = cp_logger_get_default_io_unit(logger)
836 :
837 30 : CALL get_qs_env(qs_env, dft_control=dft_control)
838 :
839 : unit_nr = cp_print_key_unit_nr(logger, dft_section, &
840 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat")
841 :
842 30 : IF (output_unit > 0) THEN
843 15 : IF (unit_nr /= output_unit) THEN
844 15 : INQUIRE (UNIT=unit_nr, NAME=filename)
845 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
846 15 : "FIELD", "The field applied is written to the file:", &
847 30 : TRIM(filename)
848 : WRITE (UNIT=unit_nr, FMT="(/,(T2,A,T40,I6))") &
849 15 : "Real time propagation step:", qs_env%sim_step
850 : ELSE
851 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED"
852 : END IF
853 :
854 60 : rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
855 :
856 15 : IF (dft_control%apply_efield_field) THEN
857 4 : WRITE (unit_nr, "(T3,A)") "Electric Field (LG) in atomic units:"
858 4 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
859 4 : to_write = field
860 11 : ELSE IF (dft_control%apply_vector_potential) THEN
861 9 : WRITE (unit_nr, "(T3,A)") "Vector potential (VG) in atomic units:"
862 36 : to_write = dft_control%rtp_control%vec_pot
863 : ELSE
864 2 : WRITE (unit_nr, "(T3,A)") "No electric field applied"
865 2 : to_write = 0._dp
866 : END IF
867 :
868 : WRITE (unit_nr, "(T5,3(A,A,E16.8,1X))") &
869 60 : (TRIM(rlab(i)), "=", to_write(i), i=1, 3)
870 : END IF
871 :
872 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
873 30 : "REAL_TIME_PROPAGATION%PRINT%FIELD")
874 :
875 30 : END SUBROUTINE print_field_applied
876 :
877 : END MODULE rt_propagation_output
|
