(git:419edc0)
Loading...
Searching...
No Matches
rt_propagation_output.F
Go to the documentation of this file.
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
16 USE cp_dbcsr_api, ONLY: &
17 dbcsr_add, dbcsr_binary_write, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
22 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
34 USE cp_fm_types, ONLY: cp_fm_create,&
45 cp_p_file,&
50 USE efield_utils, ONLY: make_field
51 USE input_constants, ONLY: ehrenfest,&
56 USE kahan_sum, ONLY: accurate_sum
57 USE kinds, ONLY: default_path_length,&
58 dp
59 USE machine, ONLY: m_flush
65 USE physcon, ONLY: femtoseconds
66 USE pw_env_types, ONLY: pw_env_get,&
68 USE pw_methods, ONLY: pw_zero
70 USE pw_types, ONLY: pw_c1d_gs_type,&
82 USE qs_rho_types, ONLY: qs_rho_get,&
92 USE rt_propagation_types, ONLY: get_rtp,&
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, &
108
109CONTAINS
110
111! **************************************************************************************************
112!> \brief ...
113!> \param qs_env ...
114!> \param run_type ...
115!> \param delta_iter ...
116!> \param used_time ...
117! **************************************************************************************************
118 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
126 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
127 LOGICAL :: new_file
128 REAL(dp) :: orthonormality, strace, tot_rho_r, trace
129 REAL(kind=dp), DIMENSION(:), POINTER :: qs_tot_rho_r
130 REAL(kind=dp), DIMENSION(:, :), POINTER :: j_int
131 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
132 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
133 TYPE(cp_logger_type), POINTER :: logger
134 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
135 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 POINTER :: sab_all, sab_orb
140 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
141 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 NULLIFY (logger, dft_control)
147
148 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 dbcsr_dist=dbcsr_dist)
159
160 rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
161
162 CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
163 n_electrons = n_electrons - dft_control%charge
164
165 CALL qs_rho_get(rho_struct=rho, tot_rho_r=qs_tot_rho_r)
166
167 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 extension=".scfLog")
171
172 IF (output_unit > 0) THEN
173 WRITE (output_unit, fmt="(/,(T3,A,T40,I5))") &
174 "Information at iteration step:", rtp%iter
175 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
176 "Total electronic density (r-space): ", &
177 tot_rho_r, &
178 tot_rho_r + &
179 REAL(n_electrons, dp)
180 WRITE (unit=output_unit, fmt="((T3,A,T59,F22.14))") &
181 "Total energy:", rtp%energy_new
182 IF (run_type == ehrenfest) &
183 WRITE (unit=output_unit, fmt="((T3,A,T61,F20.14))") &
184 "Energy difference to previous iteration step:", rtp%energy_new - rtp%energy_old
185 IF (run_type == real_time_propagation) &
186 WRITE (unit=output_unit, fmt="((T3,A,T61,F20.14))") &
187 "Energy difference to initial state:", rtp%energy_new - rtp%energy_old
188 IF (PRESENT(delta_iter)) &
189 WRITE (unit=output_unit, fmt="((T3,A,T61,E20.6))") &
190 "Convergence:", delta_iter
191 IF (rtp%converged) THEN
192 IF (run_type == real_time_propagation) &
193 WRITE (unit=output_unit, fmt="((T3,A,T61,F12.2))") &
194 "Time needed for propagation:", used_time
195 WRITE (unit=output_unit, fmt="(/,(T3,A,3X,F16.14))") &
196 "CONVERGENCE REACHED", rtp%energy_new - rtp%energy_old
197 END IF
198 END IF
199
200 IF (rtp%converged) THEN
201 IF (.NOT. rtp%linear_scaling) THEN
202 CALL get_rtp(rtp=rtp, mos_new=mos_new)
203 CALL rt_calculate_orthonormality(orthonormality, &
204 mos_new, matrix_s(1)%matrix)
205 IF (output_unit > 0) &
206 WRITE (output_unit, fmt="(/,(T3,A,T60,F20.10))") &
207 "Max deviation from orthonormalization:", orthonormality
208 END IF
209 END IF
210
211 IF (output_unit > 0) &
212 CALL m_flush(output_unit)
213 CALL cp_print_key_finished_output(output_unit, logger, rtp_section, &
214 "PRINT%PROGRAM_RUN_INFO")
215
216 IF (rtp%converged) THEN
217 dft_section => section_vals_get_subs_vals(input, "DFT")
218 IF (btest(cp_print_key_should_output(logger%iter_info, &
219 dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) &
220 CALL print_field_applied(qs_env, dft_section)
221 CALL make_moment(qs_env)
222 IF (btest(cp_print_key_should_output(logger%iter_info, &
223 dft_section, "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS"), cp_p_file)) THEN
224 CALL print_rtp_energy_components(qs_env, dft_section)
225 END IF
226 IF (.NOT. dft_control%qs_control%dftb) THEN
227 CALL write_available_results(qs_env=qs_env, rtp=rtp)
228 END IF
229
230 IF (rtp%linear_scaling) THEN
231 CALL get_rtp(rtp=rtp, rho_new=rho_new)
232 IF (btest(cp_print_key_should_output(logger%iter_info, &
233 dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART"), cp_p_file)) THEN
234 CALL write_rt_p_to_restart(rho_new, .false.)
235 END IF
236 IF (btest(cp_print_key_should_output(logger%iter_info, &
237 dft_section, "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY"), cp_p_file)) THEN
238 CALL write_rt_p_to_restart(rho_new, .true.)
239 END IF
240 IF (.NOT. dft_control%qs_control%dftb) THEN
241 !Not sure if these things could also work with dftb or not
242 IF (btest(cp_print_key_should_output(logger%iter_info, &
243 dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
244 DO spin = 1, SIZE(rho_new)/2
245 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 CALL get_rtp(rtp=rtp, mos_new=mos_new)
251 IF (.NOT. dft_control%qs_control%dftb .AND. .NOT. dft_control%qs_control%xtb) THEN
252 IF (rtp%track_imag_density) THEN
253 NULLIFY (p_im, p_xyz)
254 CALL dbcsr_allocate_matrix_set(p_xyz, 3)
255
256! Linear momentum operator
257! prepare for allocation
258 natom = SIZE(particle_set, 1)
259 ALLOCATE (first_sgf(natom))
260 ALLOCATE (last_sgf(natom))
261 CALL get_particle_set(particle_set, qs_kind_set, &
262 first_sgf=first_sgf, &
263 last_sgf=last_sgf)
264 ALLOCATE (row_blk_sizes(natom))
265 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
266 DEALLOCATE (first_sgf)
267 DEALLOCATE (last_sgf)
268
269 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 mutable_work=.true.)
275 CALL cp_dbcsr_alloc_block_from_nbl(p_xyz(1)%matrix, sab_orb)
276 CALL dbcsr_set(p_xyz(1)%matrix, 0.0_dp)
277 DO i = 2, 3
278 ALLOCATE (p_xyz(i)%matrix)
279 CALL dbcsr_copy(p_xyz(i)%matrix, p_xyz(1)%matrix, "p_xyz-"//trim(adjustl(cp_to_string(i))))
280 CALL dbcsr_set(p_xyz(i)%matrix, 0.0_dp)
281 END DO
282 CALL build_lin_mom_matrix(qs_env, p_xyz)
283 DEALLOCATE (row_blk_sizes)
284
285 nspin = SIZE(mos_new)/2
286 CALL qs_rho_get(rho, rho_ao_im=p_im)
287 ALLOCATE (j_int(nspin, 3))
288 j_int = 0.0_dp
289
290 NULLIFY (tmp_ao)
291 CALL dbcsr_init_p(tmp_ao)
292 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
293 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
294 CALL dbcsr_set(tmp_ao, 0.0_dp)
295
296 DO i = 1, 3
297 strace = 0.0_dp
298 DO spin = 1, nspin
299 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 0.0_dp, tmp_ao)
302 CALL dbcsr_trace(tmp_ao, trace)
303 strace = strace + trace
304 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 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 output_unit = cp_logger_get_default_io_unit(logger)
313 unit_nr = cp_print_key_unit_nr(logger, dft_section, &
314 "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT", extension=".dat", is_new_file=new_file)
315 IF (output_unit > 0) THEN
316 IF (nspin == 2) THEN
317 WRITE (unit=unit_nr, fmt="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
318 j_int(1, 1:3), j_int(2, 1:3)
319 ELSE
320 WRITE (unit=unit_nr, fmt="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
321 j_int(1, 1:3)
322 END IF
323 END IF
324 CALL cp_print_key_finished_output(unit_nr, logger, dft_section, &
325 "REAL_TIME_PROPAGATION%PRINT%CURRENT_INT")
326 END IF
327 DEALLOCATE (j_int)
328
329 IF (btest(cp_print_key_should_output(logger%iter_info, &
330 dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT"), cp_p_file)) THEN
331 DO spin = 1, nspin
332 CALL rt_current(qs_env, p_im(spin)%matrix, dft_section, spin, nspin)
333 END DO
334 END IF
335 CALL dbcsr_deallocate_matrix(tmp_ao)
337 END IF
338
339! projection of molecular orbitals
340 IF (dft_control%rtp_control%is_proj_mo) THEN
341 DO n_proj = 1, SIZE(dft_control%rtp_control%proj_mo_list)
342 CALL compute_and_write_proj_mo(qs_env, mos_new, &
343 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 dft_section, qs_kind_set)
349 END IF
350 END IF
351
352 rtp%energy_old = rtp%energy_new
353
354 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 "or use a smaller TIMESTEP")
357
358 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 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 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 NULLIFY (tmp_fm_struct)
383
384 CALL timeset(routinen, handle)
385
386 nspin = SIZE(mos_new)/2
387 max_alpha = 0.0_dp
388 max_beta = 0.0_dp
389 DO ispin = 1, nspin
390 re = ispin*2 - 1
391 im = ispin*2
392 ! get S*C
393 CALL cp_fm_create(svec_re, mos_new(im)%matrix_struct)
394 CALL cp_fm_create(svec_im, mos_new(im)%matrix_struct)
395 CALL cp_fm_get_info(mos_new(im), &
396 nrow_global=n, ncol_global=k)
397 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(re), &
398 svec_re, k)
399 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mos_new(im), &
400 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 context=mos_new(re)%matrix_struct%context)
406 CALL cp_fm_create(overlap_re, tmp_fm_struct)
407
408 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 svec_re, 0.0_dp, overlap_re)
412 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mos_new(im), &
413 svec_im, 1.0_dp, overlap_re)
414
415 CALL cp_fm_release(svec_re)
416 CALL cp_fm_release(svec_im)
417
418 CALL cp_fm_get_info(overlap_re, nrow_local=nrow_local, ncol_local=ncol_local, &
419 row_indices=row_indices, col_indices=col_indices)
420 DO i = 1, nrow_local
421 DO j = 1, ncol_local
422 alpha = overlap_re%local_data(i, j)
423 IF (row_indices(i) .EQ. col_indices(j)) alpha = alpha - 1.0_dp
424 max_alpha = max(max_alpha, abs(alpha))
425 END DO
426 END DO
427 CALL cp_fm_release(overlap_re)
428 END DO
429 CALL mos_new(1)%matrix_struct%para_env%max(max_alpha)
430 CALL mos_new(1)%matrix_struct%para_env%max(max_beta)
431 orthonormality = max_alpha
432
433 CALL timestop(handle)
434
435 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 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 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
462
463 NULLIFY (tmp_fm_struct)
464
465 CALL timeset(routinen, handle)
466
467 CALL get_rtp(rtp=rtp, mos_new=mos_new)
468
469 nspin = SIZE(delta_mos)/2
470 max_alpha = 0.0_dp
471
472 DO i = 1, SIZE(mos_new)
473 CALL cp_fm_scale_and_add(-one, delta_mos(i), one, mos_new(i))
474 END DO
475
476 DO ispin = 1, nspin
477 re = ispin*2 - 1
478 im = ispin*2
479
480 double_col = .true.
481 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 double_row)
487
488 CALL cp_fm_create(work, matrix_struct=newstruct)
489 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 nrow_global=nao)
493 CALL cp_fm_get_info(work, ncol_global=newdim)
494
495 CALL cp_fm_set_all(work, zero, zero)
496
497 DO icol = 1, lcol
498 work%local_data(:, icol) = delta_mos(re)%local_data(:, icol)
499 work%local_data(:, icol + lcol) = delta_mos(im)%local_data(:, icol)
500 END DO
501
502 CALL cp_dbcsr_sm_fm_multiply(matrix_s, work, work1, ncol=newdim)
503
504 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 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 double_row)
514
515 CALL cp_fm_create(work, matrix_struct=newstruct1)
516 CALL cp_fm_create(work2, matrix_struct=newstruct1)
517
518 CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(re), &
519 work1, zero, work)
520
521 CALL parallel_gemm("T", "N", nmo, newdim, nao, one, delta_mos(im), &
522 work1, zero, work2)
523
524 CALL cp_fm_get_info(work, nrow_local=lrow)
525 DO i = 1, lrow
526 DO j = 1, lcol
527 alpha = sqrt((work%local_data(i, j) + work2%local_data(i, j + lcol))**2 + &
528 (work%local_data(i, j + lcol) - work2%local_data(i, j))**2)
529 max_alpha = max(max_alpha, abs(alpha))
530 END DO
531 END DO
532
533 CALL cp_fm_release(work)
534 CALL cp_fm_release(work1)
535 CALL cp_fm_release(work2)
536 CALL cp_fm_struct_release(tmp_fm_struct)
537 CALL cp_fm_struct_release(newstruct)
538 CALL cp_fm_struct_release(newstruct1)
539
540 END DO
541
542 CALL delta_mos(1)%matrix_struct%para_env%max(max_alpha)
543 delta_eps = sqrt(max_alpha)
544
545 CALL timestop(handle)
546
547 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 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 REAL(dp), DIMENSION(:, :), POINTER :: block_values
569 TYPE(dbcsr_iterator_type) :: iter
570 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new
571 TYPE(dbcsr_type), POINTER :: tmp
572 TYPE(mp_comm_type) :: group
573
574 CALL timeset(routinen, handle)
575
576 CALL get_rtp(rtp=rtp, rho_new=rho_new)
577
578 DO i = 1, SIZE(rho_new)
579 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 DO i = 1, SIZE(delta_p)
583 !square all entries of both matrices
584 CALL dbcsr_iterator_start(iter, delta_p(i)%matrix)
585 DO WHILE (dbcsr_iterator_blocks_left(iter))
586 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
587 block_values = block_values*block_values
588 END DO
589 CALL dbcsr_iterator_stop(iter)
590 END DO
591 NULLIFY (tmp)
592 ALLOCATE (tmp)
593 CALL dbcsr_create(tmp, template=delta_p(1)%matrix, matrix_type="N")
594 DO ispin = 1, SIZE(delta_p)/2
595 CALL dbcsr_desymmetrize(delta_p(2*ispin - 1)%matrix, tmp)
596 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 max_alpha = zero
600 DO ispin = 1, SIZE(delta_p)/2
601 CALL dbcsr_iterator_start(iter, delta_p(2*ispin)%matrix)
602 DO WHILE (dbcsr_iterator_blocks_left(iter))
603 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
604 alpha = maxval(block_values)
605 IF (alpha > max_alpha) max_alpha = alpha
606 END DO
607 CALL dbcsr_iterator_stop(iter)
608 END DO
609 CALL dbcsr_get_info(delta_p(1)%matrix, group=group)
610 CALL group%max(max_alpha)
611 delta_eps = sqrt(max_alpha)
613 CALL timestop(handle)
614
615 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 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 CALL timeset(routinen, handle)
634
635 NULLIFY (dft_control)
636
637 logger => cp_get_default_logger()
638 output_unit = cp_logger_get_default_io_unit(logger)
639 CALL get_qs_env(qs_env, dft_control=dft_control)
640 IF (dft_control%qs_control%dftb) THEN
641 CALL scf_post_calculation_tb(qs_env, "DFTB", .false.)
642 ELSE IF (dft_control%qs_control%xtb) THEN
643 CALL scf_post_calculation_tb(qs_env, "xTB", .false.)
644 ELSE
645 CALL qs_scf_post_moments(qs_env%input, logger, qs_env, output_unit)
646 END IF
647 CALL timestop(handle)
648
649 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmp
669
670 CALL timeset(routinen, handle)
671
672 logger => cp_get_default_logger()
673 unit_nr = cp_logger_get_default_io_unit(logger)
674 NULLIFY (tmp)
675 CALL dbcsr_allocate_matrix_set(tmp, SIZE(rho))
676 DO i = 1, SIZE(rho)
677 CALL dbcsr_init_p(tmp(i)%matrix)
678 CALL dbcsr_create(tmp(i)%matrix, template=rho(i)%matrix)
679 CALL dbcsr_copy(tmp(i)%matrix, rho(i)%matrix)
680 END DO
681 DO ispin = 1, SIZE(rho)/2
682 re = 2*ispin - 1
683 im = 2*ispin
684 eps = max(filter_eps, 1.0e-11_dp)
685 DO WHILE (eps < 1.1_dp)
686 CALL dbcsr_filter(tmp(re)%matrix, eps)
687 occ = dbcsr_get_occupation(tmp(re)%matrix)
688 IF (unit_nr > 0) WRITE (unit_nr, fmt="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
689 ispin, " eps ", eps, " real: ", occ
690 eps = eps*10
691 END DO
692 eps = max(filter_eps, 1.0e-11_dp)
693 DO WHILE (eps < 1.1_dp)
694 CALL dbcsr_filter(tmp(im)%matrix, eps)
695 occ = dbcsr_get_occupation(tmp(im)%matrix)
696 IF (unit_nr > 0) WRITE (unit_nr, fmt="((T3,A,I1,A,F15.12,A,T61,F20.10))") "Occupation of rho spin ", &
697 ispin, " eps ", eps, " imag: ", occ
698 eps = eps*10.0_dp
699 END DO
700 END DO
702 CALL timestop(handle)
703
704 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 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 CALL timeset(routinen, handle)
726 logger => cp_get_default_logger()
727 IF (logger%para_env%is_source()) THEN
728 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
729 ELSE
730 unit_nr = -1
731 END IF
732
733 project_name = logger%iter_info%project_name
734 DO ispin = 1, SIZE(rho_new)/2
735 re = 2*ispin - 1
736 im = 2*ispin
737 IF (history) THEN
738 WRITE (file_name, '(A,I0,A)') &
739 trim(project_name)//"_LS_DM_SPIN_RE", ispin, "_"//trim(cp_iter_string(logger%iter_info))//"_RESTART.dm"
740 ELSE
741 WRITE (file_name, '(A,I0,A)') trim(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
742 END IF
743 cs_pos = dbcsr_checksum(rho_new(re)%matrix, pos=.true.)
744 IF (unit_nr > 0) THEN
745 WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//trim(file_name)//" with checksum: ", cs_pos
746 END IF
747 CALL dbcsr_binary_write(rho_new(re)%matrix, file_name)
748 IF (history) THEN
749 WRITE (file_name, '(A,I0,A)') &
750 trim(project_name)//"_LS_DM_SPIN_IM", ispin, "_"//trim(cp_iter_string(logger%iter_info))//"_RESTART.dm"
751 ELSE
752 WRITE (file_name, '(A,I0,A)') trim(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
753 END IF
754 cs_pos = dbcsr_checksum(rho_new(im)%matrix, pos=.true.)
755 IF (unit_nr > 0) THEN
756 WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//trim(file_name)//" with checksum: ", cs_pos
757 END IF
758 CALL dbcsr_binary_write(rho_new(im)%matrix, file_name)
759 END DO
760
761 CALL timestop(handle)
762
763 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 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 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 CALL timeset(routinen, handle)
799
800 logger => cp_get_default_logger()
801 CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
802 CALL qs_subsys_get(subsys, particles=particles)
803 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
804
805 NULLIFY (zero, tmp)
806 ALLOCATE (zero, tmp)
807 CALL dbcsr_create(zero, template=p_im)
808 CALL dbcsr_copy(zero, p_im)
809 CALL dbcsr_set(zero, 0.0_dp)
810 CALL dbcsr_create(tmp, template=p_im)
811 CALL dbcsr_copy(tmp, p_im)
812 IF (nspin == 1) THEN
813 CALL dbcsr_scale(tmp, 0.5_dp)
814 END IF
815 current_env%gauge = -1
816 current_env%gauge_init = .false.
817 CALL auxbas_pw_pool%create_pw(rs)
818 CALL auxbas_pw_pool%create_pw(gs)
819
820 NULLIFY (stride)
821 ALLOCATE (stride(3))
822
823 DO dir = 1, 3
824
825 CALL pw_zero(rs)
826 CALL pw_zero(gs)
827
828 CALL calculate_jrho_resp(zero, tmp, zero, zero, dir, dir, rs, gs, qs_env, current_env, retain_rsgrid=.true.)
829
830 stride = section_get_ivals(dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT%STRIDE")
831
832 IF (dir == 1) THEN
833 sdir = "-x"
834 ELSEIF (dir == 2) THEN
835 sdir = "-y"
836 ELSE
837 sdir = "-z"
838 END IF
839 WRITE (char_spin, "(I1)") spin
840
841 ext = "-SPIN-"//char_spin//sdir//".cube"
842 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 log_filename=.false., mpi_io=mpi_io)
846
847 CALL cp_pw_to_cube(rs, print_unit, "EMD current", particles=particles, stride=stride, &
848 mpi_io=mpi_io)
849
850 CALL cp_print_key_finished_output(print_unit, logger, dft_section, "REAL_TIME_PROPAGATION%PRINT%CURRENT", &
851 mpi_io=mpi_io)
852
853 END DO
854
855 CALL auxbas_pw_pool%give_back_pw(rs)
856 CALL auxbas_pw_pool%give_back_pw(gs)
857
858 CALL dbcsr_deallocate_matrix(zero)
860
861 DEALLOCATE (stride)
862
863 CALL timestop(handle)
864
865 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 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 CALL timeset(routinen, handle)
887
888 CALL get_qs_env(qs_env, scf_env=scf_env)
889 IF (rtp%linear_scaling) THEN
890 CALL write_mo_free_results(qs_env)
891 ELSE
892 CALL write_mo_free_results(qs_env)
893 CALL write_mo_dependent_results(qs_env, scf_env)
894 ! Time-dependent MO print
895 CALL write_rtp_mos_to_output_unit(qs_env, rtp)
896 CALL write_rtp_mo_cubes(qs_env, rtp)
897 END IF
898
899 CALL timestop(handle)
900
901 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 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 NULLIFY (dft_control)
925
926 logger => cp_get_default_logger()
927 output_unit = cp_logger_get_default_io_unit(logger)
928
929 CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp)
930
931 i_step = rtp%istep
932
933 unit_nr = cp_print_key_unit_nr(logger, dft_section, &
934 "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat", is_new_file=new_file)
935
936 IF (output_unit > 0) THEN
937 rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
938 IF (unit_nr /= output_unit) THEN
939 INQUIRE (unit=unit_nr, name=filename)
940 WRITE (unit=output_unit, fmt="(/,T2,A,2(/,T3,A),/)") &
941 "FIELD", "The field applied is written to the file:", &
942 trim(filename)
943 ELSE
944 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 (trim(rlab(i)), "=", dft_control%rtp_control%field(i), i=1, 3)
947 END IF
948
949 IF (new_file) THEN
950 IF (dft_control%apply_efield_field) THEN
951 WRITE (unit=unit_nr, fmt='("#",5X,A,8X,A,3(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z"
952 ELSE IF (dft_control%apply_vector_potential) THEN
953 WRITE (unit=unit_nr, fmt='("#",5X,A,8X,A,6(6X,A))') "Step Nr.", "Time[fs]", " Field X", " Field Y", " Field Z", &
954 " Vec. Pot. X", " Vec. Pot. Y", " Vec. Pot. Z"
955 END IF
956 END IF
957
958 field = 0.0_dp
959 IF (dft_control%apply_efield_field) THEN
960 CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
961 WRITE (unit=unit_nr, fmt="(I10,F16.6,3(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
962 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 ELSE IF (dft_control%apply_vector_potential) THEN
967 WRITE (unit=unit_nr, fmt="(I10,F16.6,6(F16.8,1X))") qs_env%sim_step, qs_env%sim_time*femtoseconds, &
968 dft_control%rtp_control%field(1), dft_control%rtp_control%field(2), dft_control%rtp_control%field(3), &
969 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 "REAL_TIME_PROPAGATION%PRINT%FIELD")
976
977 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 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 NULLIFY (dft_control, energy, rtp)
999
1000 logger => cp_get_default_logger()
1001 output_unit = cp_logger_get_default_io_unit(logger)
1002
1003 CALL get_qs_env(qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
1004 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 file_action="WRITE", is_new_file=new_file)
1009
1010 IF (output_unit > 0) THEN
1011 IF (unit_nr /= output_unit) THEN
1012 INQUIRE (unit=unit_nr, name=filename)
1013 WRITE (unit=output_unit, fmt="(/,T2,A,2(/,T3,A),/)") &
1014 "ENERGY_CONSTITUENTS", "Total Energy constituents written to file:", &
1015 trim(filename)
1016 ELSE
1017 WRITE (unit=output_unit, fmt="(/,T2,A)") "ENERGY_CONSTITUENTS"
1018 END IF
1019
1020 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 WRITE (unit=unit_nr, fmt='("#",5X,A,8X,A,10(6X,A))') "Step Nr.", "Time[fs]", &
1024 "Total ener.[a.u.]", "core[a.u.] ", " overlap [a.u.]", "hartree[a.u.]", " exc. [a.u.] ", &
1025 " 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 qs_env%sim_step, qs_env%sim_time*femtoseconds, &
1030 energy%total, energy%core, energy%core_overlap, energy%hartree, energy%exc, &
1031 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 "REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS")
1037
1038 END SUBROUTINE print_rtp_energy_components
1039
1040END MODULE rt_propagation_output
Define the atomic kind types and their sub types.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_binary_write(matrix, filepath)
...
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
real(kind=dp) function, public dbcsr_checksum(matrix, pos)
Calculates the checksum of a DBCSR matrix.
subroutine, public dbcsr_trace(matrix, trace)
Computes the trace of the given matrix, also known as the sum of its diagonal elements.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
basic linear algebra operations for full matrices
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_double(fmstruct, struct, context, col, row)
creates a struct with twice the number of blocks on each core. If matrix A has to be multiplied with ...
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
character(len=default_string_length) function, public cp_iter_string(iter_info, print_key, for_file)
returns the iteration string, a string that is useful to create unique filenames (once you trim it)
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
all routins needed for a nonperiodic electric field
subroutine, public make_field(dft_control, field, sim_step, sim_time)
computes the amplitude of the efield within a given envelop
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public ehrenfest
integer, parameter, public real_time_propagation
objects that represent the structure of input sections and the data contained in an input section
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_path_length
Definition kinds.F:58
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:130
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public femtoseconds
Definition physcon.F:153
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
given the response wavefunctions obtained by the application of the (rxp), p, and ((dk-dl)xp) operato...
subroutine, public calculate_jrho_resp(mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, ib, idir, current_rs, current_gs, qs_env, current_env, soft_valid, retain_rsgrid)
Calculation of the idir component of the response current density in the presence of a constant magne...
Type definitiona for linear response calculations.
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
subroutine, public write_rt_mos_to_restart(mo_array, rt_mos, particle_set, dft_section, qs_kind_set)
...
Definition qs_mo_io.F:233
Define the neighbor list data types and the corresponding functionality.
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Does all kind of post scf calculations for GPW/GAPW.
subroutine, public write_mo_free_results(qs_env)
Write QS results always available (if switched on through the print_keys) Can be called from ls_scf.
subroutine, public qs_scf_post_moments(input, logger, qs_env, output_unit)
Computes and prints electric moments.
subroutine, public write_mo_dependent_results(qs_env, scf_env)
Write QS results available if MO's are present (if switched on through the print_keys) Writes only MO...
Does all kind of post scf calculations for DFTB.
subroutine, public scf_post_calculation_tb(qs_env, tb_type, no_mos)
collects possible post - scf calculations and prints info / computes properties.
module that contains the definitions of the scf types
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Function related to MO projection in RTP calculations.
subroutine, public compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
Compute the projection of the current MO coefficients on reference ones and write the results.
Routine for the real time propagation output.
subroutine, public report_density_occupation(filter_eps, rho)
Reports the sparsity pattern of the complex density matrix.
subroutine, public rt_prop_output(qs_env, run_type, delta_iter, used_time)
...
subroutine, public rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
computes the convergence criterion for RTP and EMD
subroutine, public rt_convergence_density(rtp, delta_p, delta_eps)
computes the convergence criterion for RTP and EMD based on the density matrix
Types and set_get for real time propagation depending on runtype and diagonalization method different...
subroutine, public get_rtp(rtp, exp_h_old, exp_h_new, h_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, s_inv, s_half, s_minus_half, b_mat, c_mat, propagator_matrix, mixing, mixing_factor, s_der, dt, nsteps, sinvh, sinvh_imag, sinvb, admm_mos)
...
Routines needed for EMD.
subroutine, public write_rtp_mos_to_output_unit(qs_env, rtp)
...
subroutine, public write_rtp_mo_cubes(qs_env, rtp)
Write the time dependent amplitude of the MOs in real grid. Very close to qs_scf_post_gpw/qs_scf_post...
Provides all information about an atomic kind.
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
keeps the density in various representations, keeping track of which ones are valid.