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 to compute singles correction to RPA (RSE) 8!> \par History 9!> 08.2019 created [Vladimir Rybkin] 10!> \author Vladimir Rybkin 11! ************************************************************************************************** 12MODULE rpa_rse 13 14 USE cp_blacs_env, ONLY: cp_blacs_env_type 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 17 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 18 copy_fm_to_dbcsr,& 19 dbcsr_allocate_matrix_set 20 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add 21 USE cp_fm_diag, ONLY: choose_eigv_solver 22 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 23 cp_fm_struct_release,& 24 cp_fm_struct_type 25 USE cp_fm_types, ONLY: cp_fm_create,& 26 cp_fm_get_info,& 27 cp_fm_release,& 28 cp_fm_set_all,& 29 cp_fm_to_fm_submat,& 30 cp_fm_type 31 USE cp_gemm_interface, ONLY: cp_gemm 32 USE cp_para_types, ONLY: cp_para_env_type 33 USE dbcsr_api, ONLY: dbcsr_copy,& 34 dbcsr_create,& 35 dbcsr_init_p,& 36 dbcsr_p_type,& 37 dbcsr_release,& 38 dbcsr_set,& 39 dbcsr_type_symmetric 40 USE hfx_energy_potential, ONLY: integrate_four_center 41 USE input_section_types, ONLY: section_vals_get,& 42 section_vals_get_subs_vals,& 43 section_vals_type,& 44 section_vals_val_get 45 USE kinds, ONLY: dp 46 USE message_passing, ONLY: mp_sum 47 USE mp2_types, ONLY: mp2_type 48 USE pw_types, ONLY: pw_p_type,& 49 pw_release 50 USE qs_energy_types, ONLY: qs_energy_type 51 USE qs_environment_types, ONLY: get_qs_env,& 52 qs_environment_type 53 USE qs_ks_types, ONLY: qs_ks_env_type 54 USE qs_ks_utils, ONLY: compute_matrix_vxc 55 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 56 USE qs_rho_types, ONLY: qs_rho_get,& 57 qs_rho_type 58 USE qs_vxc, ONLY: qs_vxc_create 59 60!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads 61 62#include "./base/base_uses.f90" 63 64 IMPLICIT NONE 65 66 PRIVATE 67 68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_rse' 69 70 PUBLIC :: rse_energy 71 72CONTAINS 73 74! ************************************************************************************************** 75!> \brief Single excitations energy corrections for RPA 76!> \param qs_env ... 77!> \param mp2_env ... 78!> \param para_env ... 79!> \param dft_control ... 80!> \param mo_coeff ... 81!> \param nmo ... 82!> \param homo ... 83!> \param Eigenval ... 84!> \param Eigenval_beta ... 85!> \param homo_beta ... 86!> \param mo_coeff_beta ... 87!> \author Vladimir Rybkin, 08/2019 88! ************************************************************************************************** 89 SUBROUTINE rse_energy(qs_env, mp2_env, para_env, dft_control, & 90 mo_coeff, nmo, homo, Eigenval, & 91 Eigenval_beta, homo_beta, mo_coeff_beta) 92 TYPE(qs_environment_type), POINTER :: qs_env 93 TYPE(mp2_type), POINTER :: mp2_env 94 TYPE(cp_para_env_type), POINTER :: para_env 95 TYPE(dft_control_type), POINTER :: dft_control 96 TYPE(cp_fm_type), POINTER :: mo_coeff 97 INTEGER :: nmo, homo 98 REAL(KIND=dp), DIMENSION(:) :: Eigenval 99 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: Eigenval_beta 100 INTEGER, OPTIONAL :: homo_beta 101 TYPE(cp_fm_type), OPTIONAL, POINTER :: mo_coeff_beta 102 103 CHARACTER(LEN=*), PARAMETER :: routineN = 'rse_energy', routineP = moduleN//':'//routineN 104 105 INTEGER :: dimen, handle, i_global, iiB, ispin, & 106 j_global, jjB, n_rep_hf, ncol_local, & 107 nrow_local 108 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 109 LOGICAL :: beta, do_hfx, hfx_treat_lsd_in_core 110 REAL(KIND=dp) :: coeff, rse_corr, rse_corr_beta 111 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_diff 112 TYPE(cp_blacs_env_type), POINTER :: blacs_env 113 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 114 TYPE(cp_fm_type), POINTER :: fm_ao, fm_ao_mo, fm_P_mu_nu, & 115 fm_P_mu_nu_beta, fm_X_mo, & 116 fm_X_mo_beta, fm_XC_mo, fm_XC_mo_beta 117 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_mu_nu, matrix_s, rho_ao 118 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 119 POINTER :: sab_orb 120 TYPE(qs_energy_type), POINTER :: energy 121 TYPE(qs_rho_type), POINTER :: rho 122 TYPE(section_vals_type), POINTER :: hfx_sections, input 123 124 CALL timeset(routineN, handle) 125 126 beta = .FALSE. 127 IF (PRESENT(homo_beta) .AND. PRESENT(Eigenval_beta) & 128 .AND. PRESENT(mo_coeff_beta)) beta = .TRUE. 129 130 ! Pick the diagonal terms 131 CALL cp_fm_get_info(matrix=mo_coeff, & 132 nrow_local=nrow_local, & 133 ncol_local=ncol_local, & 134 row_indices=row_indices, & 135 col_indices=col_indices) 136 137 ! start collecting stuff 138 dimen = nmo 139 NULLIFY (input, matrix_s, blacs_env, rho, energy, sab_orb) 140 CALL get_qs_env(qs_env, & 141 input=input, & 142 matrix_s=matrix_s, & 143 blacs_env=blacs_env, & 144 rho=rho, & 145 energy=energy, & 146 sab_orb=sab_orb) 147 148 CALL qs_rho_get(rho, rho_ao=rho_ao) 149 150 ! hfx section 151 NULLIFY (hfx_sections) 152 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF") 153 CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf) 154 IF (do_hfx) THEN 155 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, & 156 i_rep_section=1) 157 END IF 158 159 ! create work array 160 NULLIFY (mat_mu_nu) 161 CALL dbcsr_allocate_matrix_set(mat_mu_nu, dft_control%nspins) 162 DO ispin = 1, dft_control%nspins 163 ALLOCATE (mat_mu_nu(ispin)%matrix) 164 CALL dbcsr_create(matrix=mat_mu_nu(ispin)%matrix, template=matrix_s(1)%matrix, name="T_mu_nu", & 165 matrix_type=dbcsr_type_symmetric, nze=0) 166 CALL cp_dbcsr_alloc_block_from_nbl(mat_mu_nu(ispin)%matrix, sab_orb) 167 CALL dbcsr_set(mat_mu_nu(ispin)%matrix, 0.0_dp) 168 END DO 169 170 ! Dense (full) matrices 171 NULLIFY (fm_P_mu_nu, fm_struct_tmp) 172 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 173 nrow_global=dimen, ncol_global=dimen) 174 CALL cp_fm_create(fm_P_mu_nu, fm_struct_tmp, name="P_mu_nu") 175 CALL cp_fm_set_all(fm_P_mu_nu, 0.0_dp) 176 IF (beta) THEN 177 CALL cp_fm_create(fm_P_mu_nu_beta, fm_struct_tmp, name="P_mu_nu_beta") 178 CALL cp_fm_set_all(fm_P_mu_nu_beta, 0.0_dp) 179 ENDIF 180 CALL cp_fm_struct_release(fm_struct_tmp) 181 182 NULLIFY (fm_X_mo, fm_XC_mo, fm_struct_tmp) 183 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 184 nrow_global=dimen, ncol_global=dimen) 185 CALL cp_fm_create(fm_X_mo, fm_struct_tmp, name="f_X_mo") 186 CALL cp_fm_create(fm_XC_mo, fm_struct_tmp, name="f_XC_mo") 187 CALL cp_fm_set_all(fm_X_mo, 0.0_dp) 188 CALL cp_fm_set_all(fm_XC_mo, 0.0_dp) 189 IF (beta) THEN 190 CALL cp_fm_create(fm_X_mo_beta, fm_struct_tmp, name="f_X_mo") 191 CALL cp_fm_create(fm_XC_mo_beta, fm_struct_tmp, name="f_XC_mo") 192 CALL cp_fm_set_all(fm_X_mo_beta, 0.0_dp) 193 CALL cp_fm_set_all(fm_XC_mo_beta, 0.0_dp) 194 ENDIF 195 CALL cp_fm_struct_release(fm_struct_tmp) 196 197 NULLIFY (fm_ao) 198 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 199 nrow_global=dimen, ncol_global=dimen) 200 CALL cp_fm_create(fm_ao, fm_struct_tmp, name="f_ao") 201 CALL cp_fm_struct_release(fm_struct_tmp) 202 CALL cp_fm_set_all(fm_ao, 0.0_dp) 203 NULLIFY (fm_ao_mo) 204 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 205 nrow_global=dimen, ncol_global=dimen) 206 CALL cp_fm_create(fm_ao_mo, fm_struct_tmp, name="f_ao_mo") 207 CALL cp_fm_struct_release(fm_struct_tmp) 208 CALL cp_fm_set_all(fm_ao_mo, 0.0_dp) 209 210 ! 211 ! Ready with preparations, do the real staff 212 ! 213 214 ! Obtain density matrix like quantity 215 216 coeff = 1.0_dp 217 IF (beta) coeff = 0.5_dp 218 CALL cp_gemm(transa='N', transb='T', m=dimen, n=dimen, k=dimen, alpha=coeff, & 219 matrix_a=mo_coeff, matrix_b=mo_coeff, beta=0.0_dp, matrix_c=fm_P_mu_nu) 220 IF (beta) CALL cp_gemm(transa='N', transb='T', m=dimen, n=dimen, k=dimen, alpha=coeff, & 221 matrix_a=mo_coeff_beta, matrix_b=mo_coeff_beta, beta=0.0_dp, matrix_c=fm_P_mu_nu_beta) 222 223 ! Calculate exact exchange contribution 224 IF (.NOT. beta) THEN 225 CALL exchange_contribution(qs_env, para_env, dimen, mo_coeff, & 226 hfx_sections, n_rep_hf, & 227 rho, mat_mu_nu, fm_P_mu_nu, & 228 fm_ao, fm_X_mo, fm_ao_mo, & 229 mp2_env%ri_mp2%free_hfx_buffer) 230 ELSE 231 CALL exchange_contribution(qs_env, para_env, dimen, mo_coeff, & 232 hfx_sections, n_rep_hf, & 233 rho, mat_mu_nu, fm_P_mu_nu_beta, & 234 fm_ao, fm_X_mo, fm_ao_mo, & 235 mp2_env%ri_mp2%free_hfx_buffer, & 236 fm_X_mo_beta, fm_P_mu_nu_beta, mo_coeff_beta) 237 ENDIF 238 239 ! Calculate DFT exchange-correlation contribution 240 IF (.NOT. beta) THEN 241 CALL xc_contribution(qs_env, fm_ao, fm_ao_mo, fm_XC_mo, mo_coeff, dimen) 242 ELSE 243 CALL xc_contribution(qs_env, fm_ao, fm_ao_mo, fm_XC_mo, mo_coeff, dimen, fm_XC_mo_beta, mo_coeff_beta) 244 ENDIF 245 246 ! Compute the correction matrix: it is stored in fm_X_mo 247 CALL cp_fm_scale_and_add(1.0_dp, fm_X_mo, -1.0_dp, fm_XC_mo) 248 249 ! Pick the diagonal terms 250 CALL cp_fm_get_info(matrix=fm_X_mo, & 251 nrow_local=nrow_local, & 252 ncol_local=ncol_local, & 253 row_indices=row_indices, & 254 col_indices=col_indices) 255 256 ALLOCATE (diag_diff(dimen)) 257 diag_diff = 0.0_dp 258 259!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 260!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff,fm_X_mo) 261 DO jjB = 1, ncol_local 262 j_global = col_indices(jjB) 263 DO iiB = 1, nrow_local 264 i_global = row_indices(iiB) 265 IF (i_global .EQ. j_global) diag_diff(i_global) = fm_X_mo%local_data(iib, jjb) 266 END DO 267 END DO 268!$OMP END PARALLEL DO 269 CALL mp_sum(diag_diff, para_env%group) 270 271 ! Compute the correction 272 rse_corr = 0.0_dp 273!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 274!$OMP REDUCTION(+: rse_corr) & 275!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff, eigenval, fm_X_mo,homo) 276 DO jjB = 1, ncol_local 277 j_global = col_indices(jjB) 278 DO iiB = 1, nrow_local 279 i_global = row_indices(iiB) 280 IF ((i_global .LE. homo) .AND. (j_global .GT. homo)) THEN 281 rse_corr = rse_corr + fm_X_mo%local_data(iib, jjb)**2.0_dp/ & 282 (eigenval(i_global) - eigenval(j_global) - diag_diff(i_global) + diag_diff(j_global)) 283 ENDIF 284 END DO 285 END DO 286!$OMP END PARALLEL DO 287 CALL mp_sum(rse_corr, para_env%group) 288 289! Beta spin 290 IF (beta) THEN 291 ! Compute the correction matrix: it is stored in fm_X_mo 292 CALL cp_fm_scale_and_add(1.0_dp, fm_X_mo_beta, -1.0_dp, fm_XC_mo_beta) 293 294 rse_corr_beta = 0.0_dp 295 ! Pick the diagonal terms 296 CALL cp_fm_get_info(matrix=fm_X_mo_beta, & 297 nrow_local=nrow_local, & 298 ncol_local=ncol_local, & 299 row_indices=row_indices, & 300 col_indices=col_indices) 301 302 diag_diff = 0.0_dp 303!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 304!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff,fm_X_mo_beta) 305 DO jjB = 1, ncol_local 306 j_global = col_indices(jjB) 307 DO iiB = 1, nrow_local 308 i_global = row_indices(iiB) 309 IF (i_global .EQ. j_global) diag_diff(i_global) = fm_X_mo_beta%local_data(iib, jjb) 310 END DO 311 END DO 312!$OMP END PARALLEL DO 313 CALL mp_sum(diag_diff, para_env%group) 314 315!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 316!$OMP REDUCTION(+: rse_corr_beta) & 317!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff,eigenval_beta,fm_X_mo_Beta,homo_beta) 318 DO jjB = 1, ncol_local 319 j_global = col_indices(jjB) 320 DO iiB = 1, nrow_local 321 i_global = row_indices(iiB) 322 IF ((i_global .LE. homo_beta) .AND. (j_global .GT. homo_beta)) THEN 323 rse_corr_beta = rse_corr_beta + fm_X_mo_beta%local_data(iib, jjb)**2.0_dp/ & 324 (eigenval_beta(i_global) - eigenval_beta(j_global) - diag_diff(i_global) + diag_diff(j_global)) 325 ENDIF 326 END DO 327 END DO 328!$OMP END PARALLEL DO 329 CALL mp_sum(rse_corr_beta, para_env%group) 330 rse_corr = 0.5_dp*(rse_corr + rse_corr_beta) 331 ENDIF 332 333 mp2_env%ri_rpa%rse_corr_diag = rse_corr 334 335 ! Nondiagonal correction 336 IF (.NOT. beta) THEN 337 CALL non_diag_rse(fm_X_mo, eigenval, dimen, homo, para_env, blacs_env, rse_corr) 338 ELSE 339 CALL non_diag_rse(fm_X_mo, eigenval, dimen, homo, para_env, blacs_env, rse_corr, & 340 fm_X_mo_beta, eigenval_beta, homo_beta) 341 ENDIF 342 343 mp2_env%ri_rpa%rse_corr = rse_corr 344 345 ! Release staff 346 DEALLOCATE (diag_diff) 347 CALL cp_fm_release(fm_X_mo) 348 CALL cp_fm_release(fm_XC_mo) 349 CALL cp_fm_release(fm_ao) 350 CALL cp_fm_release(fm_P_mu_nu) 351 CALL cp_fm_release(fm_ao_mo) 352 DO ispin = 1, dft_control%nspins 353 CALL dbcsr_release(mat_mu_nu(ispin)%matrix) 354 DEALLOCATE (mat_mu_nu(ispin)%matrix) 355 ENDDO 356 DEALLOCATE (mat_mu_nu) 357 IF (beta) THEN 358 CALL cp_fm_release(fm_X_mo_beta) 359 CALL cp_fm_release(fm_XC_mo_beta) 360 CALL cp_fm_release(fm_P_mu_nu_beta) 361 ENDIF 362 363 CALL timestop(handle) 364 365 END SUBROUTINE rse_energy 366 367! ************************************************************************************************** 368!> \brief HF exchange occupied-virtual matrix 369!> \param qs_env ... 370!> \param para_env ... 371!> \param dimen ... 372!> \param mo_coeff ... 373!> \param hfx_sections ... 374!> \param n_rep_hf ... 375!> \param rho_work ... 376!> \param mat_mu_nu ... 377!> \param fm_P_mu_nu ... 378!> \param fm_X_ao ... 379!> \param fm_X_mo ... 380!> \param fm_X_ao_mo ... 381!> \param recalc_hfx_integrals ... 382!> \param fm_X_mo_beta ... 383!> \param fm_P_mu_nu_beta ... 384!> \param mo_coeff_beta ... 385! ************************************************************************************************** 386 SUBROUTINE exchange_contribution(qs_env, para_env, dimen, mo_coeff, & 387 hfx_sections, n_rep_hf, & 388 rho_work, mat_mu_nu, fm_P_mu_nu, & 389 fm_X_ao, fm_X_mo, fm_X_ao_mo, & 390 recalc_hfx_integrals, fm_X_mo_beta, & 391 fm_P_mu_nu_beta, mo_coeff_beta) 392 TYPE(qs_environment_type), POINTER :: qs_env 393 TYPE(cp_para_env_type), POINTER :: para_env 394 INTEGER :: dimen 395 TYPE(cp_fm_type), POINTER :: mo_coeff 396 TYPE(section_vals_type), POINTER :: hfx_sections 397 INTEGER :: n_rep_hf 398 TYPE(qs_rho_type), POINTER :: rho_work 399 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_mu_nu 400 TYPE(cp_fm_type), POINTER :: fm_P_mu_nu, fm_X_ao, fm_X_mo, fm_X_ao_mo 401 LOGICAL, OPTIONAL :: recalc_hfx_integrals 402 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_X_mo_beta, fm_P_mu_nu_beta, & 403 mo_coeff_beta 404 405 CHARACTER(LEN=*), PARAMETER :: routineN = 'exchange_contribution', & 406 routineP = moduleN//':'//routineN 407 408 INTEGER :: handle, irep, is, ns 409 LOGICAL :: alpha_beta, my_recalc_hfx_integrals 410 REAL(KIND=dp) :: ehfx 411 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: P_mu_nu, rho_work_ao 412 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_2d, rho_ao_2d 413 414 CALL timeset(routineN, handle) 415 416 alpha_beta = .FALSE. 417 IF (PRESENT(mo_coeff_beta)) alpha_beta = .TRUE. 418 419 my_recalc_hfx_integrals = .FALSE. 420 IF (PRESENT(recalc_hfx_integrals)) my_recalc_hfx_integrals = recalc_hfx_integrals 421 422 CALL qs_rho_get(rho_work, rho_ao=rho_work_ao) 423 ns = SIZE(rho_work_ao) 424 NULLIFY (P_mu_nu) 425 CALL dbcsr_allocate_matrix_set(P_mu_nu, ns) 426 DO is = 1, ns 427 CALL dbcsr_init_p(P_mu_nu(is)%matrix) 428 CALL dbcsr_create(P_mu_nu(is)%matrix, template=rho_work_ao(1)%matrix) 429 CALL dbcsr_copy(P_mu_nu(is)%matrix, rho_work_ao(1)%matrix) 430 CALL dbcsr_set(P_mu_nu(is)%matrix, 0.0_dp) 431 ENDDO 432 433 ! copy fm into DBCSR 434 CALL copy_fm_to_dbcsr(fm_P_mu_nu, P_mu_nu(1)%matrix, keep_sparsity=.TRUE.) 435 436 DO irep = 1, n_rep_hf 437 rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns) 438 ns = SIZE(mat_mu_nu) 439 mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns) 440 CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, & 441 para_env, my_recalc_hfx_integrals, irep, .TRUE., & 442 ispin=1) 443 END DO 444 445 ! copy back to fm 446 CALL cp_fm_set_all(fm_X_ao, 0.0_dp) 447 CALL copy_dbcsr_to_fm(matrix=mat_mu_nu(1)%matrix, fm=fm_X_ao) 448 CALL cp_fm_set_all(fm_X_mo, 0.0_dp) 449 450 ! First index 451 CALL cp_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, & 452 mo_coeff, fm_X_ao, 0.0_dp, fm_X_ao_mo) 453 454 ! Second index 455 CALL cp_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, & 456 fm_X_ao_mo, mo_coeff, 1.0_dp, fm_X_mo) 457 458 ! Beta spin 459 IF (alpha_beta) THEN 460 CALL copy_fm_to_dbcsr(fm_P_mu_nu_beta, P_mu_nu(1)%matrix, keep_sparsity=.TRUE.) 461 462 CALL dbcsr_set(mat_mu_nu(1)%matrix, 0.0_dp) 463 464 DO irep = 1, n_rep_hf 465 rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns) 466 mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns) 467 my_recalc_hfx_integrals = .FALSE. 468 CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, & 469 para_env, my_recalc_hfx_integrals, irep, .TRUE., & 470 ispin=1) 471 END DO 472 473 ! copy back to fm 474 CALL cp_fm_set_all(fm_X_ao, 0.0_dp) 475 CALL copy_dbcsr_to_fm(matrix=mat_mu_nu(1)%matrix, fm=fm_X_ao) 476 CALL cp_fm_set_all(fm_X_mo_beta, 0.0_dp) 477 478 ! First index 479 CALL cp_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, & 480 mo_coeff_beta, fm_X_ao, 0.0_dp, fm_X_ao_mo) 481 482 ! Second index 483 CALL cp_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, & 484 fm_X_ao_mo, mo_coeff_beta, 1.0_dp, fm_X_mo_beta) 485 486 ENDIF 487 488 ! Release dbcsr objects 489 DO is = 1, SIZE(P_mu_nu) 490 CALL dbcsr_release(P_mu_nu(is)%matrix) 491 DEALLOCATE (P_mu_nu(is)%matrix) 492 ENDDO 493 DEALLOCATE (P_mu_nu) 494 495 CALL timestop(handle) 496 497 END SUBROUTINE exchange_contribution 498 499! ************************************************************************************************** 500!> \brief Exchange-correlation occupied-virtual matrix 501!> \param qs_env ... 502!> \param fm_XC_ao ... 503!> \param fm_XC_ao_mo ... 504!> \param fm_XC_mo ... 505!> \param mo_coeff ... 506!> \param dimen ... 507!> \param fm_XC_mo_beta ... 508!> \param mo_coeff_beta ... 509! ************************************************************************************************** 510 SUBROUTINE xc_contribution(qs_env, fm_XC_ao, fm_XC_ao_mo, fm_XC_mo, mo_coeff, dimen, fm_XC_mo_beta, mo_coeff_beta) 511 TYPE(qs_environment_type), POINTER :: qs_env 512 TYPE(cp_fm_type), POINTER :: fm_XC_ao, fm_XC_ao_mo, fm_XC_mo, mo_coeff 513 INTEGER :: dimen 514 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_XC_mo_beta, mo_coeff_beta 515 516 CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_contribution', & 517 routineP = moduleN//':'//routineN 518 519 INTEGER :: handle, i 520 LOGICAL :: alpha_beta 521 REAL(KIND=dp) :: exc 522 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc 523 TYPE(pw_p_type), DIMENSION(:), POINTER :: tau_rspace, v_rspace 524 TYPE(qs_ks_env_type), POINTER :: ks_env 525 TYPE(qs_rho_type), POINTER :: rho 526 TYPE(section_vals_type), POINTER :: input, xc_section 527 528 CALL timeset(routineN, handle) 529 530 alpha_beta = .FALSE. 531 IF (PRESENT(fm_XC_mo_beta) .AND. PRESENT(mo_coeff_beta)) alpha_beta = .TRUE. 532 533 NULLIFY (matrix_vxc, v_rspace, tau_rspace, input, xc_section, ks_env, & 534 rho) 535 CALL get_qs_env(qs_env, matrix_vxc=matrix_vxc, input=input, ks_env=ks_env, rho=rho) 536 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 537 538 ! Compute XC matrix in AO basis 539 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, & 540 vxc_rho=v_rspace, vxc_tau=tau_rspace, exc=exc) 541 542 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace, matrix_vxc=matrix_vxc) 543 544 CALL cp_fm_set_all(fm_XC_ao, 0.0_dp) 545 CALL copy_dbcsr_to_fm(matrix=matrix_vxc(1)%matrix, fm=fm_XC_ao) 546 CALL cp_fm_set_all(fm_XC_mo, 0.0_dp) 547 548 DO i = 1, SIZE(v_rspace) 549 CALL pw_release(v_rspace(i)%pw) 550 ENDDO 551 DEALLOCATE (v_rspace) 552 553 ! First index 554 CALL cp_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, & 555 mo_coeff, fm_XC_ao, 0.0_dp, fm_XC_ao_mo) 556 557 ! Second index 558 CALL cp_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, & 559 fm_XC_ao_mo, mo_coeff, 1.0_dp, fm_XC_mo) 560 561 ! If open shell 562 IF (alpha_beta) THEN 563 CALL cp_fm_set_all(fm_XC_ao, 0.0_dp) 564 CALL copy_dbcsr_to_fm(matrix=matrix_vxc(2)%matrix, fm=fm_XC_ao) 565 CALL cp_fm_set_all(fm_XC_mo_beta, 0.0_dp) 566 567 ! First index 568 CALL cp_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, & 569 mo_coeff_beta, fm_XC_ao, 0.0_dp, fm_XC_ao_mo) 570 571 ! Second index 572 CALL cp_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, & 573 fm_XC_ao_mo, mo_coeff_beta, 1.0_dp, fm_XC_mo_beta) 574 575 ENDIF 576 577 DO i = 1, SIZE(matrix_vxc) 578 CALL dbcsr_release(matrix_vxc(i)%matrix) 579 DEALLOCATE (matrix_vxc(i)%matrix) 580 ENDDO 581 DEALLOCATE (matrix_vxc) 582 583 CALL timestop(handle) 584 585 END SUBROUTINE xc_contribution 586 587! ************************************************************************************************** 588!> \brief ... 589!> \param fm_F_mo ... 590!> \param eigenval ... 591!> \param dimen ... 592!> \param homo ... 593!> \param para_env ... 594!> \param blacs_env ... 595!> \param rse_corr ... 596!> \param fm_F_mo_beta ... 597!> \param eigenval_beta ... 598!> \param homo_beta ... 599! ************************************************************************************************** 600 SUBROUTINE non_diag_rse(fm_F_mo, eigenval, dimen, homo, para_env, & 601 blacs_env, rse_corr, fm_F_mo_beta, eigenval_beta, homo_beta) 602 TYPE(cp_fm_type), POINTER :: fm_F_mo 603 REAL(KIND=dp), DIMENSION(:) :: Eigenval 604 INTEGER :: dimen, homo 605 TYPE(cp_para_env_type), POINTER :: para_env 606 TYPE(cp_blacs_env_type), POINTER :: blacs_env 607 REAL(KIND=dp) :: rse_corr 608 TYPE(cp_fm_type), OPTIONAL, POINTER :: fm_F_mo_Beta 609 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: Eigenval_beta 610 INTEGER, OPTIONAL :: homo_beta 611 612 CHARACTER(LEN=*), PARAMETER :: routineN = 'non_diag_rse', routineP = moduleN//':'//routineN 613 614 INTEGER :: handle, i_global, iiB, j_global, jjB, & 615 ncol_local, nrow_local, virtual, & 616 virtual_beta 617 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 618 LOGICAL :: alpha_beta 619 REAL(KIND=dp) :: rse_corr_beta 620 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_o, eig_semi_can, eig_v 621 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 622 TYPE(cp_fm_type), POINTER :: fm_F_oo, fm_F_ov, fm_F_vv, fm_O, fm_tmp, & 623 fm_U 624 625 CALL timeset(routineN, handle) 626 627 alpha_beta = .FALSE. 628 IF (PRESENT(eigenval_beta) .AND. PRESENT(homo_beta)) alpha_beta = .TRUE. 629 630 ! Add eigenvalues on the diagonal 631 CALL cp_fm_get_info(matrix=fm_F_mo, & 632 nrow_local=nrow_local, & 633 ncol_local=ncol_local, & 634 row_indices=row_indices, & 635 col_indices=col_indices) 636 637!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 638!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_mo, eigenval) 639 DO jjB = 1, ncol_local 640 j_global = col_indices(jjB) 641 DO iiB = 1, nrow_local 642 i_global = row_indices(iiB) 643 IF (i_global .EQ. j_global) THEN 644 fm_F_mo%local_data(iib, jjb) = & 645 fm_F_mo%local_data(iib, jjb) + eigenval(i_global) 646 ENDIF 647 END DO 648 END DO 649!$OMP END PARALLEL DO 650 651 IF (alpha_beta) THEN 652 ! Add eigenvalues on the diagonal 653 CALL cp_fm_get_info(matrix=fm_F_mo_beta, & 654 nrow_local=nrow_local, & 655 ncol_local=ncol_local, & 656 row_indices=row_indices, & 657 col_indices=col_indices) 658 659!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 660!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_mo_beta,eigenval_beta) 661 DO jjB = 1, ncol_local 662 j_global = col_indices(jjB) 663 DO iiB = 1, nrow_local 664 i_global = row_indices(iiB) 665 IF (i_global .EQ. j_global) fm_F_mo_beta%local_data(iib, jjb) = & 666 fm_F_mo_beta%local_data(iib, jjb) + eigenval_beta(i_global) 667 END DO 668 END DO 669!$OMP END PARALLEL DO 670 ENDIF 671 672 ! Create the occupied-occupied and virtual-virtual blocks, eigenvectors 673 NULLIFY (fm_F_oo, fm_O, fm_struct_tmp) 674 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 675 nrow_global=homo, ncol_global=homo) 676 CALL cp_fm_create(fm_F_oo, fm_struct_tmp, name="F_oo") 677 CALL cp_fm_create(fm_O, fm_struct_tmp, name="O") 678 CALL cp_fm_set_all(fm_F_oo, 0.0_dp) 679 CALL cp_fm_set_all(fm_O, 0.0_dp) 680 CALL cp_fm_struct_release(fm_struct_tmp) 681 682 CALL cp_fm_to_fm_submat(msource=fm_F_mo, mtarget=fm_F_oo, & 683 nrow=homo, ncol=homo, & 684 s_firstrow=1, s_firstcol=1, & 685 t_firstrow=1, t_firstcol=1) 686 687 virtual = dimen - homo 688 NULLIFY (fm_F_vv, fm_U, fm_struct_tmp) 689 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 690 nrow_global=virtual, ncol_global=virtual) 691 CALL cp_fm_create(fm_F_vv, fm_struct_tmp, name="F_vv") 692 CALL cp_fm_create(fm_U, fm_struct_tmp, name="U") 693 CALL cp_fm_set_all(fm_F_vv, 0.0_dp) 694 CALL cp_fm_set_all(fm_U, 0.0_dp) 695 CALL cp_fm_struct_release(fm_struct_tmp) 696 697 CALL cp_fm_to_fm_submat(msource=fm_F_mo, mtarget=fm_F_vv, & 698 nrow=virtual, ncol=virtual, & 699 s_firstrow=homo + 1, s_firstcol=homo + 1, & 700 t_firstrow=1, t_firstcol=1) 701 702 ! Diagonalize occupied-occupied and virtual-virtual matrices 703 704 ALLOCATE (eig_o(homo)) 705 ALLOCATE (eig_v(virtual)) 706 eig_v = 0.0_dp 707 eig_o = 0.0_dp 708 CALL choose_eigv_solver(fm_F_oo, fm_O, eig_o) 709 CALL choose_eigv_solver(fm_F_vv, fm_U, eig_v) 710 711 ! Collect the eigenvalues to one array 712 ALLOCATE (eig_semi_can(dimen)) 713 eig_semi_can = 0.0_dp 714 eig_semi_can(1:homo) = eig_o(:) 715 eig_semi_can(homo + 1:dimen) = eig_v(:) 716 717 ! Create occupied-virtual block 718 NULLIFY (fm_F_ov, fm_tmp, fm_struct_tmp) 719 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 720 nrow_global=homo, ncol_global=virtual) 721 CALL cp_fm_create(fm_F_ov, fm_struct_tmp, name="F_ov") 722 CALL cp_fm_create(fm_tmp, fm_struct_tmp, name="tmp") 723 CALL cp_fm_set_all(fm_F_ov, 0.0_dp) 724 CALL cp_fm_set_all(fm_tmp, 0.0_dp) 725 CALL cp_fm_struct_release(fm_struct_tmp) 726 727 CALL cp_fm_to_fm_submat(msource=fm_F_mo, mtarget=fm_F_ov, & 728 nrow=homo, ncol=virtual, & 729 s_firstrow=1, s_firstcol=homo + 1, & 730 t_firstrow=1, t_firstcol=1) 731 732 CALL cp_fm_get_info(matrix=fm_F_ov, & 733 nrow_local=nrow_local, & 734 ncol_local=ncol_local, & 735 row_indices=row_indices, & 736 col_indices=col_indices) 737 738 CALL cp_gemm(transa='N', transb='N', m=homo, n=virtual, k=homo, alpha=1.0_dp, & 739 matrix_a=fm_O, matrix_b=fm_F_ov, beta=0.0_dp, matrix_c=fm_tmp) 740 741 CALL cp_fm_set_all(fm_F_ov, 0.0_dp) 742 743 CALL cp_gemm(transa='N', transb='N', m=homo, n=virtual, k=virtual, alpha=1.0_dp, & 744 matrix_a=fm_tmp, matrix_b=fm_U, beta=0.0_dp, matrix_c=fm_F_ov) 745 746 ! Compute the correction 747 CALL cp_fm_get_info(matrix=fm_F_ov, & 748 nrow_local=nrow_local, & 749 ncol_local=ncol_local, & 750 row_indices=row_indices, & 751 col_indices=col_indices) 752 753 rse_corr = 0.0_dp 754!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 755!$OMP REDUCTION(+:rse_corr) & 756!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_ov,eig_semi_can,homo) 757 DO jjB = 1, ncol_local 758 j_global = col_indices(jjB) 759 DO iiB = 1, nrow_local 760 i_global = row_indices(iiB) 761 rse_corr = rse_corr + fm_F_ov%local_data(iib, jjb)**2.0_dp/ & 762 (eig_semi_can(i_global) - eig_semi_can(j_global + homo)) 763 END DO 764 END DO 765!$OMP END PARALLEL DO 766 CALL mp_sum(rse_corr, para_env%group) 767 768 ! Release 769 DEALLOCATE (eig_semi_can) 770 DEALLOCATE (eig_o) 771 DEALLOCATE (eig_v) 772 773 CALL cp_fm_release(fm_F_ov) 774 CALL cp_fm_release(fm_F_oo) 775 CALL cp_fm_release(fm_F_vv) 776 CALL cp_fm_release(fm_U) 777 CALL cp_fm_release(fm_O) 778 CALL cp_fm_release(fm_tmp) 779 780 ! Beta spin 781 IF (alpha_beta) THEN 782 ! Create the occupied-occupied and virtual-virtual blocks, eigenvectors 783 NULLIFY (fm_F_oo, fm_O, fm_struct_tmp) 784 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 785 nrow_global=homo_beta, ncol_global=homo_beta) 786 CALL cp_fm_create(fm_F_oo, fm_struct_tmp, name="F_oo") 787 CALL cp_fm_create(fm_O, fm_struct_tmp, name="O") 788 CALL cp_fm_set_all(fm_F_oo, 0.0_dp) 789 CALL cp_fm_set_all(fm_O, 0.0_dp) 790 CALL cp_fm_struct_release(fm_struct_tmp) 791 792 CALL cp_fm_to_fm_submat(msource=fm_F_mo_beta, mtarget=fm_F_oo, & 793 nrow=homo_beta, ncol=homo_beta, & 794 s_firstrow=1, s_firstcol=1, & 795 t_firstrow=1, t_firstcol=1) 796 virtual_beta = dimen - homo_beta 797 NULLIFY (fm_F_vv, fm_U, fm_struct_tmp) 798 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 799 nrow_global=virtual_beta, ncol_global=virtual_beta) 800 CALL cp_fm_create(fm_F_vv, fm_struct_tmp, name="F_vv") 801 CALL cp_fm_create(fm_U, fm_struct_tmp, name="U") 802 CALL cp_fm_set_all(fm_F_vv, 0.0_dp) 803 CALL cp_fm_set_all(fm_U, 0.0_dp) 804 CALL cp_fm_struct_release(fm_struct_tmp) 805 806 CALL cp_fm_to_fm_submat(msource=fm_F_mo_beta, mtarget=fm_F_vv, & 807 nrow=virtual_beta, ncol=virtual_beta, & 808 s_firstrow=homo_beta + 1, s_firstcol=homo_beta + 1, & 809 t_firstrow=1, t_firstcol=1) 810 811 ! Diagonalize occupied-occupied and virtual-virtual matrices 812 ALLOCATE (eig_o(homo_beta)) 813 ALLOCATE (eig_v(virtual_beta)) 814 eig_v = 0.0_dp 815 eig_o = 0.0_dp 816 CALL choose_eigv_solver(fm_F_oo, fm_O, eig_o) 817 CALL choose_eigv_solver(fm_F_vv, fm_U, eig_v) 818 819 ! Collect the eigenvalues to one array 820 ALLOCATE (eig_semi_can(dimen)) 821 eig_semi_can = 0.0_dp 822 eig_semi_can(1:homo_beta) = eig_o(:) 823 eig_semi_can(homo_beta + 1:dimen) = eig_v(:) 824 825 ! Create occupied-virtual block 826 NULLIFY (fm_F_ov, fm_tmp, fm_struct_tmp) 827 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 828 nrow_global=homo_beta, ncol_global=virtual_beta) 829 CALL cp_fm_create(fm_F_ov, fm_struct_tmp, name="F_ov") 830 CALL cp_fm_create(fm_tmp, fm_struct_tmp, name="tmp") 831 CALL cp_fm_set_all(fm_F_ov, 0.0_dp) 832 CALL cp_fm_set_all(fm_tmp, 0.0_dp) 833 CALL cp_fm_struct_release(fm_struct_tmp) 834 835 CALL cp_fm_to_fm_submat(msource=fm_F_mo_beta, mtarget=fm_F_ov, & 836 nrow=homo_beta, ncol=virtual_beta, & 837 s_firstrow=1, s_firstcol=homo_beta + 1, & 838 t_firstrow=1, t_firstcol=1) 839 840 CALL cp_gemm(transa='N', transb='N', m=homo_beta, n=virtual_beta, k=homo_beta, alpha=1.0_dp, & 841 matrix_a=fm_O, matrix_b=fm_F_ov, beta=0.0_dp, matrix_c=fm_tmp) 842 843 CALL cp_gemm(transa='N', transb='N', m=homo_beta, n=virtual_beta, k=virtual_beta, alpha=1.0_dp, & 844 matrix_a=fm_tmp, matrix_b=fm_U, beta=0.0_dp, matrix_c=fm_F_ov) 845 846 ! Compute the correction 847 CALL cp_fm_get_info(matrix=fm_F_ov, & 848 nrow_local=nrow_local, & 849 ncol_local=ncol_local, & 850 row_indices=row_indices, & 851 col_indices=col_indices) 852 rse_corr_beta = 0.0_dp 853!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) & 854!$OMP REDUCTION(+:rse_corr_beta) & 855!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_ov,eig_semi_can,homo_beta) 856 DO jjB = 1, ncol_local 857 j_global = col_indices(jjB) 858 DO iiB = 1, nrow_local 859 i_global = row_indices(iiB) 860 rse_corr_beta = rse_corr_beta + fm_F_ov%local_data(iib, jjb)**2.0_dp/ & 861 (eig_semi_can(i_global) - eig_semi_can(j_global + homo_beta)) 862 END DO 863 END DO 864!$OMP END PARALLEL DO 865 CALL mp_sum(rse_corr_beta, para_env%group) 866 867 ! Release 868 DEALLOCATE (eig_semi_can) 869 DEALLOCATE (eig_o) 870 DEALLOCATE (eig_v) 871 872 CALL cp_fm_release(fm_F_ov) 873 CALL cp_fm_release(fm_F_oo) 874 CALL cp_fm_release(fm_F_vv) 875 CALL cp_fm_release(fm_U) 876 CALL cp_fm_release(fm_O) 877 CALL cp_fm_release(fm_tmp) 878 879 rse_corr = 0.5_dp*(rse_corr + rse_corr_beta) 880 881 ENDIF 882 883 CALL timestop(handle) 884 885 END SUBROUTINE non_diag_rse 886 887END MODULE rpa_rse 888