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