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