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