1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for propagating the orbitals 8!> \author Florian Schiffmann (02.09) 9! ************************************************************************************************** 10MODULE rt_propagation_methods 11 USE bibliography, ONLY: Kolafa2004,& 12 cite_reference 13 USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,& 14 cp_cfm_gemm,& 15 cp_cfm_triangular_multiply 16 USE cp_cfm_types, ONLY: cp_cfm_create,& 17 cp_cfm_release,& 18 cp_cfm_type 19 USE cp_control_types, ONLY: rtp_control_type 20 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,& 21 cp_dbcsr_cholesky_invert 22 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,& 23 dbcsr_allocate_matrix_set,& 24 dbcsr_deallocate_matrix_set 25 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add 26 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 27 cp_fm_struct_double,& 28 cp_fm_struct_release,& 29 cp_fm_struct_type 30 USE cp_fm_types, ONLY: cp_fm_create,& 31 cp_fm_get_info,& 32 cp_fm_p_type,& 33 cp_fm_release,& 34 cp_fm_to_fm,& 35 cp_fm_type 36 USE cp_fm_vect, ONLY: cp_fm_vect_dealloc 37 USE cp_log_handling, ONLY: cp_get_default_logger,& 38 cp_logger_get_default_unit_nr,& 39 cp_logger_type,& 40 cp_to_string 41 USE dbcsr_api, ONLY: & 42 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, & 43 dbcsr_filter, dbcsr_frobenius_norm, dbcsr_get_block_p, dbcsr_init_p, & 44 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, & 45 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, & 46 dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type 47 USE input_constants, ONLY: do_arnoldi,& 48 do_bch,& 49 do_em,& 50 do_pade,& 51 do_taylor 52 USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz 53 USE kinds, ONLY: dp 54 USE ls_matrix_exp, ONLY: cp_complex_dbcsr_gemm_3 55 USE mathlib, ONLY: binomial 56 USE qs_energy_init, ONLY: qs_energies_init 57 USE qs_energy_types, ONLY: qs_energy_type 58 USE qs_environment_types, ONLY: get_qs_env,& 59 qs_environment_type 60 USE qs_ks_methods, ONLY: qs_ks_update_qs_env 61 USE rt_make_propagators, ONLY: propagate_arnoldi,& 62 propagate_bch,& 63 propagate_exp,& 64 propagate_exp_density 65 USE rt_propagation_output, ONLY: report_density_occupation,& 66 rt_convergence,& 67 rt_convergence_density 68 USE rt_propagation_types, ONLY: get_rtp,& 69 rt_prop_type 70 USE rt_propagation_utils, ONLY: calc_S_derivs,& 71 calc_update_rho,& 72 calc_update_rho_sparse 73#include "../base/base_uses.f90" 74 75 IMPLICIT NONE 76 77 PRIVATE 78 79 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_methods' 80 81 PUBLIC :: propagation_step, & 82 s_matrices_create, & 83 calc_sinvH, & 84 put_data_to_history 85 86CONTAINS 87 88! ************************************************************************************************** 89!> \brief performes a single propagation step a(t+Dt)=U(t+Dt,t)*a(0) 90!> and calculates the new exponential 91!> \param qs_env ... 92!> \param rtp ... 93!> \param rtp_control ... 94!> \author Florian Schiffmann (02.09) 95! ************************************************************************************************** 96 97 SUBROUTINE propagation_step(qs_env, rtp, rtp_control) 98 99 TYPE(qs_environment_type), POINTER :: qs_env 100 TYPE(rt_prop_type), POINTER :: rtp 101 TYPE(rtp_control_type), POINTER :: rtp_control 102 103 CHARACTER(len=*), PARAMETER :: routineN = 'propagation_step', & 104 routineP = moduleN//':'//routineN 105 106 INTEGER :: aspc_order, handle, i, im, re, unit_nr 107 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: delta_mos, mos_new 108 TYPE(cp_logger_type), POINTER :: logger 109 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P, H_last_iter, ks_mix, ks_mix_im, & 110 matrix_ks, matrix_ks_im, matrix_s, & 111 rho_new 112 113 CALL timeset(routineN, handle) 114 115 logger => cp_get_default_logger() 116 IF (logger%para_env%ionode) THEN 117 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 118 ELSE 119 unit_nr = -1 120 ENDIF 121 122 NULLIFY (delta_P, rho_new, delta_mos, mos_new) 123 NULLIFY (ks_mix, ks_mix_im) 124 ! get everything needed and set some values 125 CALL get_qs_env(qs_env, matrix_s=matrix_s) 126 IF (rtp%iter == 1) THEN 127 CALL qs_energies_init(qs_env, .FALSE.) 128 CALL get_qs_env(qs_env, matrix_s=matrix_s) 129 IF (.NOT. rtp_control%fixed_ions) THEN 130 CALL s_matrices_create(matrix_s, rtp) 131 END IF 132 rtp%delta_iter = 100.0_dp 133 rtp%mixing_factor = 1.0_dp 134 rtp%mixing = .FALSE. 135 aspc_order = rtp_control%aspc_order 136 CALL aspc_extrapolate(rtp, matrix_s, aspc_order) 137 IF (rtp%linear_scaling) THEN 138 CALL calc_update_rho_sparse(qs_env) 139 ELSE 140 CALL calc_update_rho(qs_env) 141 ENDIF 142 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.) 143 END IF 144 IF (.NOT. rtp_control%fixed_ions) THEN 145 CALL calc_S_derivs(qs_env) 146 END IF 147 rtp%converged = .FALSE. 148 149 IF (rtp%linear_scaling) THEN 150 ! keep temporary copy of the starting density matrix to check for convergence 151 CALL get_rtp(rtp=rtp, rho_new=rho_new) 152 NULLIFY (delta_P) 153 CALL dbcsr_allocate_matrix_set(delta_P, SIZE(rho_new)) 154 DO i = 1, SIZE(rho_new) 155 CALL dbcsr_init_p(delta_P(i)%matrix) 156 CALL dbcsr_create(delta_P(i)%matrix, template=rho_new(i)%matrix) 157 CALL dbcsr_copy(delta_P(i)%matrix, rho_new(i)%matrix) 158 END DO 159 ELSE 160 ! keep temporary copy of the starting mos to check for convergence 161 CALL get_rtp(rtp=rtp, mos_new=mos_new) 162 ALLOCATE (delta_mos(SIZE(mos_new))) 163 DO i = 1, SIZE(mos_new) 164 CALL cp_fm_create(delta_mos(i)%matrix, & 165 matrix_struct=mos_new(i)%matrix%matrix_struct, & 166 name="delta_mos"//TRIM(ADJUSTL(cp_to_string(i)))) 167 CALL cp_fm_to_fm(mos_new(i)%matrix, delta_mos(i)%matrix) 168 END DO 169 ENDIF 170 171 CALL get_qs_env(qs_env, & 172 matrix_ks=matrix_ks, & 173 matrix_ks_im=matrix_ks_im) 174 175 CALL get_rtp(rtp=rtp, H_last_iter=H_last_iter) 176 IF (rtp%mixing) THEN 177 IF (unit_nr > 0) THEN 178 WRITE (unit_nr, '(t3,a,2f16.8)') "Mixing the Hamiltonians to improve robustness, mixing factor: ", rtp%mixing_factor 179 ENDIF 180 CALL dbcsr_allocate_matrix_set(ks_mix, SIZE(matrix_ks)) 181 CALL dbcsr_allocate_matrix_set(ks_mix_im, SIZE(matrix_ks)) 182 DO i = 1, SIZE(matrix_ks) 183 CALL dbcsr_init_p(ks_mix(i)%matrix) 184 CALL dbcsr_create(ks_mix(i)%matrix, template=matrix_ks(1)%matrix) 185 CALL dbcsr_init_p(ks_mix_im(i)%matrix) 186 CALL dbcsr_create(ks_mix_im(i)%matrix, template=matrix_ks(1)%matrix) 187 ENDDO 188 DO i = 1, SIZE(matrix_ks) 189 re = 2*i - 1 190 im = 2*i 191 CALL dbcsr_add(ks_mix(i)%matrix, matrix_ks(i)%matrix, 0.0_dp, rtp%mixing_factor) 192 CALL dbcsr_add(ks_mix(i)%matrix, H_last_iter(re)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor) 193 IF (rtp%do_hfx) THEN 194 CALL dbcsr_add(ks_mix_im(i)%matrix, matrix_ks_im(i)%matrix, 0.0_dp, rtp%mixing_factor) 195 CALL dbcsr_add(ks_mix_im(i)%matrix, H_last_iter(im)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor) 196 ENDIF 197 ENDDO 198 CALL calc_SinvH(rtp, ks_mix, ks_mix_im, rtp_control) 199 DO i = 1, SIZE(matrix_ks) 200 re = 2*i - 1 201 im = 2*i 202 CALL dbcsr_copy(H_last_iter(re)%matrix, ks_mix(i)%matrix) 203 IF (rtp%do_hfx) THEN 204 CALL dbcsr_copy(H_last_iter(im)%matrix, ks_mix_im(i)%matrix) 205 ENDIF 206 ENDDO 207 CALL dbcsr_deallocate_matrix_set(ks_mix) 208 CALL dbcsr_deallocate_matrix_set(ks_mix_im) 209 ELSE 210 CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control) 211 DO i = 1, SIZE(matrix_ks) 212 re = 2*i - 1 213 im = 2*i 214 CALL dbcsr_copy(H_last_iter(re)%matrix, matrix_ks(i)%matrix) 215 IF (rtp%do_hfx) THEN 216 CALL dbcsr_copy(H_last_iter(im)%matrix, matrix_ks_im(i)%matrix) 217 ENDIF 218 ENDDO 219 ENDIF 220 221 CALL compute_propagator_matrix(rtp, rtp_control%propagator) 222 223 SELECT CASE (rtp_control%mat_exp) 224 CASE (do_pade, do_taylor) 225 IF (rtp%linear_scaling) THEN 226 CALL propagate_exp_density(rtp, rtp_control) 227 CALL calc_update_rho_sparse(qs_env) 228 ELSE 229 CALL propagate_exp(rtp, rtp_control) 230 CALL calc_update_rho(qs_env) 231 END IF 232 CASE (do_arnoldi) 233 CALL propagate_arnoldi(rtp, rtp_control) 234 CALL calc_update_rho(qs_env) 235 CASE (do_bch) 236 CALL propagate_bch(rtp, rtp_control) 237 CALL calc_update_rho_sparse(qs_env) 238 END SELECT 239 CALL step_finalize(qs_env, rtp_control, delta_mos, delta_P) 240 IF (rtp%linear_scaling) THEN 241 CALL dbcsr_deallocate_matrix_set(delta_P) 242 ELSE 243 CALL cp_fm_vect_dealloc(delta_mos) 244 ENDIF 245 246 CALL timestop(handle) 247 248 END SUBROUTINE propagation_step 249 250! ************************************************************************************************** 251!> \brief Performes all the stuff to finish the step: 252!> convergence checks 253!> copying stuff into right place for the next step 254!> updating the history for extrapolation 255!> \param qs_env ... 256!> \param rtp_control ... 257!> \param delta_mos ... 258!> \param delta_P ... 259!> \author Florian Schiffmann (02.09) 260! ************************************************************************************************** 261 262 SUBROUTINE step_finalize(qs_env, rtp_control, delta_mos, delta_P) 263 TYPE(qs_environment_type), POINTER :: qs_env 264 TYPE(rtp_control_type), POINTER :: rtp_control 265 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: delta_mos 266 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P 267 268 CHARACTER(len=*), PARAMETER :: routineN = 'step_finalize', routineP = moduleN//':'//routineN 269 270 INTEGER :: handle, i, ihist 271 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new, mos_old 272 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, matrix_ks, & 273 matrix_ks_im, rho_new, rho_old, s_mat 274 TYPE(qs_energy_type), POINTER :: energy 275 TYPE(rt_prop_type), POINTER :: rtp 276 277 CALL timeset(routineN, handle) 278 279 CALL get_qs_env(qs_env=qs_env, rtp=rtp, matrix_s=s_mat, matrix_ks=matrix_ks, matrix_ks_im=matrix_ks_im, energy=energy) 280 CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, exp_H_new=exp_H_new) 281 282 IF (rtp_control%sc_check_start .LT. rtp%iter) THEN 283 rtp%delta_iter_old = rtp%delta_iter 284 IF (rtp%linear_scaling) THEN 285 CALL rt_convergence_density(rtp, delta_P, rtp%delta_iter) 286 ELSE 287 CALL rt_convergence(rtp, s_mat(1)%matrix, delta_mos, rtp%delta_iter) 288 END IF 289 rtp%converged = (rtp%delta_iter .LT. rtp_control%eps_ener) 290 !Apply mixing if scf loop is not converging 291 292 !It would be better to redo the current step with mixixng, 293 !but currently the decision is made to use mixing from the next step on 294 IF (rtp_control%sc_check_start .LT. rtp%iter + 1) THEN 295 IF (rtp%delta_iter/rtp%delta_iter_old > 0.9) THEN 296 rtp%mixing_factor = MAX(rtp%mixing_factor/2.0_dp, 0.125_dp) 297 rtp%mixing = .TRUE. 298 ENDIF 299 ENDIF 300 END IF 301 302 IF (rtp%converged) THEN 303 IF (rtp%linear_scaling) THEN 304 CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new) 305 CALL purify_mcweeny_complex_nonorth(rho_new, s_mat, rtp%filter_eps, rtp%filter_eps_small, & 306 rtp_control%mcweeny_max_iter, rtp_control%mcweeny_eps) 307 IF (rtp_control%mcweeny_max_iter > 0) CALL calc_update_rho_sparse(qs_env) 308 CALL report_density_occupation(rtp%filter_eps, rho_new) 309 DO i = 1, SIZE(rho_new) 310 CALL dbcsr_copy(rho_old(i)%matrix, rho_new(i)%matrix) 311 END DO 312 ELSE 313 CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new) 314 DO i = 1, SIZE(mos_new) 315 CALL cp_fm_to_fm(mos_new(i)%matrix, mos_old(i)%matrix) 316 END DO 317 ENDIF 318 IF (rtp_control%propagator == do_em) CALL calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control) 319 DO i = 1, SIZE(exp_H_new) 320 CALL dbcsr_copy(exp_H_old(i)%matrix, exp_H_new(i)%matrix) 321 END DO 322 ihist = MOD(rtp%istep, rtp_control%aspc_order) + 1 323 IF (rtp_control%fixed_ions) THEN 324 CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, ihist=ihist) 325 ELSE 326 CALL put_data_to_history(rtp, rho=rho_new, mos=mos_new, s_mat=s_mat, ihist=ihist) 327 END IF 328 END IF 329 330 rtp%energy_new = energy%total 331 332 CALL timestop(handle) 333 334 END SUBROUTINE step_finalize 335 336! ************************************************************************************************** 337!> \brief computes the propagator matrix for EM/ETRS, RTP/EMD 338!> \param rtp ... 339!> \param propagator ... 340!> \author Florian Schiffmann (02.09) 341! ************************************************************************************************** 342 343 SUBROUTINE compute_propagator_matrix(rtp, propagator) 344 TYPE(rt_prop_type), POINTER :: rtp 345 INTEGER :: propagator 346 347 CHARACTER(len=*), PARAMETER :: routineN = 'compute_propagator_matrix', & 348 routineP = moduleN//':'//routineN 349 350 INTEGER :: handle, i 351 REAL(Kind=dp) :: dt, prefac 352 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, propagator_matrix 353 354 CALL timeset(routineN, handle) 355 CALL get_rtp(rtp=rtp, exp_H_new=exp_H_new, exp_H_old=exp_H_old, & 356 propagator_matrix=propagator_matrix, dt=dt) 357 358 prefac = -0.5_dp*dt 359 360 DO i = 1, SIZE(exp_H_new) 361 CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_new(i)%matrix, 0.0_dp, prefac) 362 IF (propagator == do_em) & 363 CALL dbcsr_add(propagator_matrix(i)%matrix, exp_H_old(i)%matrix, 1.0_dp, prefac) 364 END DO 365 366 CALL timestop(handle) 367 368 END SUBROUTINE compute_propagator_matrix 369 370! ************************************************************************************************** 371!> \brief computes t*S_inv*H, if needed t*Sinv*B 372!> \param rtp ... 373!> \param matrix_ks ... 374!> \param matrix_ks_im ... 375!> \param rtp_control ... 376!> \author Florian Schiffmann (02.09) 377! ************************************************************************************************** 378 379 SUBROUTINE calc_SinvH(rtp, matrix_ks, matrix_ks_im, rtp_control) 380 TYPE(rt_prop_type), POINTER :: rtp 381 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_im 382 TYPE(rtp_control_type), POINTER :: rtp_control 383 384 CHARACTER(len=*), PARAMETER :: routineN = 'calc_SinvH', routineP = moduleN//':'//routineN 385 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 386 387 INTEGER :: handle, im, ispin, re 388 REAL(dp) :: t 389 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H, SinvB, SinvH 390 TYPE(dbcsr_type) :: matrix_ks_nosym 391 TYPE(dbcsr_type), POINTER :: B_mat, S_inv, S_minus_half 392 393 CALL timeset(routineN, handle) 394 CALL get_rtp(rtp=rtp, S_inv=S_inv, S_minus_half=S_minus_half, exp_H_new=exp_H, dt=t) 395 CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks(1)%matrix, matrix_type="N") 396 DO ispin = 1, SIZE(matrix_ks) 397 re = ispin*2 - 1 398 im = ispin*2 399 CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, matrix_ks_nosym) 400 CALL dbcsr_multiply("N", "N", one, S_inv, matrix_ks_nosym, zero, exp_H(im)%matrix, & 401 filter_eps=rtp%filter_eps) 402 IF (.NOT. rtp_control%fixed_ions) THEN 403 CALL get_rtp(rtp=rtp, SinvH=SinvH) 404 CALL dbcsr_copy(SinvH(ispin)%matrix, exp_H(im)%matrix) 405 END IF 406 END DO 407 IF (.NOT. rtp_control%fixed_ions .OR. rtp%do_hfx) THEN 408 CALL get_rtp(rtp=rtp, B_mat=B_mat, SinvB=SinvB) 409 IF (rtp%do_hfx) THEN 410 DO ispin = 1, SIZE(matrix_ks) 411 re = ispin*2 - 1 412 im = ispin*2 413 CALL dbcsr_set(matrix_ks_nosym, 0.0_dp) 414 CALL dbcsr_desymmetrize(matrix_ks_im(ispin)%matrix, matrix_ks_nosym) 415 416 ! take care of the EMD case and add the velocity scaled S_derivative 417 IF (.NOT. rtp_control%fixed_ions) & 418 CALL dbcsr_add(matrix_ks_nosym, B_mat, 1.0_dp, -1.0_dp) 419 420 CALL dbcsr_multiply("N", "N", -one, S_inv, matrix_ks_nosym, zero, exp_H(re)%matrix, & 421 filter_eps=rtp%filter_eps) 422 423 IF (.NOT. rtp_control%fixed_ions) & 424 CALL dbcsr_copy(SinvB(ispin)%matrix, exp_H(re)%matrix) 425 END DO 426 ELSE 427 ! in case of pure EMD its only needed once as B is the same for both spins 428 CALL dbcsr_multiply("N", "N", one, S_inv, B_mat, zero, exp_H(1)%matrix, filter_eps=rtp%filter_eps) 429 430 CALL dbcsr_copy(SinvB(1)%matrix, exp_H(1)%matrix) 431 432 IF (SIZE(matrix_ks) == 2) CALL dbcsr_copy(exp_H(3)%matrix, exp_H(1)%matrix) 433 IF (SIZE(matrix_ks) == 2) CALL dbcsr_copy(SinvB(2)%matrix, SinvB(1)%matrix) 434 END IF 435 ELSE 436 !set real part to zero 437 DO ispin = 1, SIZE(exp_H)/2 438 re = ispin*2 - 1 439 im = ispin*2 440 CALL dbcsr_set(exp_H(re)%matrix, zero) 441 ENDDO 442 END IF 443 CALL dbcsr_release(matrix_ks_nosym) 444 CALL timestop(handle) 445 END SUBROUTINE calc_SinvH 446 447! ************************************************************************************************** 448!> \brief calculates the needed overlap-like matrices 449!> depending on the way the exponential is calculated, only S^-1 is needed 450!> \param s_mat ... 451!> \param rtp ... 452!> \author Florian Schiffmann (02.09) 453! ************************************************************************************************** 454 455 SUBROUTINE s_matrices_create(s_mat, rtp) 456 457 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_mat 458 TYPE(rt_prop_type), POINTER :: rtp 459 460 CHARACTER(len=*), PARAMETER :: routineN = 's_matrices_create', & 461 routineP = moduleN//':'//routineN 462 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 463 464 INTEGER :: handle 465 TYPE(dbcsr_type), POINTER :: S_half, S_inv, S_minus_half 466 467 CALL timeset(routineN, handle) 468 469 CALL get_rtp(rtp=rtp, S_inv=S_inv) 470 471 IF (rtp%linear_scaling) THEN 472 CALL get_rtp(rtp=rtp, S_half=S_half, S_minus_half=S_minus_half) 473 CALL matrix_sqrt_Newton_Schulz(S_half, S_minus_half, s_mat(1)%matrix, rtp%filter_eps, & 474 rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter) 475 CALL dbcsr_multiply("N", "N", one, S_minus_half, S_minus_half, zero, S_inv, & 476 filter_eps=rtp%filter_eps) 477 ELSE 478 CALL dbcsr_copy(S_inv, s_mat(1)%matrix) 479 CALL cp_dbcsr_cholesky_decompose(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, & 480 blacs_env=rtp%ao_ao_fmstruct%context) 481 CALL cp_dbcsr_cholesky_invert(S_inv, para_env=rtp%ao_ao_fmstruct%para_env, & 482 blacs_env=rtp%ao_ao_fmstruct%context, upper_to_full=.TRUE.) 483 ENDIF 484 485 CALL timestop(handle) 486 END SUBROUTINE s_matrices_create 487 488! ************************************************************************************************** 489!> \brief Calculates the frobenius norm of a complex matrix represented by two real matrices 490!> \param frob_norm ... 491!> \param mat_re ... 492!> \param mat_im ... 493!> \author Samuel Andermatt (04.14) 494! ************************************************************************************************** 495 496 SUBROUTINE complex_frobenius_norm(frob_norm, mat_re, mat_im) 497 498 REAL(KIND=dp), INTENT(out) :: frob_norm 499 TYPE(dbcsr_type), POINTER :: mat_re, mat_im 500 501 CHARACTER(len=*), PARAMETER :: routineN = 'complex_frobenius_norm', & 502 routineP = moduleN//':'//routineN 503 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 504 505 INTEGER :: col_atom, handle, row_atom 506 LOGICAL :: found 507 REAL(dp), DIMENSION(:), POINTER :: block_values, block_values2 508 TYPE(dbcsr_iterator_type) :: iter 509 TYPE(dbcsr_type), POINTER :: tmp 510 511 CALL timeset(routineN, handle) 512 513 NULLIFY (tmp) 514 ALLOCATE (tmp) 515 CALL dbcsr_create(tmp, template=mat_re) 516 !make sure the tmp has the same sparsity pattern as the real and the complex part combined 517 CALL dbcsr_add(tmp, mat_re, zero, one) 518 CALL dbcsr_add(tmp, mat_im, zero, one) 519 CALL dbcsr_set(tmp, zero) 520 !calculate the hadamard product 521 CALL dbcsr_iterator_start(iter, tmp) 522 DO WHILE (dbcsr_iterator_blocks_left(iter)) 523 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values) 524 CALL dbcsr_get_block_p(mat_re, row_atom, col_atom, block_values2, found=found) 525 IF (found) THEN 526 block_values = block_values2*block_values2 527 ENDIF 528 CALL dbcsr_get_block_p(mat_im, row_atom, col_atom, block_values2, found=found) 529 IF (found) THEN 530 block_values = block_values + block_values2*block_values2 531 ENDIF 532 block_values = SQRT(block_values) 533 END DO 534 CALL dbcsr_iterator_stop(iter) 535 frob_norm = dbcsr_frobenius_norm(tmp) 536 537 CALL dbcsr_deallocate_matrix(tmp) 538 539 CALL timestop(handle) 540 541 END SUBROUTINE complex_frobenius_norm 542 543! ************************************************************************************************** 544!> \brief Does McWeeny for complex matrices in the non-orthogonal basis 545!> \param P ... 546!> \param s_mat ... 547!> \param eps ... 548!> \param eps_small ... 549!> \param max_iter ... 550!> \param threshold ... 551!> \author Samuel Andermatt (04.14) 552! ************************************************************************************************** 553 554 SUBROUTINE purify_mcweeny_complex_nonorth(P, s_mat, eps, eps_small, max_iter, threshold) 555 556 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: P, s_mat 557 REAL(KIND=dp), INTENT(in) :: eps, eps_small 558 INTEGER, INTENT(in) :: max_iter 559 REAL(KIND=dp), INTENT(in) :: threshold 560 561 CHARACTER(len=*), PARAMETER :: routineN = 'purify_mcweeny_complex_nonorth', & 562 routineP = moduleN//':'//routineN 563 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 564 565 INTEGER :: handle, i, im, imax, ispin, re, unit_nr 566 REAL(KIND=dp) :: frob_norm 567 TYPE(cp_logger_type), POINTER :: logger 568 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: PS, PSP, tmp 569 570 CALL timeset(routineN, handle) 571 572 logger => cp_get_default_logger() 573 IF (logger%para_env%ionode) THEN 574 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 575 ELSE 576 unit_nr = -1 577 ENDIF 578 579 NULLIFY (tmp, PS, PSP) 580 CALL dbcsr_allocate_matrix_set(tmp, SIZE(P)) 581 CALL dbcsr_allocate_matrix_set(PSP, SIZE(P)) 582 CALL dbcsr_allocate_matrix_set(PS, SIZE(P)) 583 DO i = 1, SIZE(P) 584 CALL dbcsr_init_p(PS(i)%matrix) 585 CALL dbcsr_create(PS(i)%matrix, template=P(1)%matrix) 586 CALL dbcsr_init_p(PSP(i)%matrix) 587 CALL dbcsr_create(PSP(i)%matrix, template=P(1)%matrix) 588 CALL dbcsr_init_p(tmp(i)%matrix) 589 CALL dbcsr_create(tmp(i)%matrix, template=P(1)%matrix) 590 ENDDO 591 IF (SIZE(P) == 2) THEN 592 CALL dbcsr_scale(P(1)%matrix, one/2) 593 CALL dbcsr_scale(P(2)%matrix, one/2) 594 ENDIF 595 DO ispin = 1, SIZE(P)/2 596 re = 2*ispin - 1 597 im = 2*ispin 598 imax = MAX(max_iter, 1) !if max_iter is 0 then only the deviation from idempotency needs to be calculated 599 DO i = 1, imax 600 CALL dbcsr_multiply("N", "N", one, P(re)%matrix, s_mat(1)%matrix, & 601 zero, PS(re)%matrix, filter_eps=eps_small) 602 CALL dbcsr_multiply("N", "N", one, P(im)%matrix, s_mat(1)%matrix, & 603 zero, PS(im)%matrix, filter_eps=eps_small) 604 CALL cp_complex_dbcsr_gemm_3("N", "N", one, PS(re)%matrix, PS(im)%matrix, & 605 P(re)%matrix, P(im)%matrix, zero, PSP(re)%matrix, PSP(im)%matrix, & 606 filter_eps=eps_small) 607 CALL dbcsr_copy(tmp(re)%matrix, PSP(re)%matrix) 608 CALL dbcsr_copy(tmp(im)%matrix, PSP(im)%matrix) 609 CALL dbcsr_add(tmp(re)%matrix, P(re)%matrix, 1.0_dp, -1.0_dp) 610 CALL dbcsr_add(tmp(im)%matrix, P(im)%matrix, 1.0_dp, -1.0_dp) 611 CALL complex_frobenius_norm(frob_norm, tmp(re)%matrix, tmp(im)%matrix) 612 IF (unit_nr .GT. 0) WRITE (unit_nr, '(t3,a,2f16.8)') "Deviation from idempotency: ", frob_norm 613 IF (frob_norm .GT. threshold .AND. max_iter > 0) THEN 614 CALL dbcsr_copy(P(re)%matrix, PSP(re)%matrix) 615 CALL dbcsr_copy(P(im)%matrix, PSP(im)%matrix) 616 CALL cp_complex_dbcsr_gemm_3("N", "N", -2.0_dp, PS(re)%matrix, PS(im)%matrix, & 617 PSP(re)%matrix, PSP(im)%matrix, 3.0_dp, P(re)%matrix, P(im)%matrix, & 618 filter_eps=eps_small) 619 CALL dbcsr_filter(P(re)%matrix, eps) 620 CALL dbcsr_filter(P(im)%matrix, eps) 621 !make sure P is exactly hermitian 622 CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix) 623 CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2) 624 CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix) 625 CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2) 626 ELSE 627 EXIT 628 END IF 629 END DO 630 !make sure P is hermitian 631 CALL dbcsr_transposed(tmp(re)%matrix, P(re)%matrix) 632 CALL dbcsr_add(P(re)%matrix, tmp(re)%matrix, one/2, one/2) 633 CALL dbcsr_transposed(tmp(im)%matrix, P(im)%matrix) 634 CALL dbcsr_add(P(im)%matrix, tmp(im)%matrix, one/2, -one/2) 635 END DO 636 IF (SIZE(P) == 2) THEN 637 CALL dbcsr_scale(P(1)%matrix, one*2) 638 CALL dbcsr_scale(P(2)%matrix, one*2) 639 ENDIF 640 CALL dbcsr_deallocate_matrix_set(tmp) 641 CALL dbcsr_deallocate_matrix_set(PS) 642 CALL dbcsr_deallocate_matrix_set(PSP) 643 644 CALL timestop(handle) 645 646 END SUBROUTINE purify_mcweeny_complex_nonorth 647 648! ************************************************************************************************** 649!> \brief ... 650!> \param rtp ... 651!> \param matrix_s ... 652!> \param aspc_order ... 653! ************************************************************************************************** 654 SUBROUTINE aspc_extrapolate(rtp, matrix_s, aspc_order) 655 TYPE(rt_prop_type), POINTER :: rtp 656 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 657 INTEGER, INTENT(in) :: aspc_order 658 659 CHARACTER(len=*), PARAMETER :: routineN = 'aspc_extrapolate', & 660 routineP = moduleN//':'//routineN 661 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), & 662 czero = (0.0_dp, 0.0_dp) 663 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 664 665 INTEGER :: handle, i, iaspc, icol_local, ihist, & 666 imat, k, kdbl, n, naspc, ncol_local, & 667 nmat 668 REAL(KIND=dp) :: alpha 669 TYPE(cp_cfm_type), POINTER :: cfm_tmp, cfm_tmp1, csc 670 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new 671 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: mo_hist 672 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_new 673 TYPE(cp_fm_type), POINTER :: fm_tmp, fm_tmp1, fm_tmp2 674 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, s_hist 675 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_hist 676 677 NULLIFY (rho_hist) 678 CALL timeset(routineN, handle) 679 CALL cite_reference(Kolafa2004) 680 681 IF (rtp%linear_scaling) THEN 682 CALL get_rtp(rtp=rtp, rho_new=rho_new) 683 ELSE 684 CALL get_rtp(rtp=rtp, mos_new=mos_new) 685 ENDIF 686 687 naspc = MIN(rtp%istep, aspc_order) 688 IF (rtp%linear_scaling) THEN 689 nmat = SIZE(rho_new) 690 rho_hist => rtp%history%rho_history 691 DO imat = 1, nmat 692 DO iaspc = 1, naspc 693 alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* & 694 binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1) 695 ihist = MOD(rtp%istep - iaspc, aspc_order) + 1 696 IF (iaspc == 1) THEN 697 CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, zero, alpha) 698 ELSE 699 CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, one, alpha) 700 END IF 701 END DO 702 END DO 703 ELSE 704 mo_hist => rtp%history%mo_history 705 nmat = SIZE(mos_new) 706 DO imat = 1, nmat 707 DO iaspc = 1, naspc 708 alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* & 709 binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1) 710 ihist = MOD(rtp%istep - iaspc, aspc_order) + 1 711 IF (iaspc == 1) THEN 712 CALL cp_fm_scale_and_add(zero, mos_new(imat)%matrix, alpha, mo_hist(imat, ihist)%matrix) 713 ELSE 714 CALL cp_fm_scale_and_add(one, mos_new(imat)%matrix, alpha, mo_hist(imat, ihist)%matrix) 715 END IF 716 END DO 717 END DO 718 719 mo_hist => rtp%history%mo_history 720 s_hist => rtp%history%s_history 721 DO i = 1, SIZE(mos_new)/2 722 NULLIFY (matrix_struct, matrix_struct_new, csc, fm_tmp, fm_tmp1, fm_tmp2, cfm_tmp, cfm_tmp1) 723 724 CALL cp_fm_struct_double(matrix_struct, & 725 mos_new(2*i)%matrix%matrix_struct, & 726 mos_new(2*i)%matrix%matrix_struct%context, & 727 .TRUE., .FALSE.) 728 729 CALL cp_fm_create(fm_tmp, matrix_struct) 730 CALL cp_fm_create(fm_tmp1, matrix_struct) 731 CALL cp_fm_create(fm_tmp2, mos_new(2*i)%matrix%matrix_struct) 732 CALL cp_cfm_create(cfm_tmp, mos_new(2*i)%matrix%matrix_struct) 733 CALL cp_cfm_create(cfm_tmp1, mos_new(2*i)%matrix%matrix_struct) 734 735 CALL cp_fm_get_info(fm_tmp, & 736 ncol_global=kdbl) 737 738 CALL cp_fm_get_info(mos_new(2*i)%matrix, & 739 nrow_global=n, & 740 ncol_global=k, & 741 ncol_local=ncol_local) 742 743 CALL cp_fm_struct_create(matrix_struct_new, & 744 template_fmstruct=mos_new(2*i)%matrix%matrix_struct, & 745 nrow_global=k, & 746 ncol_global=k) 747 CALL cp_cfm_create(csc, matrix_struct_new) 748 749 CALL cp_fm_struct_release(matrix_struct_new) 750 CALL cp_fm_struct_release(matrix_struct) 751 752 ! first the most recent 753 754! reorthogonalize vectors 755 DO icol_local = 1, ncol_local 756 fm_tmp%local_data(:, icol_local) = mos_new(2*i - 1)%matrix%local_data(:, icol_local) 757 fm_tmp%local_data(:, icol_local + ncol_local) = mos_new(2*i)%matrix%local_data(:, icol_local) 758 END DO 759 760 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, fm_tmp, fm_tmp1, kdbl) 761 762 DO icol_local = 1, ncol_local 763 cfm_tmp%local_data(:, icol_local) = CMPLX(fm_tmp1%local_data(:, icol_local), & 764 fm_tmp1%local_data(:, icol_local + ncol_local), dp) 765 cfm_tmp1%local_data(:, icol_local) = CMPLX(mos_new(2*i - 1)%matrix%local_data(:, icol_local), & 766 mos_new(2*i)%matrix%local_data(:, icol_local), dp) 767 END DO 768 CALL cp_cfm_gemm('C', 'N', k, k, n, cone, cfm_tmp1, cfm_tmp, czero, csc) 769 CALL cp_cfm_cholesky_decompose(csc) 770 CALL cp_cfm_triangular_multiply(csc, cfm_tmp1, n_cols=k, side='R', invert_tr=.TRUE.) 771 DO icol_local = 1, ncol_local 772 mos_new(2*i - 1)%matrix%local_data(:, icol_local) = REAL(cfm_tmp1%local_data(:, icol_local), dp) 773 mos_new(2*i)%matrix%local_data(:, icol_local) = AIMAG(cfm_tmp1%local_data(:, icol_local)) 774 END DO 775 776! deallocate work matrices 777 CALL cp_cfm_release(csc) 778 CALL cp_fm_release(fm_tmp) 779 CALL cp_fm_release(fm_tmp) 780 CALL cp_fm_release(fm_tmp1) 781 CALL cp_fm_release(fm_tmp2) 782 CALL cp_cfm_release(cfm_tmp) 783 CALL cp_cfm_release(cfm_tmp1) 784 END DO 785 786 END IF 787 788 CALL timestop(handle) 789 790 END SUBROUTINE aspc_extrapolate 791 792! ************************************************************************************************** 793!> \brief ... 794!> \param rtp ... 795!> \param mos ... 796!> \param rho ... 797!> \param s_mat ... 798!> \param ihist ... 799! ************************************************************************************************** 800 SUBROUTINE put_data_to_history(rtp, mos, rho, s_mat, ihist) 801 TYPE(rt_prop_type), POINTER :: rtp 802 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos 803 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho 804 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 805 POINTER :: s_mat 806 INTEGER :: ihist 807 808 CHARACTER(len=*), PARAMETER :: routineN = 'put_data_to_history', & 809 routineP = moduleN//':'//routineN 810 811 INTEGER :: i 812 813 IF (rtp%linear_scaling) THEN 814 DO i = 1, SIZE(rho) 815 CALL dbcsr_copy(rtp%history%rho_history(i, ihist)%matrix, rho(i)%matrix) 816 END DO 817 ELSE 818 DO i = 1, SIZE(mos) 819 CALL cp_fm_to_fm(mos(i)%matrix, rtp%history%mo_history(i, ihist)%matrix) 820 END DO 821 IF (PRESENT(s_mat)) THEN 822 IF (ASSOCIATED(rtp%history%s_history(ihist)%matrix)) THEN ! the sparsity might be different 823 ! (future struct:check) 824 CALL dbcsr_deallocate_matrix(rtp%history%s_history(ihist)%matrix) 825 END IF 826 ALLOCATE (rtp%history%s_history(ihist)%matrix) 827 CALL dbcsr_copy(rtp%history%s_history(ihist)%matrix, s_mat(1)%matrix) 828 END IF 829 END IF 830 831 END SUBROUTINE put_data_to_history 832 833END MODULE rt_propagation_methods 834