1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6MODULE qs_tddfpt2_operators 7 USE admm_types, ONLY: admm_type 8 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,& 9 cp_dbcsr_sm_fm_multiply 10 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 11 cp_fm_scale_and_add 12 USE cp_fm_struct, ONLY: cp_fm_struct_type 13 USE cp_fm_types, ONLY: cp_fm_create,& 14 cp_fm_get_info,& 15 cp_fm_p_type,& 16 cp_fm_release,& 17 cp_fm_to_fm,& 18 cp_fm_type 19 USE cp_gemm_interface, ONLY: cp_gemm 20 USE dbcsr_api, ONLY: dbcsr_p_type,& 21 dbcsr_set 22 USE hfx_admm_utils, ONLY: tddft_hfx_matrix 23 USE kinds, ONLY: dp 24 USE pw_env_types, ONLY: pw_env_get,& 25 pw_env_type 26 USE pw_methods, ONLY: pw_axpy,& 27 pw_transfer,& 28 pw_zero 29 USE pw_poisson_methods, ONLY: pw_poisson_solve 30 USE pw_poisson_types, ONLY: pw_poisson_type 31 USE pw_pool_types, ONLY: pw_pool_create_pw,& 32 pw_pool_give_back_pw,& 33 pw_pool_type 34 USE pw_types, ONLY: REALDATA3D,& 35 REALSPACE,& 36 pw_p_type,& 37 pw_release,& 38 pw_type 39 USE qs_environment_types, ONLY: get_qs_env,& 40 qs_environment_type 41 USE qs_rho_types, ONLY: qs_rho_get,& 42 qs_rho_type 43 USE qs_tddfpt2_types, ONLY: full_kernel_env_type,& 44 tddfpt_ground_state_mos 45 USE xc, ONLY: xc_calc_2nd_deriv,& 46 xc_vxc_pw_create 47 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 48 xc_rho_set_type,& 49 xc_rho_set_update 50#include "./base/base_uses.f90" 51 52 IMPLICIT NONE 53 54 PRIVATE 55 56 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_operators' 57 58 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 59 ! number of first derivative components (3: d/dx, d/dy, d/dz) 60 INTEGER, PARAMETER, PRIVATE :: nderivs = 3 61 INTEGER, PARAMETER, PRIVATE :: maxspins = 2 62 63 PUBLIC :: tddfpt_apply_energy_diff, tddfpt_apply_coulomb, tddfpt_apply_xc, tddfpt_apply_hfx 64 65! ************************************************************************************************** 66 67CONTAINS 68 69! ************************************************************************************************** 70!> \brief Apply orbital energy difference term: 71!> Aop_evects(spin,state) += KS(spin) * evects(spin,state) - 72!> S * evects(spin,state) * diag(evals_occ(spin)) 73!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit) 74!> \param evects trial vectors C_{1,i} 75!> \param S_evects S * C_{1,i} 76!> \param gs_mos molecular orbitals optimised for the ground state (only occupied orbital 77!> energies [component %evals_occ] are needed) 78!> \param matrix_ks Kohn-Sham matrix 79!> \par History 80!> * 05.2016 initialise all matrix elements in one go [Sergey Chulkov] 81!> * 03.2017 renamed from tddfpt_init_energy_diff(), altered prototype [Sergey Chulkov] 82!> \note Based on the subroutine p_op_l1() which was originally created by 83!> Thomas Chassaing on 08.2002. 84! ************************************************************************************************** 85 SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks) 86 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects, S_evects 87 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 88 INTENT(in) :: gs_mos 89 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks 90 91 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_energy_diff', & 92 routineP = moduleN//':'//routineN 93 94 INTEGER :: handle, ispin, ivect, nactive, nao, & 95 nspins, nvects 96 TYPE(cp_fm_struct_type), POINTER :: matrix_struct 97 TYPE(cp_fm_type), POINTER :: hevec 98 99 CALL timeset(routineN, handle) 100 101 nspins = SIZE(evects, 1) 102 nvects = SIZE(evects, 2) 103 104 DO ispin = 1, nspins 105 CALL cp_fm_get_info(matrix=evects(ispin, 1)%matrix, matrix_struct=matrix_struct, & 106 nrow_global=nao, ncol_global=nactive) 107 NULLIFY (hevec) 108 CALL cp_fm_create(hevec, matrix_struct) 109 110 DO ivect = 1, nvects 111 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, evects(ispin, ivect)%matrix, & 112 Aop_evects(ispin, ivect)%matrix, ncol=nactive, & 113 alpha=1.0_dp, beta=1.0_dp) 114 115 IF (ASSOCIATED(gs_mos(ispin)%evals_occ_matrix)) THEN 116 ! orbital energy correction: evals_occ_matrix is not a diagonal matrix 117 CALL cp_gemm('N', 'N', nao, nactive, nactive, 1.0_dp, & 118 S_evects(ispin, ivect)%matrix, gs_mos(ispin)%evals_occ_matrix, & 119 0.0_dp, hevec) 120 ELSE 121 CALL cp_fm_to_fm(S_evects(ispin, ivect)%matrix, hevec) 122 CALL cp_fm_column_scale(hevec, gs_mos(ispin)%evals_occ) 123 END IF 124 125 ! KS * C1 - S * C1 * occupied_orbital_energies 126 CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect)%matrix, -1.0_dp, hevec) 127 END DO 128 CALL cp_fm_release(hevec) 129 END DO 130 131 CALL timestop(handle) 132 133 END SUBROUTINE tddfpt_apply_energy_diff 134 135! ************************************************************************************************** 136!> \brief Update v_rspace by adding coulomb term. 137!> \param A_ia_rspace action of TDDFPT operator on the trial vector expressed in a plane wave 138!> representation (modified on exit) 139!> \param rho_ia_g response density in reciprocal space for the given trial vector 140!> \param pw_env plain wave environment 141!> \param work_v_gspace work reciprocal-space grid to store Coulomb potential (modified on exit) 142!> \param work_v_rspace work real-space grid to store Coulomb potential (modified on exit) 143!> \par History 144!> * 05.2016 compute all coulomb terms in one go [Sergey Chulkov] 145!> * 03.2017 proceed excited states sequentially; minimise the number of conversions between 146!> DBCSR and FM matrices [Sergey Chulkov] 147!> * 06.2018 return the action expressed in the plane wave representation instead of the one 148!> in the atomic basis set representation 149!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by 150!> Mohamed Fawzi on 10.2002. 151! ************************************************************************************************** 152 SUBROUTINE tddfpt_apply_coulomb(A_ia_rspace, rho_ia_g, pw_env, work_v_gspace, work_v_rspace) 153 TYPE(pw_p_type), DIMENSION(:), POINTER :: A_ia_rspace 154 TYPE(pw_type), POINTER :: rho_ia_g 155 TYPE(pw_env_type), POINTER :: pw_env 156 TYPE(pw_p_type), INTENT(inout) :: work_v_gspace, work_v_rspace 157 158 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_coulomb', & 159 routineP = moduleN//':'//routineN 160 161 INTEGER :: handle, ispin, nspins 162 REAL(kind=dp) :: alpha, pair_energy 163 TYPE(pw_poisson_type), POINTER :: poisson_env 164 165 CALL timeset(routineN, handle) 166 167 nspins = SIZE(A_ia_rspace) 168 CALL pw_env_get(pw_env, poisson_env=poisson_env) 169 170 IF (nspins > 1) THEN 171 alpha = 1.0_dp 172 ELSE 173 ! spin-restricted case: alpha == 2 due to singlet state. 174 ! In case of triplet states alpha == 0, so we should not call this subroutine at all. 175 alpha = 2.0_dp 176 END IF 177 178 CALL pw_poisson_solve(poisson_env, rho_ia_g, pair_energy, work_v_gspace%pw) 179 CALL pw_transfer(work_v_gspace%pw, work_v_rspace%pw) 180 181 ! (i a || j b) = ( i_alpha a_alpha + i_beta a_beta || j_alpha b_alpha + j_beta b_beta) = 182 ! tr (Cj_alpha^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_alpha) + 183 ! tr (Cj_beta^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_beta) 184 DO ispin = 1, nspins 185 CALL pw_axpy(work_v_rspace%pw, A_ia_rspace(ispin)%pw, alpha) 186 END DO 187 188 CALL timestop(handle) 189 190 END SUBROUTINE tddfpt_apply_coulomb 191 192! ************************************************************************************************** 193!> \brief Driver routine for applying fxc (analyic vs. finite difference for testing 194!> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave 195!> representation (modified on exit) 196!> \param kernel_env kernel environment 197!> \param rho_ia_struct response density for the given trial vector 198!> \param is_rks_triplets indicates that the triplet excited states calculation using 199!> spin-unpolarised molecular orbitals has been requested 200!> \param pw_env plain wave environment 201!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation 202!> potential with respect to the response density (modified on exit) 203! ************************************************************************************************** 204 SUBROUTINE tddfpt_apply_xc(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, pw_env, work_v_xc) 205 TYPE(pw_p_type), DIMENSION(:), POINTER :: A_ia_rspace 206 TYPE(full_kernel_env_type), INTENT(in) :: kernel_env 207 TYPE(qs_rho_type), POINTER :: rho_ia_struct 208 LOGICAL, INTENT(in) :: is_rks_triplets 209 TYPE(pw_env_type), POINTER :: pw_env 210 TYPE(pw_p_type), DIMENSION(:), POINTER :: work_v_xc 211 212 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc', & 213 routineP = moduleN//':'//routineN 214 215 IF (kernel_env%deriv2_analytic) THEN 216 CALL tddfpt_apply_xc_analytic(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, & 217 pw_env, work_v_xc) 218 ELSE 219 CALL tddfpt_apply_xc_fd(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, & 220 pw_env, work_v_xc) 221 END IF 222 223 END SUBROUTINE tddfpt_apply_xc 224 225! ************************************************************************************************** 226!> \brief Update A_ia_munu by adding exchange-correlation term. 227!> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave 228!> representation (modified on exit) 229!> \param kernel_env kernel environment 230!> \param rho_ia_struct response density for the given trial vector 231!> \param is_rks_triplets indicates that the triplet excited states calculation using 232!> spin-unpolarised molecular orbitals has been requested 233!> \param pw_env plain wave environment 234!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation 235!> potential with respect to the response density (modified on exit) 236!> \par History 237!> * 05.2016 compute all kernel terms in one go [Sergey Chulkov] 238!> * 03.2017 proceed excited states sequentially; minimise the number of conversions between 239!> DBCSR and FM matrices [Sergey Chulkov] 240!> * 06.2018 return the action expressed in the plane wave representation instead of the one 241!> in the atomic basis set representation 242!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by 243!> Mohamed Fawzi on 10.2002. 244! ************************************************************************************************** 245 SUBROUTINE tddfpt_apply_xc_analytic(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, pw_env, work_v_xc) 246 TYPE(pw_p_type), DIMENSION(:), POINTER :: A_ia_rspace 247 TYPE(full_kernel_env_type), INTENT(in) :: kernel_env 248 TYPE(qs_rho_type), POINTER :: rho_ia_struct 249 LOGICAL, INTENT(in) :: is_rks_triplets 250 TYPE(pw_env_type), POINTER :: pw_env 251 TYPE(pw_p_type), DIMENSION(:), POINTER :: work_v_xc 252 253 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc_analytic', & 254 routineP = moduleN//':'//routineN 255 256 INTEGER :: handle, ispin, nspins 257 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g2, rho_ia_r, & 258 rho_ia_r2, tau_ia_r, tau_ia_r2 259 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 260 261 CALL timeset(routineN, handle) 262 263 nspins = SIZE(A_ia_rspace) 264 CALL qs_rho_get(rho_ia_struct, rho_g=rho_ia_g, rho_r=rho_ia_r, tau_r=tau_ia_r) 265 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 266 267 IF (debug_this_module) THEN 268 CPASSERT(SIZE(rho_ia_g) == nspins) 269 CPASSERT(SIZE(rho_ia_r) == nspins) 270 CPASSERT((.NOT. ASSOCIATED(tau_ia_r)) .OR. SIZE(tau_ia_r) == nspins) 271 CPASSERT((.NOT. is_rks_triplets) .OR. nspins == 1) 272 END IF 273 274 NULLIFY (tau_ia_r2) 275 IF (is_rks_triplets) THEN 276 ALLOCATE (rho_ia_r2(2)) 277 ALLOCATE (rho_ia_g2(2)) 278 rho_ia_r2(1)%pw => rho_ia_r(1)%pw 279 rho_ia_r2(2)%pw => rho_ia_r(1)%pw 280 rho_ia_g2(1)%pw => rho_ia_g(1)%pw 281 rho_ia_g2(2)%pw => rho_ia_g(1)%pw 282 283 IF (ASSOCIATED(tau_ia_r)) THEN 284 ALLOCATE (tau_ia_r2(2)) 285 tau_ia_r2(1)%pw => tau_ia_r(1)%pw 286 tau_ia_r2(2)%pw => tau_ia_r(1)%pw 287 END IF 288 ELSE 289 ALLOCATE (rho_ia_r2(nspins)) 290 ALLOCATE (rho_ia_g2(nspins)) 291 DO ispin = 1, nspins 292 rho_ia_r2(ispin)%pw => rho_ia_r(ispin)%pw 293 rho_ia_g2(ispin)%pw => rho_ia_g(ispin)%pw 294 END DO 295 296 IF (ASSOCIATED(tau_ia_r)) THEN 297 ALLOCATE (tau_ia_r2(nspins)) 298 DO ispin = 1, nspins 299 tau_ia_r2(ispin)%pw => tau_ia_r(ispin)%pw 300 END DO 301 END IF 302 END IF 303 304 DO ispin = 1, nspins 305 CALL pw_zero(work_v_xc(ispin)%pw) 306 END DO 307 308 CALL xc_rho_set_update(rho_set=kernel_env%xc_rho1_set, rho_r=rho_ia_r2, rho_g=rho_ia_g2, tau=tau_ia_r2, & 309 needs=kernel_env%xc_rho1_cflags, xc_deriv_method_id=kernel_env%deriv_method_id, & 310 xc_rho_smooth_id=kernel_env%rho_smooth_id, pw_pool=auxbas_pw_pool) 311 312 CALL xc_calc_2nd_deriv(v_xc=work_v_xc, deriv_set=kernel_env%xc_deriv_set, rho_set=kernel_env%xc_rho_set, & 313 rho1_set=kernel_env%xc_rho1_set, pw_pool=auxbas_pw_pool, & 314 xc_section=kernel_env%xc_section, gapw=.FALSE., tddfpt_fac=kernel_env%beta) 315 316 DO ispin = 1, nspins 317 ! pw2 = pw2 + alpha * pw1 318 CALL pw_axpy(work_v_xc(ispin)%pw, A_ia_rspace(ispin)%pw, kernel_env%alpha) 319 END DO 320 321 DEALLOCATE (rho_ia_g2, rho_ia_r2) 322 323 CALL timestop(handle) 324 325 END SUBROUTINE tddfpt_apply_xc_analytic 326 327! ************************************************************************************************** 328!> \brief Update A_ia_munu by adding exchange-correlation term using finite difference methods. 329!> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave 330!> representation (modified on exit) 331!> \param kernel_env kernel environment 332!> \param rho_ia_struct response density for the given trial vector 333!> \param is_rks_triplets indicates that the triplet excited states calculation using 334!> spin-unpolarised molecular orbitals has been requested 335!> \param pw_env plain wave environment 336!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation 337!> potential with respect to the response density (modified on exit) 338! ************************************************************************************************** 339 SUBROUTINE tddfpt_apply_xc_fd(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, pw_env, work_v_xc) 340 TYPE(pw_p_type), DIMENSION(:), POINTER :: A_ia_rspace 341 TYPE(full_kernel_env_type), INTENT(in) :: kernel_env 342 TYPE(qs_rho_type), POINTER :: rho_ia_struct 343 LOGICAL, INTENT(in) :: is_rks_triplets 344 TYPE(pw_env_type), POINTER :: pw_env 345 TYPE(pw_p_type), DIMENSION(:), POINTER :: work_v_xc 346 347 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc_fd', & 348 routineP = moduleN//':'//routineN 349 REAL(KIND=dp), PARAMETER :: h = 0.001_dp 350 351 INTEGER :: handle, ispin, nspins 352 LOGICAL :: lsd, singlet, triplet 353 REAL(KIND=dp) :: exc 354 REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc 355 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho, rhoa, rhob 356 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1, rho_g, rho_r, tau, vxc_rho_1, & 357 vxc_rho_2, vxc_rho_3, vxc_rho_4, & 358 vxc_tau 359 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 360 TYPE(xc_rho_set_type), POINTER :: rho_set 361 362 CALL timeset(routineN, handle) 363 364 nspins = SIZE(A_ia_rspace) 365 CALL qs_rho_get(rho_ia_struct, rho_r=rho1, tau_r=tau) 366 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 367 DO ispin = 1, nspins 368 CALL pw_zero(work_v_xc(ispin)%pw) 369 END DO 370 rho_set => kernel_env%xc_rho_set 371 372 singlet = .FALSE. 373 triplet = .FALSE. 374 lsd = .FALSE. 375 IF (nspins == 1 .AND. .NOT. is_rks_triplets) THEN 376 singlet = .TRUE. 377 ELSE IF (nspins == 1 .AND. is_rks_triplets) THEN 378 triplet = .TRUE. 379 ELSE IF (nspins == 2) THEN 380 lsd = .TRUE. 381 ELSE 382 CPABORT("illegal options") 383 END IF 384 385 IF (ASSOCIATED(tau)) THEN 386 CPABORT("Tau (meta) functionals not possible") 387 END IF 388 NULLIFY (vxc_tau) 389 390 IF (singlet) THEN 391 NULLIFY (vxc_rho_1, vxc_rho_2, rho_g) 392 ALLOCATE (rho_r(1)) 393 NULLIFY (rho_r(1)%pw) 394 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(1)%pw, in_space=REALSPACE, use_data=REALDATA3D) 395 CALL xc_rho_set_get(rho_set, rho=rho) 396 rho_r(1)%pw%cr3d(:, :, :) = rho(:, :, :) + 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 397 CALL xc_vxc_pw_create(vxc_rho_1, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 398 auxbas_pw_pool, .FALSE., virial_xc) 399 rho_r(1)%pw%cr3d(:, :, :) = rho(:, :, :) - 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 400 CALL xc_vxc_pw_create(vxc_rho_2, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 401 auxbas_pw_pool, .FALSE., virial_xc) 402 work_v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :) - vxc_rho_2(1)%pw%cr3d(:, :, :))/h 403 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(1)%pw) 404 DEALLOCATE (rho_r) 405 CALL pw_axpy(work_v_xc(1)%pw, A_ia_rspace(1)%pw, kernel_env%alpha) 406 CALL pw_release(vxc_rho_1(1)%pw) 407 CALL pw_release(vxc_rho_2(1)%pw) 408 DEALLOCATE (vxc_rho_1, vxc_rho_2) 409 ELSE IF (triplet) THEN 410 NULLIFY (vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4, rho_g) 411 ALLOCATE (rho_r(2)) 412 DO ispin = 1, 2 413 NULLIFY (rho_r(ispin)%pw) 414 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(ispin)%pw, in_space=REALSPACE, use_data=REALDATA3D) 415 END DO 416 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob) 417 ! K(alpha,alpha) 418 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) + 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 419 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) 420 CALL xc_vxc_pw_create(vxc_rho_1, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 421 auxbas_pw_pool, .FALSE., virial_xc) 422 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) - 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 423 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) 424 CALL xc_vxc_pw_create(vxc_rho_2, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 425 auxbas_pw_pool, .FALSE., virial_xc) 426 work_v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :) - vxc_rho_2(1)%pw%cr3d(:, :, :))/h 427 ! K(alpha,beta) 428 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) 429 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) + 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 430 CALL xc_vxc_pw_create(vxc_rho_3, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 431 auxbas_pw_pool, .FALSE., virial_xc) 432 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) 433 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) - 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 434 CALL xc_vxc_pw_create(vxc_rho_4, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 435 auxbas_pw_pool, .FALSE., virial_xc) 436 work_v_xc(1)%pw%cr3d(:, :, :) = work_v_xc(1)%pw%cr3d(:, :, :) - & 437 (vxc_rho_3(1)%pw%cr3d(:, :, :) - vxc_rho_4(1)%pw%cr3d(:, :, :))/h 438 DO ispin = 1, 2 439 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(ispin)%pw) 440 END DO 441 DEALLOCATE (rho_r) 442 CALL pw_axpy(work_v_xc(1)%pw, A_ia_rspace(1)%pw, kernel_env%alpha) 443 DO ispin = 1, 2 444 CALL pw_release(vxc_rho_1(ispin)%pw) 445 CALL pw_release(vxc_rho_2(ispin)%pw) 446 CALL pw_release(vxc_rho_3(ispin)%pw) 447 CALL pw_release(vxc_rho_4(ispin)%pw) 448 END DO 449 DEALLOCATE (vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4) 450 ELSE IF (lsd) THEN 451 NULLIFY (vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4, rho_g) 452 ALLOCATE (rho_r(2)) 453 DO ispin = 1, 2 454 NULLIFY (rho_r(ispin)%pw) 455 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(ispin)%pw, in_space=REALSPACE, use_data=REALDATA3D) 456 END DO 457 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob) 458 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) + 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 459 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) + 0.5_dp*h*rho1(2)%pw%cr3d(:, :, :) 460 CALL xc_vxc_pw_create(vxc_rho_1, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 461 auxbas_pw_pool, .FALSE., virial_xc) 462 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) - 0.5_dp*h*rho1(1)%pw%cr3d(:, :, :) 463 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) - 0.5_dp*h*rho1(2)%pw%cr3d(:, :, :) 464 CALL xc_vxc_pw_create(vxc_rho_2, vxc_tau, exc, rho_r, rho_g, tau, kernel_env%xc_section, & 465 auxbas_pw_pool, .FALSE., virial_xc) 466 work_v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :) - vxc_rho_2(1)%pw%cr3d(:, :, :))/h 467 work_v_xc(2)%pw%cr3d(:, :, :) = (vxc_rho_1(2)%pw%cr3d(:, :, :) - vxc_rho_2(2)%pw%cr3d(:, :, :))/h 468 DO ispin = 1, 2 469 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(ispin)%pw) 470 END DO 471 DEALLOCATE (rho_r) 472 CALL pw_axpy(work_v_xc(1)%pw, A_ia_rspace(1)%pw, kernel_env%alpha) 473 CALL pw_axpy(work_v_xc(2)%pw, A_ia_rspace(2)%pw, kernel_env%alpha) 474 DO ispin = 1, 2 475 CALL pw_release(vxc_rho_1(ispin)%pw) 476 CALL pw_release(vxc_rho_2(ispin)%pw) 477 END DO 478 DEALLOCATE (vxc_rho_1, vxc_rho_2) 479 END IF 480 481 CALL timestop(handle) 482 483 END SUBROUTINE tddfpt_apply_xc_fd 484 485! ************************************************************************************************** 486!> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term. 487!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit) 488!> \param evects trial vectors 489!> \param gs_mos molecular orbitals optimised for the ground state (only occupied 490!> molecular orbitals [component %mos_occ] are needed) 491!> \param do_admm perform auxiliary density matrix method calculations 492!> \param qs_env Quickstep environment 493!> \param work_rho_ia_ao work sparse matrix with shape [nao x nao] distributed globally 494!> to store response density (modified on exit) 495!> \param work_hmat work sparse matrix with shape [nao x nao] distributed globally 496!> (modified on exit) 497!> \param wfm_rho_orb work dense matrix with shape [nao x nao] distributed globally 498!> (modified on exit) 499!> \par History 500!> * 05.2016 compute all exact-exchange terms in one go [Sergey Chulkov] 501!> * 03.2017 code related to ADMM correction is now moved to tddfpt_apply_admm_correction() 502!> in order to compute this correction within parallel groups [Sergey Chulkov] 503!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by 504!> Mohamed Fawzi on 10.2002. 505! ************************************************************************************************** 506 SUBROUTINE tddfpt_apply_hfx(Aop_evects, evects, gs_mos, do_admm, qs_env, & 507 work_rho_ia_ao, work_hmat, wfm_rho_orb) 508 TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects 509 TYPE(tddfpt_ground_state_mos), DIMENSION(:), & 510 INTENT(in) :: gs_mos 511 LOGICAL, INTENT(in) :: do_admm 512 TYPE(qs_environment_type), POINTER :: qs_env 513 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: work_rho_ia_ao, work_hmat 514 TYPE(cp_fm_type), POINTER :: wfm_rho_orb 515 516 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_hfx', & 517 routineP = moduleN//':'//routineN 518 519 INTEGER :: handle, ispin, ivect, nao, nao_aux, & 520 nspins, nvects 521 INTEGER, DIMENSION(maxspins) :: nactive 522 REAL(kind=dp) :: alpha 523 TYPE(admm_type), POINTER :: admm_env 524 525 CALL timeset(routineN, handle) 526 527 nspins = SIZE(evects, 1) 528 nvects = SIZE(evects, 2) 529 530 IF (nspins > 1) THEN 531 alpha = 2.0_dp 532 ELSE 533 alpha = 4.0_dp 534 END IF 535 536 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao) 537 DO ispin = 1, nspins 538 CALL cp_fm_get_info(evects(ispin, 1)%matrix, ncol_global=nactive(ispin)) 539 END DO 540 541 IF (do_admm) THEN 542 CALL get_qs_env(qs_env, admm_env=admm_env) 543 CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux) 544 END IF 545 546 ! some stuff from qs_ks_build_kohn_sham_matrix 547 ! TO DO: add SIC support 548 DO ivect = 1, nvects 549 DO ispin = 1, nspins 550 CALL cp_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, gs_mos(ispin)%mos_occ, & 551 evects(ispin, ivect)%matrix, 0.0_dp, wfm_rho_orb) 552 CALL cp_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect)%matrix, & 553 gs_mos(ispin)%mos_occ, 1.0_dp, wfm_rho_orb) 554 555 CALL dbcsr_set(work_hmat(ispin)%matrix, 0.0_dp) 556 IF (do_admm) THEN 557 CALL cp_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, & 558 wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb) 559 CALL cp_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%A, admm_env%work_aux_orb, & 560 0.0_dp, admm_env%work_aux_aux) 561 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.) 562 ELSE 563 CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.) 564 END IF 565 END DO 566 567 CALL tddft_hfx_matrix(work_hmat, work_rho_ia_ao, qs_env) 568 569 IF (do_admm) THEN 570 DO ispin = 1, nspins 571 CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, & 572 ncol=nao, alpha=1.0_dp, beta=0.0_dp) 573 574 CALL cp_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, & 575 admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb) 576 577 CALL cp_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, & 578 gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect)%matrix) 579 END DO 580 ELSE 581 DO ispin = 1, nspins 582 CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, gs_mos(ispin)%mos_occ, & 583 Aop_evects(ispin, ivect)%matrix, ncol=nactive(ispin), & 584 alpha=alpha, beta=1.0_dp) 585 END DO 586 END IF 587 END DO 588 589 CALL timestop(handle) 590 591 END SUBROUTINE tddfpt_apply_hfx 592 593END MODULE qs_tddfpt2_operators 594