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 calculate 2c- and 3c-integrals for RI with GPW 8!> \par History 9!> 07.2019 separated from mp2_integrals.F [Frederick Stein] 10! ************************************************************************************************** 11MODULE mp2_eri_gpw 12 USE atomic_kind_types, ONLY: atomic_kind_type 13 USE basis_set_types, ONLY: gto_basis_set_type 14 USE cell_types, ONLY: cell_type,& 15 pbc 16 USE cp_control_types, ONLY: dft_control_type 17 USE cp_fm_types, ONLY: cp_fm_type 18 USE cp_para_types, ONLY: cp_para_env_type 19 USE dbcsr_api, ONLY: dbcsr_p_type,& 20 dbcsr_set 21 USE gaussian_gridlevels, ONLY: gaussian_gridlevel 22 USE input_constants, ONLY: do_potential_coulomb,& 23 do_potential_long 24 USE kinds, ONLY: dp 25 USE message_passing, ONLY: mp_sum 26 USE orbital_pointers, ONLY: ncoset 27 USE particle_types, ONLY: particle_type 28 USE pw_env_methods, ONLY: pw_env_create,& 29 pw_env_rebuild 30 USE pw_env_types, ONLY: pw_env_get,& 31 pw_env_release,& 32 pw_env_type 33 USE pw_methods, ONLY: pw_gauss_damp,& 34 pw_scale,& 35 pw_transfer 36 USE pw_poisson_methods, ONLY: pw_poisson_solve 37 USE pw_poisson_types, ONLY: pw_poisson_type 38 USE pw_pool_types, ONLY: pw_pool_create_pw,& 39 pw_pool_give_back_pw,& 40 pw_pool_type 41 USE pw_types, ONLY: COMPLEXDATA1D,& 42 REALDATA3D,& 43 REALSPACE,& 44 RECIPROCALSPACE,& 45 pw_p_type 46 USE qs_collocate_density, ONLY: calculate_wavefunction 47 USE qs_environment_types, ONLY: get_qs_env,& 48 qs_environment_type 49 USE qs_integrate_potential, ONLY: integrate_pgf_product_rspace,& 50 integrate_v_rspace 51 USE qs_kind_types, ONLY: get_qs_kind,& 52 qs_kind_type 53 USE qs_ks_types, ONLY: qs_ks_env_type 54 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 55 USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,& 56 realspace_grid_p_type,& 57 rs_grid_release,& 58 rs_grid_retain 59 USE rs_pw_interface, ONLY: potential_pw2rs 60 USE task_list_methods, ONLY: generate_qs_task_list 61 USE task_list_types, ONLY: allocate_task_list,& 62 deallocate_task_list,& 63 task_list_type 64 65!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num 66#include "./base/base_uses.f90" 67 68 IMPLICIT NONE 69 70 PRIVATE 71 72 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_eri_gpw' 73 74 PUBLIC :: mp2_eri_2c_integrate_gpw, mp2_eri_3c_integrate_gpw, calc_potential_gpw, cleanup_gpw, prepare_gpw 75 76CONTAINS 77 78! ************************************************************************************************** 79!> \brief ... 80!> \param mo_coeff ... 81!> \param psi_L ... 82!> \param rho_g ... 83!> \param atomic_kind_set ... 84!> \param qs_kind_set ... 85!> \param cell ... 86!> \param dft_control ... 87!> \param particle_set ... 88!> \param pw_env_sub ... 89!> \param external_vector ... 90!> \param poisson_env ... 91!> \param rho_r ... 92!> \param pot_g ... 93!> \param ri_metric ... 94!> \param omega_metric ... 95!> \param mat_munu ... 96!> \param qs_env ... 97!> \param task_list_sub ... 98! ************************************************************************************************** 99 SUBROUTINE mp2_eri_3c_integrate_gpw(mo_coeff, psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, & 100 pw_env_sub, external_vector, poisson_env, rho_r, pot_g, ri_metric, omega_metric, mat_munu, & 101 qs_env, task_list_sub) 102 103 TYPE(cp_fm_type), POINTER :: mo_coeff 104 TYPE(pw_p_type), INTENT(INOUT) :: psi_L, rho_g 105 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 106 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 107 TYPE(cell_type), POINTER :: cell 108 TYPE(dft_control_type), POINTER :: dft_control 109 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 110 TYPE(pw_env_type), POINTER :: pw_env_sub 111 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: external_vector 112 TYPE(pw_poisson_type), POINTER :: poisson_env 113 TYPE(pw_p_type), INTENT(INOUT) :: rho_r, pot_g 114 INTEGER, INTENT(IN) :: ri_metric 115 REAL(KIND=dp), INTENT(IN) :: omega_metric 116 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu 117 TYPE(qs_environment_type), POINTER :: qs_env 118 TYPE(task_list_type), POINTER :: task_list_sub 119 120 CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_eri_3c_integrate_gpw', & 121 routineP = moduleN//':'//routineN 122 123 INTEGER :: handle 124 125 CALL timeset(routineN, handle) 126 127 ! pseudo psi_L 128 CALL calculate_wavefunction(mo_coeff, 1, psi_L, rho_g, atomic_kind_set, & 129 qs_kind_set, cell, dft_control, particle_set, pw_env_sub, & 130 basis_type="RI_AUX", & 131 external_vector=external_vector) 132 133 rho_r%pw%cr3d = psi_L%pw%cr3d 134 135 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, ri_metric, omega_metric) 136 137 ! and finally (K|mu nu) 138 CALL dbcsr_set(mat_munu%matrix, 0.0_dp) 139 CALL integrate_v_rspace(rho_r, hmat=mat_munu, qs_env=qs_env, & 140 calculate_forces=.FALSE., compute_tau=.FALSE., gapw=.FALSE., & 141 pw_env_external=pw_env_sub, task_list_external=task_list_sub) 142 143 CALL timestop(handle) 144 145 END SUBROUTINE mp2_eri_3c_integrate_gpw 146 147! ************************************************************************************************** 148!> \brief ... 149!> \param qs_env ... 150!> \param para_env_sub ... 151!> \param dimen_RI ... 152!> \param mo_coeff ... 153!> \param my_group_L_start ... 154!> \param my_group_L_end ... 155!> \param ri_metric ... 156!> \param natom ... 157!> \param omega ... 158!> \param sab_orb_sub ... 159!> \param L_local_col ... 160!> \param kind_of ... 161! ************************************************************************************************** 162 SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, dimen_RI, mo_coeff, my_group_L_start, my_group_L_end, ri_metric, & 163 natom, omega, sab_orb_sub, L_local_col, kind_of) 164 165 TYPE(qs_environment_type), POINTER :: qs_env 166 TYPE(cp_para_env_type), POINTER :: para_env_sub 167 INTEGER, INTENT(IN) :: dimen_RI 168 TYPE(cp_fm_type), POINTER :: mo_coeff 169 INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, & 170 ri_metric, natom 171 REAL(KIND=dp), INTENT(IN) :: omega 172 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 173 POINTER :: sab_orb_sub 174 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: L_local_col 175 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of 176 177 CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_eri_2c_integrate_gpw', & 178 routineP = moduleN//':'//routineN 179 180 INTEGER :: dir, handle, handle2, i, i_counter, iatom, igrid_level, ikind, ipgf, iset, lb(3), & 181 LLL, location(3), na1, na2, ncoa, nseta, offset, sgfa, tp(3), ub(3) 182 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa 183 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa 184 LOGICAL :: map_it_here 185 REAL(KIND=dp) :: cutoff_old, rab2, relative_cutoff_old 186 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector 187 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab 188 REAL(KIND=dp), DIMENSION(3) :: ra, rab 189 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a 190 REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, rpgfa, sphi_a, zeta 191 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 192 TYPE(cell_type), POINTER :: cell 193 TYPE(dft_control_type), POINTER :: dft_control 194 TYPE(gto_basis_set_type), POINTER :: basis_set_a 195 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 196 TYPE(pw_env_type), POINTER :: pw_env_sub 197 TYPE(pw_p_type) :: pot_g, psi_L, rho_g, rho_r 198 TYPE(pw_poisson_type), POINTER :: poisson_env 199 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 200 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 201 TYPE(realspace_grid_desc_p_type), DIMENSION(:), & 202 POINTER :: rs_descs 203 TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_v 204 TYPE(task_list_type), POINTER :: task_list_sub 205 206 CALL timeset(routineN, handle) 207 208 CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, & 209 auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub) 210 211 ALLOCATE (wf_vector(dimen_RI)) 212 213 CALL get_qs_env(qs_env, cell=cell, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, particle_set=particle_set) 214 215 L_local_col = 0.0_dp 216 217 i_counter = 0 218 DO LLL = my_group_L_start, my_group_L_end 219 i_counter = i_counter + 1 220 221 wf_vector = 0.0_dp 222 wf_vector(LLL) = 1.0_dp 223 224 ! pseudo psi_L 225 CALL calculate_wavefunction(mo_coeff, 1, psi_L, rho_g, atomic_kind_set, & 226 qs_kind_set, cell, dft_control, particle_set, pw_env_sub, & 227 basis_type="RI_AUX", & 228 external_vector=wf_vector) 229 230 CALL timeset(routineN//"_pot_lm", handle2) 231 rho_r%pw%cr3d = psi_L%pw%cr3d 232 233 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, ri_metric, omega) 234 235 NULLIFY (rs_v) 236 NULLIFY (rs_descs) 237 CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v) 238 DO i = 1, SIZE(rs_v) 239 CALL rs_grid_retain(rs_v(i)%rs_grid) 240 END DO 241 CALL potential_pw2rs(rs_v, rho_r, pw_env_sub) 242 243 CALL timestop(handle2) 244 245 ! integrate the little bastards 246 offset = 0 247 DO iatom = 1, natom 248 ikind = kind_of(iatom) 249 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX") 250 251 first_sgfa => basis_set_a%first_sgf 252 la_max => basis_set_a%lmax 253 la_min => basis_set_a%lmin 254 npgfa => basis_set_a%npgf 255 nseta = basis_set_a%nset 256 nsgfa => basis_set_a%nsgf_set 257 rpgfa => basis_set_a%pgf_radius 258 set_radius_a => basis_set_a%set_radius 259 sphi_a => basis_set_a%sphi 260 zeta => basis_set_a%zet 261 262 ra(:) = pbc(particle_set(iatom)%r, cell) 263 rab = 0.0_dp 264 rab2 = 0.0_dp 265 266 DO iset = 1, nseta 267 ncoa = npgfa(iset)*ncoset(la_max(iset)) 268 sgfa = first_sgfa(1, iset) 269 270 ALLOCATE (I_tmp2(ncoa, 1)) 271 I_tmp2 = 0.0_dp 272 ALLOCATE (I_ab(nsgfa(iset), 1)) 273 I_ab = 0.0_dp 274 275 igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, MINVAL(zeta(:, iset))) 276 map_it_here = .FALSE. 277 IF (.NOT. ALL(rs_v(igrid_level)%rs_grid%desc%perd == 1)) THEN 278 DO dir = 1, 3 279 ! bounds of local grid (i.e. removing the 'wings'), if periodic 280 tp(dir) = FLOOR(DOT_PRODUCT(cell%h_inv(dir, :), ra)*rs_v(igrid_level)%rs_grid%desc%npts(dir)) 281 tp(dir) = MODULO(tp(dir), rs_v(igrid_level)%rs_grid%desc%npts(dir)) 282 IF (rs_v(igrid_level)%rs_grid%desc%perd(dir) .NE. 1) THEN 283 lb(dir) = rs_v(igrid_level)%rs_grid%lb_local(dir) + rs_v(igrid_level)%rs_grid%desc%border 284 ub(dir) = rs_v(igrid_level)%rs_grid%ub_local(dir) - rs_v(igrid_level)%rs_grid%desc%border 285 ELSE 286 lb(dir) = rs_v(igrid_level)%rs_grid%lb_local(dir) 287 ub(dir) = rs_v(igrid_level)%rs_grid%ub_local(dir) 288 ENDIF 289 ! distributed grid, only map if it is local to the grid 290 location(dir) = tp(dir) + rs_v(igrid_level)%rs_grid%desc%lb(dir) 291 ENDDO 292 IF (lb(1) <= location(1) .AND. location(1) <= ub(1) .AND. & 293 lb(2) <= location(2) .AND. location(2) <= ub(2) .AND. & 294 lb(3) <= location(3) .AND. location(3) <= ub(3)) THEN 295 map_it_here = .TRUE. 296 ENDIF 297 ELSE 298 ! not distributed, just a round-robin distribution over the full set of CPUs 299 IF (MODULO(offset, para_env_sub%num_pe) == para_env_sub%mepos) map_it_here = .TRUE. 300 ENDIF 301 302 offset = offset + nsgfa(iset) 303 304 IF (map_it_here) THEN 305 DO ipgf = 1, npgfa(iset) 306 sgfa = first_sgfa(1, iset) 307 308 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 309 na2 = ipgf*ncoset(la_max(iset)) 310 igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset)) 311 312 CALL integrate_pgf_product_rspace( & 313 la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), & 314 lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, & 315 ra=ra, rab=rab, rab2=rab2, & 316 rsgrid=rs_v(igrid_level)%rs_grid, & 317 cell=cell, & 318 cube_info=pw_env_sub%cube_info(igrid_level), & 319 hab=I_tmp2, & 320 o1=na1 - 1, & 321 o2=0, & 322 map_consistent=.TRUE., & 323 eps_gvg_rspace=dft_control%qs_control%eps_gvg_rspace, & 324 calculate_forces=.FALSE.) 325 END DO 326 327 CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, & 328 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 329 I_tmp2(1, 1), SIZE(I_tmp2, 1), & 330 1.0_dp, I_ab(1, 1), SIZE(I_ab, 1)) 331 332 L_local_col(offset - nsgfa(iset) + 1:offset, i_counter) = I_ab(1:nsgfa(iset), 1) 333 END IF 334 335 DEALLOCATE (I_tmp2) 336 DEALLOCATE (I_ab) 337 338 END DO 339 END DO 340 341 DO i = 1, SIZE(rs_v) 342 CALL rs_grid_release(rs_v(i)%rs_grid) 343 END DO 344 345 END DO 346 347 DEALLOCATE (wf_vector) 348 349 CALL mp_sum(L_local_col, para_env_sub%group) 350 351 CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, pw_env_sub, & 352 task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L) 353 354 CALL timestop(handle) 355 356 END SUBROUTINE 357 358! ************************************************************************************************** 359!> \brief ... 360!> \param qs_env ... 361!> \param dft_control ... 362!> \param e_cutoff_old ... 363!> \param cutoff_old ... 364!> \param relative_cutoff_old ... 365!> \param para_env_sub ... 366!> \param pw_env_sub ... 367!> \param auxbas_pw_pool ... 368!> \param poisson_env ... 369!> \param task_list_sub ... 370!> \param rho_r ... 371!> \param rho_g ... 372!> \param pot_g ... 373!> \param psi_L ... 374!> \param sab_orb_sub ... 375! ************************************************************************************************** 376 SUBROUTINE prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, & 377 auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub) 378 TYPE(qs_environment_type), POINTER :: qs_env 379 TYPE(dft_control_type), POINTER :: dft_control 380 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 381 INTENT(OUT) :: e_cutoff_old 382 REAL(KIND=dp), INTENT(OUT) :: cutoff_old, relative_cutoff_old 383 TYPE(cp_para_env_type), POINTER :: para_env_sub 384 TYPE(pw_env_type), POINTER :: pw_env_sub 385 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 386 TYPE(pw_poisson_type), POINTER :: poisson_env 387 TYPE(task_list_type), POINTER :: task_list_sub 388 TYPE(pw_p_type), INTENT(OUT) :: rho_r, rho_g, pot_g, psi_L 389 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 390 INTENT(IN), POINTER :: sab_orb_sub 391 392 CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_gpw', routineP = moduleN//':'//routineN 393 394 INTEGER :: handle, i_multigrid, n_multigrid 395 LOGICAL :: skip_load_balance_distributed 396 REAL(KIND=dp) :: progression_factor 397 TYPE(qs_ks_env_type), POINTER :: ks_env 398 399 CALL timeset(routineN, handle) 400 401 CALL get_qs_env(qs_env, dft_control=dft_control, ks_env=ks_env) 402 403 ! hack hack hack XXXXXXXXXXXXXXX rebuilds the pw_en with the new cutoffs 404 progression_factor = dft_control%qs_control%progression_factor 405 n_multigrid = SIZE(dft_control%qs_control%e_cutoff) 406 ALLOCATE (e_cutoff_old(n_multigrid)) 407 e_cutoff_old(:) = dft_control%qs_control%e_cutoff 408 cutoff_old = dft_control%qs_control%cutoff 409 410 dft_control%qs_control%cutoff = qs_env%mp2_env%mp2_gpw%cutoff*0.5_dp 411 dft_control%qs_control%e_cutoff(1) = dft_control%qs_control%cutoff 412 DO i_multigrid = 2, n_multigrid 413 dft_control%qs_control%e_cutoff(i_multigrid) = dft_control%qs_control%e_cutoff(i_multigrid - 1) & 414 /progression_factor 415 END DO 416 417 relative_cutoff_old = dft_control%qs_control%relative_cutoff 418 dft_control%qs_control%relative_cutoff = qs_env%mp2_env%mp2_gpw%relative_cutoff*0.5_dp 419 420 ! a pw_env 421 NULLIFY (pw_env_sub) 422 CALL pw_env_create(pw_env_sub) 423 CALL pw_env_rebuild(pw_env_sub, qs_env, para_env_sub) 424 425 CALL pw_env_get(pw_env_sub, auxbas_pw_pool=auxbas_pw_pool, & 426 poisson_env=poisson_env) 427 ! hack hack hack XXXXXXXXXXXXXXX 428 429 ! now we need a task list, hard code skip_load_balance_distributed 430 NULLIFY (task_list_sub) 431 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed 432 CALL allocate_task_list(task_list_sub) 433 CALL generate_qs_task_list(ks_env, task_list_sub, & 434 reorder_rs_grid_ranks=.TRUE., soft_valid=.FALSE., & 435 skip_load_balance_distributed=skip_load_balance_distributed, & 436 pw_env_external=pw_env_sub, sab_orb_external=sab_orb_sub) 437 438 ! get some of the grids ready 439 NULLIFY (rho_r%pw, rho_g%pw, pot_g%pw, psi_L%pw) 440 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r%pw, & 441 use_data=REALDATA3D, & 442 in_space=REALSPACE) 443 CALL pw_pool_create_pw(auxbas_pw_pool, rho_g%pw, & 444 use_data=COMPLEXDATA1D, & 445 in_space=RECIPROCALSPACE) 446 CALL pw_pool_create_pw(auxbas_pw_pool, pot_g%pw, & 447 use_data=COMPLEXDATA1D, & 448 in_space=RECIPROCALSPACE) 449 CALL pw_pool_create_pw(auxbas_pw_pool, psi_L%pw, & 450 use_data=REALDATA3D, & 451 in_space=REALSPACE) 452 453 ! run the FFT once, to set up buffers and to take into account the memory 454 rho_r%pw%cr3d = 0.0D0 455 CALL pw_transfer(rho_r%pw, rho_g%pw) 456 457 CALL timestop(handle) 458 459 END SUBROUTINE prepare_gpw 460 461! ************************************************************************************************** 462!> \brief ... 463!> \param qs_env ... 464!> \param e_cutoff_old ... 465!> \param cutoff_old ... 466!> \param relative_cutoff_old ... 467!> \param pw_env_sub ... 468!> \param task_list_sub ... 469!> \param auxbas_pw_pool ... 470!> \param rho_r ... 471!> \param rho_g ... 472!> \param pot_g ... 473!> \param psi_L ... 474! ************************************************************************************************** 475 SUBROUTINE cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, pw_env_sub, & 476 task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L) 477 TYPE(qs_environment_type), POINTER :: qs_env 478 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 479 INTENT(IN) :: e_cutoff_old 480 REAL(KIND=dp), INTENT(IN) :: cutoff_old, relative_cutoff_old 481 TYPE(pw_env_type), POINTER :: pw_env_sub 482 TYPE(task_list_type), POINTER :: task_list_sub 483 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 484 TYPE(pw_p_type), INTENT(INOUT) :: rho_r, rho_g, pot_g, psi_L 485 486 CHARACTER(LEN=*), PARAMETER :: routineN = 'cleanup_gpw', routineP = moduleN//':'//routineN 487 488 INTEGER :: handle 489 TYPE(dft_control_type), POINTER :: dft_control 490 491 CALL timeset(routineN, handle) 492 493 ! and now free the whole lot 494 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r%pw) 495 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_g%pw) 496 CALL pw_pool_give_back_pw(auxbas_pw_pool, pot_g%pw) 497 CALL pw_pool_give_back_pw(auxbas_pw_pool, psi_L%pw) 498 499 CALL deallocate_task_list(task_list_sub) 500 CALL pw_env_release(pw_env_sub) 501 502 CALL get_qs_env(qs_env, dft_control=dft_control) 503 504 ! restore the initial value of the cutoff 505 dft_control%qs_control%e_cutoff = e_cutoff_old 506 dft_control%qs_control%cutoff = cutoff_old 507 dft_control%qs_control%relative_cutoff = relative_cutoff_old 508 509 CALL timestop(handle) 510 511 END SUBROUTINE cleanup_gpw 512 513! ************************************************************************************************** 514!> \brief ... 515!> \param rho_r ... 516!> \param rho_g ... 517!> \param poisson_env ... 518!> \param pot_g ... 519!> \param potential_type ... 520!> \param omega ... 521! ************************************************************************************************** 522 SUBROUTINE calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_type, omega) 523 TYPE(pw_p_type), INTENT(IN) :: rho_r, rho_g 524 TYPE(pw_poisson_type), POINTER :: poisson_env 525 TYPE(pw_p_type), INTENT(IN) :: pot_g 526 INTEGER, INTENT(IN), OPTIONAL :: potential_type 527 REAL(KIND=dp), INTENT(IN), OPTIONAL :: omega 528 529 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_potential_gpw', & 530 routineP = moduleN//':'//routineN 531 532 INTEGER :: handle, my_potential_type 533 534 CALL timeset(routineN, handle) 535 536 my_potential_type = do_potential_coulomb 537 IF (PRESENT(potential_type)) THEN 538 my_potential_type = potential_type 539 IF (my_potential_type /= do_potential_coulomb .AND. .NOT. PRESENT(omega)) THEN 540 CPABORT("Need omega when longrange potential is requested.") 541 END IF 542 END IF 543 544 ! in case we do Coulomb metric for RI, we need the Coulomb operator, but for RI with the 545 ! overlap metric, we do not need the Coulomb operator 546 IF (my_potential_type == do_potential_coulomb .OR. my_potential_type == do_potential_long) THEN 547 CALL pw_transfer(rho_r%pw, rho_g%pw) 548 CALL pw_poisson_solve(poisson_env, rho_g%pw, vhartree=pot_g%pw) 549 IF (my_potential_type == do_potential_long) CALL pw_gauss_damp(pot_g%pw, omega) 550 CALL pw_transfer(pot_g%pw, rho_r%pw) 551 END IF 552 553 CALL pw_scale(rho_r%pw, rho_r%pw%pw_grid%dvol) 554 CALL timestop(handle) 555 556 END SUBROUTINE calc_potential_gpw 557 558END MODULE mp2_eri_gpw 559