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 calculating a complex matrix exponential. 8!> \author Florian Schiffmann (02.09) 9! ************************************************************************************************** 10 11MODULE rt_make_propagators 12 13 USE cp_control_types, ONLY: rtp_control_type 14 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 15 copy_fm_to_dbcsr,& 16 cp_dbcsr_sm_fm_multiply 17 USE cp_fm_types, ONLY: cp_fm_create,& 18 cp_fm_get_info,& 19 cp_fm_p_type,& 20 cp_fm_to_fm 21 USE cp_fm_vect, ONLY: cp_fm_vect_dealloc 22 USE dbcsr_api, ONLY: dbcsr_copy,& 23 dbcsr_create,& 24 dbcsr_deallocate_matrix,& 25 dbcsr_p_type,& 26 dbcsr_scale,& 27 dbcsr_type 28 USE input_constants, ONLY: do_etrs,& 29 do_pade,& 30 do_taylor 31 USE kinds, ONLY: dp 32 USE ls_matrix_exp, ONLY: bch_expansion_complex_propagator,& 33 bch_expansion_imaginary_propagator,& 34 cp_complex_dbcsr_gemm_3,& 35 taylor_full_complex_dbcsr,& 36 taylor_only_imaginary_dbcsr 37 USE matrix_exp, ONLY: arnoldi,& 38 exp_pade_full_complex,& 39 exp_pade_only_imaginary,& 40 taylor_full_complex,& 41 taylor_only_imaginary 42 USE rt_propagation_types, ONLY: get_rtp,& 43 rt_prop_type 44#include "../base/base_uses.f90" 45 46 IMPLICIT NONE 47 48 PRIVATE 49 50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_make_propagators' 51 52 PUBLIC :: propagate_exp, & 53 propagate_arnoldi, & 54 compute_exponential, & 55 compute_exponential_sparse, & 56 propagate_exp_density, & 57 propagate_bch 58 59CONTAINS 60! ************************************************************************************************** 61!> \brief performs propagations if explicit matrix exponentials are used 62!> ETRS: exp(i*H(t+dt)*dt/2)*exp(i*H(t)*dt/2)*MOS 63!> EM: exp[-idt/2H(t+dt/2)*MOS 64!> \param rtp ... 65!> \param rtp_control ... 66!> \author Florian Schiffmann (02.09) 67! ************************************************************************************************** 68 69 SUBROUTINE propagate_exp(rtp, rtp_control) 70 71 TYPE(rt_prop_type), POINTER :: rtp 72 TYPE(rtp_control_type), POINTER :: rtp_control 73 74 CHARACTER(len=*), PARAMETER :: routineN = 'propagate_exp', routineP = moduleN//':'//routineN 75 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 76 77 INTEGER :: handle, i, im, nmo, re 78 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new, mos_next, mos_old 79 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, propagator_matrix 80 81 CALL timeset(routineN, handle) 82 83 CALL get_rtp(rtp=rtp, propagator_matrix=propagator_matrix, mos_old=mos_old, mos_new=mos_new, & 84 mos_next=mos_next, exp_H_new=exp_H_new, exp_H_old=exp_H_old) 85 86 ! Only compute exponential if a new propagator matrix is available 87 CALL compute_exponential(exp_H_new, propagator_matrix, rtp_control, rtp) 88 89 DO i = 1, SIZE(mos_new)/2 90 re = 2*i - 1 91 im = 2*i 92 93 CALL cp_fm_get_info(mos_new(re)%matrix, ncol_global=nmo) 94 !Save some work by computing the first half of the propagation only once in case of ETRS 95 !For EM this matrix has to be the initial matrix, thus a copy is enough 96 IF (rtp%iter == 1) THEN 97 IF (rtp_control%propagator == do_etrs) THEN 98 CALL cp_dbcsr_sm_fm_multiply(exp_H_old(re)%matrix, mos_old(re)%matrix, & 99 mos_next(re)%matrix, nmo, alpha=one, beta=zero) 100 CALL cp_dbcsr_sm_fm_multiply(exp_H_old(im)%matrix, mos_old(im)%matrix, & 101 mos_next(re)%matrix, nmo, alpha=-one, beta=one) 102 CALL cp_dbcsr_sm_fm_multiply(exp_H_old(re)%matrix, mos_old(im)%matrix, & 103 mos_next(im)%matrix, nmo, alpha=one, beta=zero) 104 CALL cp_dbcsr_sm_fm_multiply(exp_H_old(im)%matrix, mos_old(re)%matrix, & 105 mos_next(im)%matrix, nmo, alpha=one, beta=one) 106 ELSE 107 CALL cp_fm_to_fm(mos_old(re)%matrix, mos_next(re)%matrix) 108 CALL cp_fm_to_fm(mos_old(im)%matrix, mos_next(im)%matrix) 109 END IF 110 END IF 111 CALL cp_dbcsr_sm_fm_multiply(exp_H_new(re)%matrix, mos_next(re)%matrix, & 112 mos_new(re)%matrix, nmo, alpha=one, beta=zero) 113 CALL cp_dbcsr_sm_fm_multiply(exp_H_new(im)%matrix, mos_next(im)%matrix, & 114 mos_new(re)%matrix, nmo, alpha=-one, beta=one) 115 CALL cp_dbcsr_sm_fm_multiply(exp_H_new(re)%matrix, mos_next(im)%matrix, & 116 mos_new(im)%matrix, nmo, alpha=one, beta=zero) 117 CALL cp_dbcsr_sm_fm_multiply(exp_H_new(im)%matrix, mos_next(re)%matrix, & 118 mos_new(im)%matrix, nmo, alpha=one, beta=one) 119 END DO 120 121 CALL timestop(handle) 122 123 END SUBROUTINE propagate_exp 124 125! ************************************************************************************************** 126!> \brief Propagation of the density matrix instead of the atomic orbitals 127!> via a matrix exponential 128!> \param rtp ... 129!> \param rtp_control ... 130!> \author Samuel Andermatt (02.2014) 131! ************************************************************************************************** 132 133 SUBROUTINE propagate_exp_density(rtp, rtp_control) 134 135 TYPE(rt_prop_type), POINTER :: rtp 136 TYPE(rtp_control_type), POINTER :: rtp_control 137 138 CHARACTER(len=*), PARAMETER :: routineN = 'propagate_exp_density', & 139 routineP = moduleN//':'//routineN 140 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 141 142 INTEGER :: handle, i, im, re 143 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_new, exp_H_old, propagator_matrix, & 144 rho_new, rho_next, rho_old 145 TYPE(dbcsr_type), POINTER :: tmp_im, tmp_re 146 147 CALL timeset(routineN, handle) 148 149 CALL get_rtp(rtp=rtp, propagator_matrix=propagator_matrix, exp_H_new=exp_H_new, & 150 exp_H_old=exp_H_old, rho_old=rho_old, rho_new=rho_new, rho_next=rho_next) 151 152 CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp) 153 154 !I could store these matrices in the type 155 NULLIFY (tmp_re) 156 ALLOCATE (tmp_re) 157 CALL dbcsr_create(tmp_re, template=propagator_matrix(1)%matrix, matrix_type="N") 158 NULLIFY (tmp_im) 159 ALLOCATE (tmp_im) 160 CALL dbcsr_create(tmp_im, template=propagator_matrix(1)%matrix, matrix_type="N") 161 162 DO i = 1, SIZE(exp_H_new)/2 163 re = 2*i - 1 164 im = 2*i 165 !Save some work by computing the first half of the propagation only once in case of ETRS 166 !For EM this matrix has to be the initial matrix, thus a copy is enough 167 IF (rtp%iter == 1) THEN 168 IF (rtp_control%propagator == do_etrs) THEN 169 CALL cp_complex_dbcsr_gemm_3("N", "N", one, exp_H_old(re)%matrix, exp_H_old(im)%matrix, & 170 rho_old(re)%matrix, rho_old(im)%matrix, zero, tmp_re, tmp_im, filter_eps=rtp%filter_eps) 171 CALL cp_complex_dbcsr_gemm_3("N", "C", one, tmp_re, tmp_im, exp_H_old(re)%matrix, exp_H_old(im)%matrix, & 172 zero, rho_next(re)%matrix, rho_next(im)%matrix, filter_eps=rtp%filter_eps) 173 ELSE 174 CALL dbcsr_copy(rho_next(re)%matrix, rho_old(re)%matrix) 175 CALL dbcsr_copy(rho_next(im)%matrix, rho_old(im)%matrix) 176 ENDIF 177 END IF 178 CALL cp_complex_dbcsr_gemm_3("N", "N", one, exp_H_new(re)%matrix, exp_H_new(im)%matrix, & 179 rho_next(re)%matrix, rho_next(im)%matrix, zero, tmp_re, tmp_im, filter_eps=rtp%filter_eps) 180 CALL cp_complex_dbcsr_gemm_3("N", "C", one, tmp_re, tmp_im, exp_H_new(re)%matrix, exp_H_new(im)%matrix, & 181 zero, rho_new(re)%matrix, rho_new(im)%matrix, filter_eps=rtp%filter_eps) 182 END DO 183 184 CALL dbcsr_deallocate_matrix(tmp_re) 185 CALL dbcsr_deallocate_matrix(tmp_im) 186 187 CALL timestop(handle) 188 189 END SUBROUTINE propagate_exp_density 190 191! ************************************************************************************************** 192!> \brief computes U_prop*MOs using arnoldi subspace algorithm 193!> \param rtp ... 194!> \param rtp_control ... 195!> \author Florian Schiffmann (02.09) 196! ************************************************************************************************** 197 198 SUBROUTINE propagate_arnoldi(rtp, rtp_control) 199 TYPE(rt_prop_type), POINTER :: rtp 200 TYPE(rtp_control_type), POINTER :: rtp_control 201 202 CHARACTER(len=*), PARAMETER :: routineN = 'propagate_arnoldi', & 203 routineP = moduleN//':'//routineN 204 205 INTEGER :: handle, i, im, ispin, nspin, re 206 REAL(dp) :: eps_arnoldi, t 207 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new, mos_next, mos_old, & 208 propagator_matrix_fm 209 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: propagator_matrix 210 211 CALL timeset(routineN, handle) 212 213 CALL get_rtp(rtp=rtp, dt=t, mos_new=mos_new, mos_old=mos_old, & 214 mos_next=mos_next, propagator_matrix=propagator_matrix) 215 216 nspin = SIZE(mos_new)/2 217 eps_arnoldi = rtp_control%eps_exp 218 ! except for the first step the further propagated mos_next 219 ! must be copied on mos_old so that we have the half propagated mos 220 ! ready on mos_old and only need to perform the second half propagatioon 221 IF (rtp_control%propagator == do_etrs .AND. rtp%iter == 1) THEN 222 DO i = 1, SIZE(mos_new) 223 CALL cp_fm_to_fm(mos_next(i)%matrix, mos_old(i)%matrix) 224 END DO 225 END IF 226 227 NULLIFY (propagator_matrix_fm) 228 ALLOCATE (propagator_matrix_fm(SIZE(propagator_matrix))) 229 DO i = 1, SIZE(propagator_matrix) 230 CALL cp_fm_create(propagator_matrix_fm(i)%matrix, & 231 matrix_struct=rtp%ao_ao_fmstruct, & 232 name="prop_fm") 233 CALL copy_dbcsr_to_fm(propagator_matrix(i)%matrix, propagator_matrix_fm(i)%matrix) 234 END DO 235 236 DO ispin = 1, nspin 237 re = ispin*2 - 1 238 im = ispin*2 239 IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN 240 CALL arnoldi(mos_old(re:im), mos_new(re:im), & 241 eps_arnoldi, Him=propagator_matrix_fm(im)%matrix, & 242 mos_next=mos_next(re:im), narn_old=rtp%narn_old) 243 ELSE 244 CALL arnoldi(mos_old(re:im), mos_new(re:im), & 245 eps_arnoldi, Hre=propagator_matrix_fm(re)%matrix, & 246 Him=propagator_matrix_fm(im)%matrix, mos_next=mos_next(re:im), & 247 narn_old=rtp%narn_old) 248 END IF 249 END DO 250 251! DO i=1,SIZE(propagator_matrix) 252! CALL copy_fm_to_dbcsr(propagator_matrix_fm(i)%matrix,propagator_matrix(i)%matrix) 253! END DO 254 CALL cp_fm_vect_dealloc(propagator_matrix_fm) 255 256 CALL timestop(handle) 257 258 END SUBROUTINE propagate_arnoldi 259 260! ************************************************************************************************** 261!> \brief Propagation using the Baker-Campbell-Hausdorff expansion, 262!> currently only works for rtp 263!> \param rtp ... 264!> \param rtp_control ... 265!> \author Samuel Andermatt (02.2014) 266! ************************************************************************************************** 267 268 SUBROUTINE propagate_bch(rtp, rtp_control) 269 270 TYPE(rt_prop_type), POINTER :: rtp 271 TYPE(rtp_control_type), POINTER :: rtp_control 272 273 CHARACTER(len=*), PARAMETER :: routineN = 'propagate_bch', routineP = moduleN//':'//routineN 274 275 INTEGER :: handle, im, ispin, re 276 REAL(dp) :: dt 277 REAL(KIND=dp) :: prefac 278 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: exp_H_old, propagator_matrix, rho_new, & 279 rho_next, rho_old 280 281 CALL timeset(routineN, handle) 282 283 CALL get_rtp(rtp=rtp, propagator_matrix=propagator_matrix, rho_old=rho_old, rho_new=rho_new, & 284 rho_next=rho_next) 285 286 DO ispin = 1, SIZE(propagator_matrix)/2 287 re = 2*ispin - 1 288 im = 2*ispin 289 290 IF (rtp%iter == 1) THEN 291 ! For EM I have to copy rho_old onto rho_next and for ETRS, 292 ! this is the first term of the series of commutators that result in rho_next 293 CALL dbcsr_copy(rho_next(re)%matrix, rho_old(re)%matrix) 294 CALL dbcsr_copy(rho_next(im)%matrix, rho_old(im)%matrix) 295 IF (rtp_control%propagator == do_etrs) THEN 296 !since we never calculated the matrix exponential the old matrix exponential stores the unscalled propagator 297 CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, dt=dt) 298 prefac = -0.5_dp*dt 299 CALL dbcsr_scale(exp_H_old(im)%matrix, prefac) 300 IF (.NOT. rtp%do_hfx .AND. rtp_control%fixed_ions) THEN 301 CALL bch_expansion_imaginary_propagator( & 302 exp_H_old(im)%matrix, rho_next(re)%matrix, rho_next(im)%matrix, & 303 rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp) 304 ELSE 305 CALL dbcsr_scale(exp_H_old(re)%matrix, prefac) 306 CALL bch_expansion_complex_propagator( & 307 exp_H_old(re)%matrix, exp_H_old(im)%matrix, rho_next(re)%matrix, rho_next(im)%matrix, & 308 rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp) 309 ENDIF 310 END IF 311 END IF 312 CALL dbcsr_copy(rho_new(re)%matrix, rho_next(re)%matrix) 313 CALL dbcsr_copy(rho_new(im)%matrix, rho_next(im)%matrix) 314 IF (.NOT. rtp%do_hfx .AND. rtp_control%fixed_ions) THEN 315 CALL bch_expansion_imaginary_propagator( & 316 propagator_matrix(im)%matrix, rho_new(re)%matrix, rho_new(im)%matrix, & 317 rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp) 318 ELSE 319 CALL bch_expansion_complex_propagator( & 320 propagator_matrix(re)%matrix, propagator_matrix(im)%matrix, rho_new(re)%matrix, rho_new(im)%matrix, & 321 rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp) 322 ENDIF 323 324 END DO 325 326 CALL timestop(handle) 327 328 END SUBROUTINE propagate_bch 329 330! ************************************************************************************************** 331!> \brief decides which type of exponential has to be computed 332!> \param propagator ... 333!> \param propagator_matrix ... 334!> \param rtp_control ... 335!> \param rtp ... 336!> \author Florian Schiffmann (02.09) 337! ************************************************************************************************** 338 339 SUBROUTINE compute_exponential(propagator, propagator_matrix, rtp_control, rtp) 340 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: propagator, propagator_matrix 341 TYPE(rtp_control_type), POINTER :: rtp_control 342 TYPE(rt_prop_type), POINTER :: rtp 343 344 CHARACTER(len=*), PARAMETER :: routineN = 'compute_exponential', & 345 routineP = moduleN//':'//routineN 346 347 INTEGER :: i, im, ispin, re 348 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: propagator_fm, propagator_matrix_fm 349 350 NULLIFY (propagator_fm) 351 ALLOCATE (propagator_fm(SIZE(propagator))) 352 NULLIFY (propagator_matrix_fm) 353 ALLOCATE (propagator_matrix_fm(SIZE(propagator_matrix))) 354 DO i = 1, SIZE(propagator) 355 CALL cp_fm_create(propagator_fm(i)%matrix, & 356 matrix_struct=rtp%ao_ao_fmstruct, & 357 name="prop_fm") 358 CALL copy_dbcsr_to_fm(propagator(i)%matrix, propagator_fm(i)%matrix) 359 CALL cp_fm_create(propagator_matrix_fm(i)%matrix, & 360 matrix_struct=rtp%ao_ao_fmstruct, & 361 name="prop_mat_fm") 362 CALL copy_dbcsr_to_fm(propagator_matrix(i)%matrix, propagator_matrix_fm(i)%matrix) 363 END DO 364 365 DO ispin = 1, SIZE(propagator)/2 366 re = 2*ispin - 1 367 im = 2*ispin 368 369 SELECT CASE (rtp_control%mat_exp) 370 371 CASE (do_taylor) 372 IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN 373 CALL taylor_only_imaginary(propagator_fm(re:im), propagator_matrix_fm(im)%matrix, & 374 rtp%orders(1, ispin), rtp%orders(2, ispin)) 375 ELSE 376 CALL taylor_full_complex(propagator_fm(re:im), propagator_matrix_fm(re)%matrix, propagator_matrix_fm(im)%matrix, & 377 rtp%orders(1, ispin), rtp%orders(2, ispin)) 378 END IF 379 CASE (do_pade) 380 IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN 381 CALL exp_pade_only_imaginary(propagator_fm(re:im), propagator_matrix_fm(im)%matrix, & 382 rtp%orders(1, ispin), rtp%orders(2, ispin)) 383 ELSE 384 CALL exp_pade_full_complex(propagator_fm(re:im), propagator_matrix_fm(re)%matrix, propagator_matrix_fm(im)%matrix, & 385 rtp%orders(1, ispin), rtp%orders(2, ispin)) 386 END IF 387 END SELECT 388 END DO 389 390 DO i = 1, SIZE(propagator) 391 CALL copy_fm_to_dbcsr(propagator_fm(i)%matrix, propagator(i)%matrix) 392 CALL copy_fm_to_dbcsr(propagator_matrix_fm(i)%matrix, propagator_matrix(i)%matrix) 393 END DO 394 CALL cp_fm_vect_dealloc(propagator_fm) 395 CALL cp_fm_vect_dealloc(propagator_matrix_fm) 396 397 END SUBROUTINE compute_exponential 398 399! ************************************************************************************************** 400!> \brief Sparse versions of the matrix exponentials 401!> \param propagator ... 402!> \param propagator_matrix ... 403!> \param rtp_control ... 404!> \param rtp ... 405!> \author Samuel Andermatt (02.14) 406! ************************************************************************************************** 407 408 SUBROUTINE compute_exponential_sparse(propagator, propagator_matrix, rtp_control, rtp) 409 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: propagator, propagator_matrix 410 TYPE(rtp_control_type), POINTER :: rtp_control 411 TYPE(rt_prop_type), POINTER :: rtp 412 413 CHARACTER(len=*), PARAMETER :: routineN = 'compute_exponential_sparse', & 414 routineP = moduleN//':'//routineN 415 416 INTEGER :: handle, im, ispin, re 417 418 CALL timeset(routineN, handle) 419 420 DO ispin = 1, SIZE(propagator)/2 421 re = 2*ispin - 1 422 im = 2*ispin 423 IF (rtp_control%fixed_ions .AND. .NOT. rtp%do_hfx) THEN 424 CALL taylor_only_imaginary_dbcsr(propagator(re:im), propagator_matrix(im)%matrix, & 425 rtp%orders(1, ispin), rtp%orders(2, ispin), rtp%filter_eps) 426 ELSE 427 CALL taylor_full_complex_dbcsr(propagator(re:im), propagator_matrix(re)%matrix, propagator_matrix(im)%matrix, & 428 rtp%orders(1, ispin), rtp%orders(2, ispin), rtp%filter_eps) 429 END IF 430 END DO 431 432 CALL timestop(handle) 433 434 END SUBROUTINE compute_exponential_sparse 435 436END MODULE rt_make_propagators 437