(git:34ef472)
mp2_eri_gpw.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Routines to calculate 2c- and 3c-integrals for RI with GPW
10 !> \par History
11 !> 07.2019 separated from mp2_integrals.F [Frederick Stein]
12 ! **************************************************************************************************
15  USE atomic_kind_types, ONLY: atomic_kind_type
16  USE basis_set_types, ONLY: gto_basis_set_type
17  USE cell_types, ONLY: cell_type,&
18  pbc
19  USE cp_control_types, ONLY: dft_control_type
20  USE cp_fm_types, ONLY: cp_fm_type
21  USE dbcsr_api, ONLY: dbcsr_p_type,&
22  dbcsr_set
30  USE kinds, ONLY: dp
31  USE libint_2c_3c, ONLY: libint_potential_type
32  USE mathconstants, ONLY: fourpi
33  USE message_passing, ONLY: mp_para_env_type
34  USE orbital_pointers, ONLY: ncoset
35  USE particle_types, ONLY: particle_type
36  USE pw_env_methods, ONLY: pw_env_create,&
38  USE pw_env_types, ONLY: pw_env_get,&
40  pw_env_type
41  USE pw_methods, ONLY: &
42  pw_compl_gauss_damp, pw_copy, pw_derive, pw_gauss_damp, pw_gauss_damp_mix, pw_integral_ab, &
44  pw_scale, pw_transfer, pw_truncated, pw_zero
45  USE pw_poisson_methods, ONLY: pw_poisson_solve
46  USE pw_poisson_types, ONLY: pw_poisson_type
47  USE pw_pool_types, ONLY: pw_pool_type
48  USE pw_types, ONLY: pw_c1d_gs_type,&
49  pw_r3d_rs_type
53  USE qs_environment_types, ONLY: get_qs_env,&
54  qs_environment_type
55  USE qs_force_types, ONLY: qs_force_type
56  USE qs_integrate_potential, ONLY: integrate_pgf_product,&
57  integrate_v_rspace
58  USE qs_kind_types, ONLY: get_qs_kind,&
60  qs_kind_type
61  USE qs_ks_types, ONLY: qs_ks_env_type
62  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
64  realspace_grid_desc_p_type,&
65  realspace_grid_type
70  task_list_type
71 
72 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
73 #include "./base/base_uses.f90"
74 
75  IMPLICIT NONE
76 
77  PRIVATE
78 
79  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_eri_gpw'
80 
84 
85 CONTAINS
86 
87 ! **************************************************************************************************
88 !> \brief ...
89 !> \param mo_coeff ...
90 !> \param psi_L ...
91 !> \param rho_g ...
92 !> \param atomic_kind_set ...
93 !> \param qs_kind_set ...
94 !> \param cell ...
95 !> \param dft_control ...
96 !> \param particle_set ...
97 !> \param pw_env_sub ...
98 !> \param external_vector ...
99 !> \param poisson_env ...
100 !> \param rho_r ...
101 !> \param pot_g ...
102 !> \param potential_parameter ...
103 !> \param mat_munu ...
104 !> \param qs_env ...
105 !> \param task_list_sub ...
106 ! **************************************************************************************************
107  SUBROUTINE mp2_eri_3c_integrate_gpw(mo_coeff, psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
108  pw_env_sub, external_vector, poisson_env, rho_r, pot_g, &
109  potential_parameter, mat_munu, qs_env, task_list_sub)
110 
111  TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
112  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: psi_l
113  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
114  TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
115  POINTER :: atomic_kind_set
116  TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
117  POINTER :: qs_kind_set
118  TYPE(cell_type), INTENT(IN), POINTER :: cell
119  TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
120  TYPE(particle_type), DIMENSION(:), INTENT(IN), &
121  POINTER :: particle_set
122  TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
123  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: external_vector
124  TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
125  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
126  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
127  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
128  TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu
129  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
130  TYPE(task_list_type), INTENT(IN), POINTER :: task_list_sub
131 
132  CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_eri_3c_integrate_gpw'
133 
134  INTEGER :: handle
135 
136  CALL timeset(routinen, handle)
137 
138  ! pseudo psi_L
139  CALL calculate_wavefunction(mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
140  qs_kind_set, cell, dft_control, particle_set, pw_env_sub, &
141  basis_type="RI_AUX", &
142  external_vector=external_vector)
143 
144  CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
145 
146  ! and finally (K|mu nu)
147  CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
148  CALL integrate_v_rspace(rho_r, hmat=mat_munu, qs_env=qs_env, &
149  calculate_forces=.false., compute_tau=.false., gapw=.false., &
150  pw_env_external=pw_env_sub, task_list_external=task_list_sub)
151 
152  CALL timestop(handle)
153 
154  END SUBROUTINE mp2_eri_3c_integrate_gpw
155 
156 ! **************************************************************************************************
157 !> \brief Integrates the potential of an RI function
158 !> \param qs_env ...
159 !> \param para_env_sub ...
160 !> \param my_group_L_start ...
161 !> \param my_group_L_end ...
162 !> \param natom ...
163 !> \param potential_parameter ...
164 !> \param sab_orb_sub ...
165 !> \param L_local_col ...
166 !> \param kind_of ...
167 ! **************************************************************************************************
168  SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, &
169  natom, potential_parameter, sab_orb_sub, L_local_col, kind_of)
170 
171  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
172  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
173  INTEGER, INTENT(IN) :: my_group_l_start, my_group_l_end, natom
174  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
175  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
176  INTENT(IN), POINTER :: sab_orb_sub
177  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: l_local_col
178  INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
179 
180  CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_eri_2c_integrate_gpw'
181 
182  INTEGER :: dir, handle, handle2, i_counter, iatom, igrid_level, ikind, ipgf, iset, lb(3), &
183  lll, location(3), max_nseta, na1, na2, ncoa, nseta, offset, sgfa, tp(3), ub(3)
184  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: offset_2d
185  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
186  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
187  LOGICAL :: map_it_here, use_subpatch
188  REAL(kind=dp) :: cutoff_old, radius, relative_cutoff_old
189  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old
190  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: i_ab
191  REAL(kind=dp), DIMENSION(3) :: ra
192  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a
193  REAL(kind=dp), DIMENSION(:, :), POINTER :: i_tmp2, rpgfa, sphi_a, zeta
194  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
195  TYPE(cell_type), POINTER :: cell
196  TYPE(dft_control_type), POINTER :: dft_control
197  TYPE(gto_basis_set_type), POINTER :: basis_set_a
198  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
199  TYPE(pw_c1d_gs_type) :: pot_g, rho_g
200  TYPE(pw_env_type), POINTER :: pw_env_sub
201  TYPE(pw_poisson_type), POINTER :: poisson_env
202  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
203  TYPE(pw_r3d_rs_type) :: psi_l, rho_r
204  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
205  TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
206  POINTER :: rs_descs
207  TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
208  TYPE(task_list_type), POINTER :: task_list_sub
209 
210  CALL timeset(routinen, handle)
211 
212  CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
213  auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_l, sab_orb_sub)
214 
215  CALL get_qs_env(qs_env, cell=cell, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, particle_set=particle_set)
216 
217  l_local_col = 0.0_dp
218 
219  i_counter = 0
220  DO lll = my_group_l_start, my_group_l_end
221  i_counter = i_counter + 1
222 
223  ! pseudo psi_L
224  CALL collocate_single_gaussian(psi_l, rho_g, atomic_kind_set, &
225  qs_kind_set, cell, dft_control, particle_set, pw_env_sub, &
226  required_function=lll, basis_type="RI_AUX")
227 
228  CALL timeset(routinen//"_pot_lm", handle2)
229 
230  CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
231 
232  NULLIFY (rs_v)
233  NULLIFY (rs_descs)
234  CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
235  CALL potential_pw2rs(rs_v, rho_r, pw_env_sub)
236 
237  CALL timestop(handle2)
238 
239  offset = 0
240  ! Prepare offset ahead of OMP parallel loop
241  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxnset=max_nseta, basis_type="RI_AUX")
242  ALLOCATE (offset_2d(max_nseta, natom))
243 
244  DO iatom = 1, natom
245  ikind = kind_of(iatom)
246  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
247  nseta = basis_set_a%nset
248  nsgfa => basis_set_a%nsgf_set
249  DO iset = 1, nseta
250  offset = offset + nsgfa(iset)
251  offset_2d(iset, iatom) = offset
252  END DO
253  END DO
254 
255  ! integrate the little bastards
256  !$OMP PARALLEL DO DEFAULT(NONE) &
257  !$OMP SHARED(natom, particle_set, cell, pw_env_sub, rs_v, offset_2d, &
258  !$OMP qs_kind_set, ncoset, para_env_sub, dft_control, i_counter, &
259  !$OMP kind_of, l_local_col) &
260  !$OMP PRIVATE(iatom, ikind, basis_set_a, first_sgfa, la_max, la_min, npgfa, &
261  !$OMP nseta, nsgfa, rpgfa, set_radius_a, sphi_a, zeta, &
262  !$OMP ra, iset, ncoa, I_tmp2, I_ab, igrid_level, &
263  !$OMP map_it_here, dir, tp, lb, ub, location, ipgf, &
264  !$OMP sgfa, na1, na2, radius, offset, use_subpatch)
265  DO iatom = 1, natom
266  ikind = kind_of(iatom)
267  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
268 
269  first_sgfa => basis_set_a%first_sgf
270  la_max => basis_set_a%lmax
271  la_min => basis_set_a%lmin
272  npgfa => basis_set_a%npgf
273  nseta = basis_set_a%nset
274  nsgfa => basis_set_a%nsgf_set
275  rpgfa => basis_set_a%pgf_radius
276  set_radius_a => basis_set_a%set_radius
277  sphi_a => basis_set_a%sphi
278  zeta => basis_set_a%zet
279 
280  ra(:) = pbc(particle_set(iatom)%r, cell)
281 
282  DO iset = 1, nseta
283  ncoa = npgfa(iset)*ncoset(la_max(iset))
284  sgfa = first_sgfa(1, iset)
285 
286  ALLOCATE (i_tmp2(ncoa, 1))
287  i_tmp2 = 0.0_dp
288  ALLOCATE (i_ab(nsgfa(iset), 1))
289  i_ab = 0.0_dp
290 
291  offset = offset_2d(iset, iatom)
292 
293  igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, minval(zeta(:, iset)))
294  use_subpatch = .NOT. all(rs_v(igrid_level)%desc%perd == 1)
295 
296  IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, &
297  offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
298  DO ipgf = 1, npgfa(iset)
299  sgfa = first_sgfa(1, iset)
300  na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
301  na2 = ipgf*ncoset(la_max(iset))
302  igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))
303 
304  radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
305  lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
306  zetp=zeta(ipgf, iset), &
307  eps=dft_control%qs_control%eps_gvg_rspace, &
308  prefactor=1.0_dp, cutoff=1.0_dp)
309 
310  CALL integrate_pgf_product( &
311  la_max=la_max(iset), zeta=zeta(ipgf, iset), la_min=la_min(iset), &
312  lb_max=0, zetb=0.0_dp, lb_min=0, &
313  ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
314  rsgrid=rs_v(igrid_level), &
315  hab=i_tmp2, &
316  o1=na1 - 1, &
317  o2=0, &
318  radius=radius, &
319  calculate_forces=.false., &
320  use_subpatch=use_subpatch, &
321  subpatch_pattern=0)
322  END DO
323 
324  CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, &
325  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
326  i_tmp2(1, 1), SIZE(i_tmp2, 1), &
327  1.0_dp, i_ab(1, 1), SIZE(i_ab, 1))
328 
329  l_local_col(offset - nsgfa(iset) + 1:offset, i_counter) = i_ab(1:nsgfa(iset), 1)
330  END IF
331 
332  DEALLOCATE (i_tmp2)
333  DEALLOCATE (i_ab)
334 
335  END DO
336  END DO
337  !$OMP END PARALLEL DO
338  DEALLOCATE (offset_2d)
339 
340  END DO
341 
342  CALL para_env_sub%sum(l_local_col)
343 
344  CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
345  task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_l)
346 
347  CALL timestop(handle)
348 
349  END SUBROUTINE
350 
351 ! **************************************************************************************************
352 !> \brief Integrates the potential of a RI function obtaining the forces and stress tensor
353 !> \param rho_r ...
354 !> \param LLL ...
355 !> \param matrix ...
356 !> \param rho_g ...
357 !> \param atomic_kind_set ...
358 !> \param qs_kind_set ...
359 !> \param particle_set ...
360 !> \param cell ...
361 !> \param pw_env_sub ...
362 !> \param poisson_env ...
363 !> \param pot_g ...
364 !> \param potential_parameter ...
365 !> \param use_virial ...
366 !> \param rho_g_copy ...
367 !> \param dvg ...
368 !> \param kind_of ...
369 !> \param atom_of_kind ...
370 !> \param G_PQ_local ...
371 !> \param force ...
372 !> \param h_stress ...
373 !> \param para_env_sub ...
374 !> \param dft_control ...
375 !> \param psi_L ...
376 !> \param factor ...
377 ! **************************************************************************************************
378  SUBROUTINE integrate_potential_forces_2c(rho_r, LLL, matrix, rho_g, atomic_kind_set, &
379  qs_kind_set, particle_set, cell, pw_env_sub, poisson_env, pot_g, &
380  potential_parameter, use_virial, rho_g_copy, dvg, &
381  kind_of, atom_of_kind, G_PQ_local, force, h_stress, para_env_sub, &
382  dft_control, psi_L, factor)
383 
384  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
385  INTEGER, INTENT(IN) :: lll
386  TYPE(cp_fm_type), INTENT(IN) :: matrix
387  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
388  TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
389  POINTER :: atomic_kind_set
390  TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
391  POINTER :: qs_kind_set
392  TYPE(particle_type), DIMENSION(:), INTENT(IN), &
393  POINTER :: particle_set
394  TYPE(cell_type), INTENT(IN), POINTER :: cell
395  TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
396  TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
397  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
398  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
399  LOGICAL, INTENT(IN) :: use_virial
400  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g_copy, dvg(3)
401  INTEGER, DIMENSION(:), INTENT(IN) :: kind_of, atom_of_kind
402  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: g_pq_local
403  TYPE(qs_force_type), DIMENSION(:), INTENT(IN), &
404  POINTER :: force
405  REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
406  TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
407  TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
408  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: psi_l
409  REAL(kind=dp), INTENT(IN) :: factor
410 
411  CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_potential_forces_2c'
412 
413  INTEGER :: handle, handle2
414 
415  CALL timeset(routinen, handle)
416 
417  ! calculate potential associated to the single aux function
418  CALL timeset(routinen//"_wf_pot", handle2)
419  ! pseudo psi_L
420  CALL pw_zero(rho_r)
421  CALL pw_zero(rho_g)
422  CALL collocate_single_gaussian(rho_r, rho_g, atomic_kind_set, &
423  qs_kind_set, cell, dft_control, particle_set, &
424  pw_env_sub, required_function=lll, basis_type='RI_AUX')
425  IF (use_virial) THEN
426  CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter, dvg)
427  ELSE
428  CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
429  END IF
430  CALL timestop(handle2)
431 
432  IF (use_virial) THEN
433  ! make a copy of the density in G space
434  CALL pw_copy(rho_g, rho_g_copy)
435 
436  ! add the volume contribution to the virial due to
437  ! the (P|Q) integrals, first we put the full gamme_PQ
438  ! pseudo wave-function into grid in order to calculate the
439  ! hartree potential derivatives
440  CALL pw_zero(psi_l)
441  CALL pw_zero(rho_g)
442  CALL calculate_wavefunction(matrix, 1, psi_l, rho_g, atomic_kind_set, &
443  qs_kind_set, cell, dft_control, particle_set, pw_env_sub, &
444  basis_type="RI_AUX", &
445  external_vector=0.5_dp*factor*g_pq_local)
446  ! transfer to reciprocal space and calculate potential
447  CALL calc_potential_gpw(psi_l, rho_g, poisson_env, pot_g, potential_parameter, no_transfer=.true.)
448  ! update virial with volume term (first calculate hartree like energy (diagonal part of the virial))
449  CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
450  END IF
451 
452  ! integrate the potential of the single gaussian and update
453  ! 2-center forces with Gamma_PQ
454  CALL integrate_potential(pw_env_sub, rho_r, kind_of, atom_of_kind, particle_set, qs_kind_set, &
455  -0.25_dp*factor*g_pq_local, cell, force, use_virial, h_stress, para_env_sub, dft_control)
456 
457  CALL timestop(handle)
458  END SUBROUTINE integrate_potential_forces_2c
459 
460 ! **************************************************************************************************
461 !> \brief Takes the precomputed potential of an RI wave-function and determines matrix element and
462 !> gradients with product of Gaussians
463 !> \param mat_munu ...
464 !> \param rho_r ...
465 !> \param matrix_P_munu ...
466 !> \param qs_env ...
467 !> \param pw_env_sub ...
468 !> \param task_list_sub ...
469 ! **************************************************************************************************
470  SUBROUTINE integrate_potential_forces_3c_1c(mat_munu, rho_r, matrix_P_munu, qs_env, pw_env_sub, &
471  task_list_sub)
472 
473  TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu
474  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_r
475  TYPE(dbcsr_p_type), INTENT(IN) :: matrix_p_munu
476  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
477  TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
478  TYPE(task_list_type), INTENT(INOUT), POINTER :: task_list_sub
479 
480  CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_potential_forces_3c_1c'
481 
482  INTEGER :: handle
483 
484  CALL timeset(routinen, handle)
485 
486  ! integrate the potential of the single gaussian and update
487  ! 3-center forces
488  CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
489  CALL integrate_v_rspace(rho_r, hmat=mat_munu, pmat=matrix_p_munu, &
490  qs_env=qs_env, calculate_forces=.true., compute_tau=.false., gapw=.false., &
491  pw_env_external=pw_env_sub, &
492  task_list_external=task_list_sub)
493 
494  CALL timestop(handle)
495  END SUBROUTINE integrate_potential_forces_3c_1c
496 
497 ! **************************************************************************************************
498 !> \brief Integrates potential of two Gaussians to a potential
499 !> \param matrix_P_munu ...
500 !> \param rho_r ...
501 !> \param rho_g ...
502 !> \param task_list_sub ...
503 !> \param pw_env_sub ...
504 !> \param potential_parameter ...
505 !> \param ks_env ...
506 !> \param poisson_env ...
507 !> \param pot_g ...
508 !> \param use_virial ...
509 !> \param rho_g_copy ...
510 !> \param dvg ...
511 !> \param h_stress ...
512 !> \param para_env_sub ...
513 !> \param kind_of ...
514 !> \param atom_of_kind ...
515 !> \param qs_kind_set ...
516 !> \param particle_set ...
517 !> \param cell ...
518 !> \param LLL ...
519 !> \param force ...
520 !> \param dft_control ...
521 ! **************************************************************************************************
522  SUBROUTINE integrate_potential_forces_3c_2c(matrix_P_munu, rho_r, rho_g, task_list_sub, pw_env_sub, &
523  potential_parameter, &
524  ks_env, poisson_env, pot_g, use_virial, rho_g_copy, dvg, &
525  h_stress, para_env_sub, kind_of, atom_of_kind, &
526  qs_kind_set, particle_set, cell, LLL, force, dft_control)
527 
528  TYPE(dbcsr_p_type), INTENT(IN) :: matrix_p_munu
529  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
530  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
531  TYPE(task_list_type), INTENT(IN), POINTER :: task_list_sub
532  TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
533  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
534  TYPE(qs_ks_env_type), INTENT(IN), POINTER :: ks_env
535  TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
536  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
537  LOGICAL, INTENT(IN) :: use_virial
538  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g_copy
539  TYPE(pw_c1d_gs_type), INTENT(IN) :: dvg(3)
540  REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
541  TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
542  INTEGER, DIMENSION(:), INTENT(IN) :: kind_of, atom_of_kind
543  TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
544  POINTER :: qs_kind_set
545  TYPE(particle_type), DIMENSION(:), INTENT(IN), &
546  POINTER :: particle_set
547  TYPE(cell_type), INTENT(IN), POINTER :: cell
548  INTEGER, INTENT(IN) :: lll
549  TYPE(qs_force_type), DIMENSION(:), INTENT(IN), &
550  POINTER :: force
551  TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
552 
553  CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_potential_forces_3c_2c'
554 
555  INTEGER :: atom_a, handle, handle2, iatom, &
556  igrid_level, ikind, iorb, ipgf, iset, &
557  na1, na2, ncoa, nseta, offset, sgfa
558  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
559  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
560  LOGICAL :: map_it_here, skip_shell, use_subpatch
561  REAL(kind=dp) :: radius
562  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: i_ab
563  REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
564  REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
565  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a
566  REAL(kind=dp), DIMENSION(:, :), POINTER :: i_tmp2, pab, rpgfa, sphi_a, zeta
567  TYPE(gto_basis_set_type), POINTER :: basis_set_a
568  TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
569  POINTER :: rs_descs
570  TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
571 
572  CALL timeset(routinen, handle)
573 
574  ! put the gamma density on grid
575  CALL timeset(routinen//"_Gpot", handle2)
576  CALL pw_zero(rho_r)
577  CALL pw_zero(rho_g)
578  CALL calculate_rho_elec(matrix_p=matrix_p_munu%matrix, &
579  rho=rho_r, &
580  rho_gspace=rho_g, &
581  task_list_external=task_list_sub, &
582  pw_env_external=pw_env_sub, &
583  ks_env=ks_env)
584  ! calculate associated hartree potential
585  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
586  CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
587  CALL timestop(handle2)
588 
589  IF (use_virial) CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
590 
591  ! integrate potential with auxiliary basis function derivatives
592  NULLIFY (rs_v)
593  NULLIFY (rs_descs)
594  CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
595  CALL potential_pw2rs(rs_v, rho_r, pw_env_sub)
596 
597  offset = 0
598  DO iatom = 1, SIZE(kind_of)
599  ikind = kind_of(iatom)
600  atom_a = atom_of_kind(iatom)
601  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
602  basis_type="RI_AUX")
603 
604  first_sgfa => basis_set_a%first_sgf
605  la_max => basis_set_a%lmax
606  la_min => basis_set_a%lmin
607  npgfa => basis_set_a%npgf
608  nseta = basis_set_a%nset
609  nsgfa => basis_set_a%nsgf_set
610  rpgfa => basis_set_a%pgf_radius
611  set_radius_a => basis_set_a%set_radius
612  sphi_a => basis_set_a%sphi
613  zeta => basis_set_a%zet
614 
615  ra(:) = pbc(particle_set(iatom)%r, cell)
616 
617  force_a(:) = 0.0_dp
618  force_b(:) = 0.0_dp
619  IF (use_virial) THEN
620  my_virial_a = 0.0_dp
621  my_virial_b = 0.0_dp
622  END IF
623 
624  DO iset = 1, nseta
625  ncoa = npgfa(iset)*ncoset(la_max(iset))
626  sgfa = first_sgfa(1, iset)
627 
628  ALLOCATE (i_tmp2(ncoa, 1))
629  i_tmp2 = 0.0_dp
630  ALLOCATE (i_ab(nsgfa(iset), 1))
631  i_ab = 0.0_dp
632  ALLOCATE (pab(ncoa, 1))
633  pab = 0.0_dp
634 
635  skip_shell = .true.
636  DO iorb = 1, nsgfa(iset)
637  IF (iorb + offset == lll) THEN
638  i_ab(iorb, 1) = 1.0_dp
639  skip_shell = .false.
640  END IF
641  END DO
642 
643  IF (skip_shell) THEN
644  offset = offset + nsgfa(iset)
645  DEALLOCATE (i_tmp2)
646  DEALLOCATE (i_ab)
647  DEALLOCATE (pab)
648  cycle
649  END IF
650 
651  CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
652  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
653  i_ab(1, 1), SIZE(i_ab, 1), &
654  0.0_dp, pab(1, 1), SIZE(pab, 1))
655  DEALLOCATE (i_ab)
656 
657  igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, minval(zeta(:, iset)))
658  map_it_here = .false.
659  use_subpatch = .NOT. all(rs_v(igrid_level)%desc%perd == 1)
660 
661  IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
662  DO ipgf = 1, npgfa(iset)
663  na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
664  na2 = ipgf*ncoset(la_max(iset))
665  igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))
666 
667  radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
668  lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
669  zetp=zeta(ipgf, iset), &
670  eps=dft_control%qs_control%eps_gvg_rspace, &
671  prefactor=1.0_dp, cutoff=1.0_dp)
672 
673  CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), &
674  lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, &
675  ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
676  rsgrid=rs_v(igrid_level), &
677  hab=i_tmp2, &
678  pab=pab, &
679  o1=na1 - 1, &
680  o2=0, &
681  radius=radius, &
682  calculate_forces=.true., &
683  force_a=force_a, force_b=force_b, &
684  use_virial=use_virial, &
685  my_virial_a=my_virial_a, &
686  my_virial_b=my_virial_b, &
687  use_subpatch=use_subpatch, &
688  subpatch_pattern=0)
689  END DO
690  END IF
691 
692  DEALLOCATE (i_tmp2)
693  DEALLOCATE (pab)
694 
695  offset = offset + nsgfa(iset)
696 
697  END DO
698 
699  force(ikind)%rho_elec(:, atom_a) = &
700  force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b(:)
701  IF (use_virial) THEN
702  h_stress = h_stress + my_virial_a + my_virial_b
703  END IF
704  END DO
705 
706  CALL timestop(handle)
707 
708  END SUBROUTINE integrate_potential_forces_3c_2c
709 
710 ! **************************************************************************************************
711 !> \brief Calculates stress tensor contribution from the operator
712 !> \param rho_g_copy ...
713 !> \param pot_g ...
714 !> \param rho_g ...
715 !> \param dvg ...
716 !> \param h_stress ...
717 !> \param potential_parameter ...
718 !> \param para_env_sub ...
719 ! **************************************************************************************************
720  SUBROUTINE virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
721 
722  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_g_copy, pot_g
723  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
724  TYPE(pw_c1d_gs_type), INTENT(IN) :: dvg(3)
725  REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
726  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
727  TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
728 
729  CHARACTER(LEN=*), PARAMETER :: routinen = 'virial_gpw_potential'
730 
731  INTEGER :: alpha, beta, handle
732  INTEGER, DIMENSION(3) :: comp
733  REAL(kind=dp) :: e_hartree
734 
735  ! add the volume contribution
736  CALL timeset(routinen, handle)
737  e_hartree = 0.0_dp
738  e_hartree = pw_integral_ab(rho_g_copy, pot_g)
739 
740  DO alpha = 1, 3
741  comp = 0
742  comp(alpha) = 1
743  CALL pw_copy(pot_g, rho_g)
744  CALL pw_derive(rho_g, comp)
745  CALL factor_virial_gpw(rho_g, potential_parameter)
746  h_stress(alpha, alpha) = h_stress(alpha, alpha) - e_hartree/real(para_env_sub%num_pe, dp)
747  DO beta = alpha, 3
748  h_stress(alpha, beta) = h_stress(alpha, beta) &
749  - 2.0_dp*pw_integral_ab(rho_g, dvg(beta))/fourpi/real(para_env_sub%num_pe, dp)
750  h_stress(beta, alpha) = h_stress(alpha, beta)
751  END DO
752  END DO
753 
754  CALL timestop(handle)
755  END SUBROUTINE virial_gpw_potential
756 
757 ! **************************************************************************************************
758 !> \brief Multiplies a potential in g space with the factor ln(V(g)/Vc(g))' with Vc(g) being the
759 !> Fourier-transformed of the Coulomg potential
760 !> \param pw ...
761 !> \param potential_parameter parameters of potential V(g)
762 ! **************************************************************************************************
763  SUBROUTINE factor_virial_gpw(pw, potential_parameter)
764  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
765  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
766 
767  SELECT CASE (potential_parameter%potential_type)
768  CASE (do_potential_coulomb)
769  RETURN
770  CASE (do_potential_long)
771  CALL pw_log_deriv_gauss(pw, potential_parameter%omega)
772  CASE (do_potential_short)
773  CALL pw_log_deriv_compl_gauss(pw, potential_parameter%omega)
774  CASE (do_potential_mix_cl)
775  CALL pw_log_deriv_mix_cl(pw, potential_parameter%omega, &
776  potential_parameter%scale_coulomb, potential_parameter%scale_longrange)
778  CALL pw_log_deriv_trunc(pw, potential_parameter%cutoff_radius)
779  CASE (do_potential_id)
780  CALL pw_zero(pw)
781  CASE DEFAULT
782  cpabort("Unknown potential type")
783  END SELECT
784 
785  END SUBROUTINE factor_virial_gpw
786 
787 ! **************************************************************************************************
788 !> \brief Integrate potential of RI function with RI function
789 !> \param pw_env_sub ...
790 !> \param pot_r ...
791 !> \param kind_of ...
792 !> \param atom_of_kind ...
793 !> \param particle_set ...
794 !> \param qs_kind_set ...
795 !> \param G_PQ_local ...
796 !> \param cell ...
797 !> \param force ...
798 !> \param use_virial ...
799 !> \param h_stress ...
800 !> \param para_env_sub ...
801 !> \param dft_control ...
802 ! **************************************************************************************************
803  SUBROUTINE integrate_potential(pw_env_sub, pot_r, kind_of, atom_of_kind, particle_set, qs_kind_set, &
804  G_PQ_local, cell, force, use_virial, h_stress, para_env_sub, dft_control)
805 
806  TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
807  TYPE(pw_r3d_rs_type), INTENT(IN) :: pot_r
808  INTEGER, DIMENSION(:), INTENT(IN) :: kind_of, atom_of_kind
809  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
810  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
811  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: g_pq_local
812  TYPE(cell_type), INTENT(IN), POINTER :: cell
813  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
814  LOGICAL, INTENT(IN) :: use_virial
815  REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
816  TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
817  TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
818 
819  CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_potential'
820 
821  INTEGER :: atom_a, handle, iatom, igrid_level, &
822  ikind, ipgf, iset, na1, na2, ncoa, &
823  nseta, offset, sgfa
824  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
825  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
826  LOGICAL :: use_subpatch
827  REAL(kind=dp) :: radius
828  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: i_ab
829  REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
830  REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
831  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a
832  REAL(kind=dp), DIMENSION(:, :), POINTER :: i_tmp2, pab, rpgfa, sphi_a, zeta
833  TYPE(gto_basis_set_type), POINTER :: basis_set_a
834  TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
835  POINTER :: rs_descs
836  TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
837 
838  CALL timeset(routinen, handle)
839 
840  NULLIFY (rs_v)
841  NULLIFY (rs_descs)
842  CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
843  CALL potential_pw2rs(rs_v, pot_r, pw_env_sub)
844 
845  offset = 0
846  DO iatom = 1, SIZE(kind_of)
847  ikind = kind_of(iatom)
848  atom_a = atom_of_kind(iatom)
849  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
850  basis_type="RI_AUX")
851 
852  first_sgfa => basis_set_a%first_sgf
853  la_max => basis_set_a%lmax
854  la_min => basis_set_a%lmin
855  npgfa => basis_set_a%npgf
856  nseta = basis_set_a%nset
857  nsgfa => basis_set_a%nsgf_set
858  rpgfa => basis_set_a%pgf_radius
859  set_radius_a => basis_set_a%set_radius
860  sphi_a => basis_set_a%sphi
861  zeta => basis_set_a%zet
862 
863  ra(:) = pbc(particle_set(iatom)%r, cell)
864 
865  force_a(:) = 0.0_dp
866  force_b(:) = 0.0_dp
867  IF (use_virial) THEN
868  my_virial_a = 0.0_dp
869  my_virial_b = 0.0_dp
870  END IF
871 
872  DO iset = 1, nseta
873  ncoa = npgfa(iset)*ncoset(la_max(iset))
874  sgfa = first_sgfa(1, iset)
875 
876  ALLOCATE (i_tmp2(ncoa, 1))
877  i_tmp2 = 0.0_dp
878  ALLOCATE (i_ab(nsgfa(iset), 1))
879  i_ab = 0.0_dp
880  ALLOCATE (pab(ncoa, 1))
881  pab = 0.0_dp
882 
883  i_ab(1:nsgfa(iset), 1) = -4.0_dp*g_pq_local(offset + 1:offset + nsgfa(iset))
884 
885  CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
886  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
887  i_ab(1, 1), SIZE(i_ab, 1), &
888  0.0_dp, pab(1, 1), SIZE(pab, 1))
889 
890  i_ab = 0.0_dp
891 
892  igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, minval(zeta(:, iset)))
893  use_subpatch = .NOT. all(rs_v(igrid_level)%desc%perd == 1)
894 
895  IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
896  DO ipgf = 1, npgfa(iset)
897  na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
898  na2 = ipgf*ncoset(la_max(iset))
899  igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))
900 
901  radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
902  lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
903  zetp=zeta(ipgf, iset), &
904  eps=dft_control%qs_control%eps_gvg_rspace, &
905  prefactor=1.0_dp, cutoff=1.0_dp)
906 
907  CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), &
908  lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, &
909  ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
910  rsgrid=rs_v(igrid_level), &
911  hab=i_tmp2, &
912  pab=pab, &
913  o1=na1 - 1, &
914  o2=0, &
915  radius=radius, &
916  calculate_forces=.true., &
917  force_a=force_a, force_b=force_b, &
918  use_virial=use_virial, &
919  my_virial_a=my_virial_a, &
920  my_virial_b=my_virial_b, &
921  use_subpatch=use_subpatch, &
922  subpatch_pattern=0)
923 
924  END DO
925 
926  END IF
927 
928  DEALLOCATE (i_tmp2)
929  DEALLOCATE (i_ab)
930  DEALLOCATE (pab)
931 
932  offset = offset + nsgfa(iset)
933 
934  END DO
935 
936  force(ikind)%rho_elec(:, atom_a) = &
937  force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b
938  IF (use_virial) THEN
939  h_stress = h_stress + my_virial_a + my_virial_b
940  END IF
941  END DO
942 
943  CALL timestop(handle)
944 
945  END SUBROUTINE
946 
947 ! **************************************************************************************************
948 !> \brief Prepares GPW calculation for RI-MP2/RI-RPA
949 !> \param qs_env ...
950 !> \param dft_control ...
951 !> \param e_cutoff_old ...
952 !> \param cutoff_old ...
953 !> \param relative_cutoff_old ...
954 !> \param para_env_sub ...
955 !> \param pw_env_sub ...
956 !> \param auxbas_pw_pool ...
957 !> \param poisson_env ...
958 !> \param task_list_sub ...
959 !> \param rho_r ...
960 !> \param rho_g ...
961 !> \param pot_g ...
962 !> \param psi_L ...
963 !> \param sab_orb_sub ...
964 ! **************************************************************************************************
965  SUBROUTINE prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
966  auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
967  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
968  TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
969  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
970  INTENT(OUT) :: e_cutoff_old
971  REAL(kind=dp), INTENT(OUT) :: cutoff_old, relative_cutoff_old
972  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
973  TYPE(pw_env_type), POINTER :: pw_env_sub
974  TYPE(pw_pool_type), INTENT(IN), POINTER :: auxbas_pw_pool
975  TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
976  TYPE(task_list_type), POINTER :: task_list_sub
977  TYPE(pw_r3d_rs_type), INTENT(OUT) :: rho_r
978  TYPE(pw_c1d_gs_type), INTENT(OUT) :: rho_g, pot_g
979  TYPE(pw_r3d_rs_type), INTENT(OUT) :: psi_l
980  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
981  INTENT(IN), POINTER :: sab_orb_sub
982 
983  CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_gpw'
984 
985  INTEGER :: handle, i_multigrid, n_multigrid
986  LOGICAL :: skip_load_balance_distributed
987  REAL(kind=dp) :: progression_factor
988  TYPE(qs_ks_env_type), POINTER :: ks_env
989 
990  CALL timeset(routinen, handle)
991 
992  CALL get_qs_env(qs_env, dft_control=dft_control, ks_env=ks_env)
993 
994  ! hack hack hack XXXXXXXXXXXXXXX rebuilds the pw_en with the new cutoffs
995  progression_factor = dft_control%qs_control%progression_factor
996  n_multigrid = SIZE(dft_control%qs_control%e_cutoff)
997  ALLOCATE (e_cutoff_old(n_multigrid))
998  e_cutoff_old(:) = dft_control%qs_control%e_cutoff
999  cutoff_old = dft_control%qs_control%cutoff
1000 
1001  dft_control%qs_control%cutoff = qs_env%mp2_env%mp2_gpw%cutoff*0.5_dp
1002  dft_control%qs_control%e_cutoff(1) = dft_control%qs_control%cutoff
1003  DO i_multigrid = 2, n_multigrid
1004  dft_control%qs_control%e_cutoff(i_multigrid) = dft_control%qs_control%e_cutoff(i_multigrid - 1) &
1005  /progression_factor
1006  END DO
1007 
1008  relative_cutoff_old = dft_control%qs_control%relative_cutoff
1009  dft_control%qs_control%relative_cutoff = qs_env%mp2_env%mp2_gpw%relative_cutoff*0.5_dp
1010 
1011  ! a pw_env
1012  NULLIFY (pw_env_sub)
1013  CALL pw_env_create(pw_env_sub)
1014  CALL pw_env_rebuild(pw_env_sub, qs_env, para_env_sub)
1015 
1016  CALL pw_env_get(pw_env_sub, auxbas_pw_pool=auxbas_pw_pool, &
1017  poisson_env=poisson_env)
1018  ! hack hack hack XXXXXXXXXXXXXXX
1019 
1020  ! now we need a task list, hard code skip_load_balance_distributed
1021  NULLIFY (task_list_sub)
1022  skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
1023  CALL allocate_task_list(task_list_sub)
1024  CALL generate_qs_task_list(ks_env, task_list_sub, &
1025  reorder_rs_grid_ranks=.true., soft_valid=.false., &
1026  skip_load_balance_distributed=skip_load_balance_distributed, &
1027  pw_env_external=pw_env_sub, sab_orb_external=sab_orb_sub)
1028 
1029  ! get some of the grids ready
1030  CALL auxbas_pw_pool%create_pw(rho_r)
1031  CALL auxbas_pw_pool%create_pw(rho_g)
1032  CALL auxbas_pw_pool%create_pw(pot_g)
1033  CALL auxbas_pw_pool%create_pw(psi_l)
1034 
1035  ! run the FFT once, to set up buffers and to take into account the memory
1036  CALL pw_zero(rho_r)
1037  CALL pw_transfer(rho_r, rho_g)
1038 
1039  CALL timestop(handle)
1040 
1041  END SUBROUTINE prepare_gpw
1042 
1043 ! **************************************************************************************************
1044 !> \brief Cleanup GPW integration for RI-MP2/RI-RPA
1045 !> \param qs_env ...
1046 !> \param e_cutoff_old ...
1047 !> \param cutoff_old ...
1048 !> \param relative_cutoff_old ...
1049 !> \param para_env_sub ...
1050 !> \param pw_env_sub ...
1051 !> \param task_list_sub ...
1052 !> \param auxbas_pw_pool ...
1053 !> \param rho_r ...
1054 !> \param rho_g ...
1055 !> \param pot_g ...
1056 !> \param psi_L ...
1057 ! **************************************************************************************************
1058  SUBROUTINE cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
1059  task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
1060  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1061  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1062  INTENT(IN) :: e_cutoff_old
1063  REAL(kind=dp), INTENT(IN) :: cutoff_old, relative_cutoff_old
1064  TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
1065  TYPE(pw_env_type), POINTER :: pw_env_sub
1066  TYPE(task_list_type), POINTER :: task_list_sub
1067  TYPE(pw_pool_type), INTENT(IN), POINTER :: auxbas_pw_pool
1068  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
1069  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g, pot_g
1070  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: psi_l
1071 
1072  CHARACTER(LEN=*), PARAMETER :: routinen = 'cleanup_gpw'
1073 
1074  INTEGER :: handle
1075  TYPE(dft_control_type), POINTER :: dft_control
1076 
1077  CALL timeset(routinen, handle)
1078 
1079  ! and now free the whole lot
1080  CALL auxbas_pw_pool%give_back_pw(rho_r)
1081  CALL auxbas_pw_pool%give_back_pw(rho_g)
1082  CALL auxbas_pw_pool%give_back_pw(pot_g)
1083  CALL auxbas_pw_pool%give_back_pw(psi_l)
1084 
1085  CALL deallocate_task_list(task_list_sub)
1086  CALL pw_env_release(pw_env_sub, para_env=para_env_sub)
1087 
1088  CALL get_qs_env(qs_env, dft_control=dft_control)
1089 
1090  ! restore the initial value of the cutoff
1091  dft_control%qs_control%e_cutoff = e_cutoff_old
1092  dft_control%qs_control%cutoff = cutoff_old
1093  dft_control%qs_control%relative_cutoff = relative_cutoff_old
1094 
1095  CALL timestop(handle)
1096 
1097  END SUBROUTINE cleanup_gpw
1098 
1099 ! **************************************************************************************************
1100 !> \brief Calculates potential from a given density in g-space
1101 !> \param pot_r on output: potential in r space
1102 !> \param rho_g on input: rho in g space
1103 !> \param poisson_env ...
1104 !> \param pot_g on output: potential in g space
1105 !> \param potential_parameter Potential parameters, if not provided, assume Coulomb potential
1106 !> \param dvg in output: first derivatives of the corresponding Coulomb potential
1107 !> \param no_transfer whether NOT to transform potential from g-space to r-space (default: do it)
1108 ! **************************************************************************************************
1109  SUBROUTINE calc_potential_gpw(pot_r, rho_g, poisson_env, pot_g, potential_parameter, dvg, no_transfer)
1110  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pot_r
1111  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_g
1112  TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
1113  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
1114  TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
1115  TYPE(pw_c1d_gs_type), DIMENSION(3), &
1116  INTENT(INOUT), OPTIONAL :: dvg
1117  LOGICAL, INTENT(IN), OPTIONAL :: no_transfer
1118 
1119  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_potential_gpw'
1120 
1121  INTEGER :: comp(3), handle, idir, my_potential_type
1122  LOGICAL :: my_transfer
1123 
1124  CALL timeset(routinen, handle)
1125 
1126  my_potential_type = do_potential_coulomb
1127  IF (PRESENT(potential_parameter)) THEN
1128  my_potential_type = potential_parameter%potential_type
1129  END IF
1130 
1131  my_transfer = .true.
1132  IF (PRESENT(no_transfer)) my_transfer = .NOT. no_transfer
1133 
1134  ! in case we do Coulomb metric for RI, we need the Coulomb operator, but for RI with the
1135  ! overlap metric, we do not need the Coulomb operator
1136  IF (my_potential_type /= do_potential_id) THEN
1137  IF (PRESENT(dvg)) THEN
1138  CALL pw_poisson_solve(poisson_env, rho_g, vhartree=pot_g, dvhartree=dvg)
1139  ELSE
1140  CALL pw_poisson_solve(poisson_env, rho_g, vhartree=pot_g)
1141  END IF
1142  IF (my_potential_type == do_potential_long) CALL pw_gauss_damp(pot_g, potential_parameter%omega)
1143  IF (my_potential_type == do_potential_short) CALL pw_compl_gauss_damp(pot_g, potential_parameter%omega)
1144  IF (my_potential_type == do_potential_truncated) CALL pw_truncated(pot_g, potential_parameter%cutoff_radius)
1145  IF (my_potential_type == do_potential_mix_cl) CALL pw_gauss_damp_mix(pot_g, potential_parameter%omega, &
1146  potential_parameter%scale_coulomb, &
1147  potential_parameter%scale_longrange)
1148  IF (my_transfer) CALL pw_transfer(pot_g, pot_r)
1149  ELSE
1150  ! If we use an overlap metric, make sure to use the correct potential=density on output
1151  CALL pw_copy(rho_g, pot_g)
1152  IF (PRESENT(dvg)) THEN
1153  DO idir = 1, 3
1154  CALL pw_copy(pot_g, dvg(idir))
1155  comp = 0
1156  comp(idir) = 1
1157  CALL pw_derive(dvg(idir), comp)
1158  END DO
1159  END IF
1160  IF (my_transfer) CALL pw_transfer(rho_g, pot_r)
1161  END IF
1162 
1163  IF (my_transfer) CALL pw_scale(pot_r, pot_r%pw_grid%dvol)
1164  CALL timestop(handle)
1165 
1166  END SUBROUTINE calc_potential_gpw
1167 
1168 END MODULE mp2_eri_gpw
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
All kind of helpful little routines.
Definition: ao_util.F:14
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition: ao_util.F:209
Define the atomic kind types and their sub types.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_potential_mix_cl
integer, parameter, public do_potential_truncated
integer, parameter, public do_potential_id
integer, parameter, public do_potential_coulomb
integer, parameter, public do_potential_short
integer, parameter, public do_potential_long
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition: libint_2c_3c.F:14
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public fourpi
Interface to the message passing library MPI.
Routines to calculate 2c- and 3c-integrals for RI with GPW.
Definition: mp2_eri_gpw.F:13
subroutine, public cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
Cleanup GPW integration for RI-MP2/RI-RPA.
Definition: mp2_eri_gpw.F:1060
subroutine, public integrate_potential_forces_3c_1c(mat_munu, rho_r, matrix_P_munu, qs_env, pw_env_sub, task_list_sub)
Takes the precomputed potential of an RI wave-function and determines matrix element and gradients wi...
Definition: mp2_eri_gpw.F:472
subroutine, public integrate_potential_forces_2c(rho_r, LLL, matrix, rho_g, atomic_kind_set, qs_kind_set, particle_set, cell, pw_env_sub, poisson_env, pot_g, potential_parameter, use_virial, rho_g_copy, dvg, kind_of, atom_of_kind, G_PQ_local, force, h_stress, para_env_sub, dft_control, psi_L, factor)
Integrates the potential of a RI function obtaining the forces and stress tensor.
Definition: mp2_eri_gpw.F:383
subroutine, public virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
Calculates stress tensor contribution from the operator.
Definition: mp2_eri_gpw.F:721
subroutine, public integrate_potential_forces_3c_2c(matrix_P_munu, rho_r, rho_g, task_list_sub, pw_env_sub, potential_parameter, ks_env, poisson_env, pot_g, use_virial, rho_g_copy, dvg, h_stress, para_env_sub, kind_of, atom_of_kind, qs_kind_set, particle_set, cell, LLL, force, dft_control)
Integrates potential of two Gaussians to a potential.
Definition: mp2_eri_gpw.F:527
subroutine, public prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
Prepares GPW calculation for RI-MP2/RI-RPA.
Definition: mp2_eri_gpw.F:967
subroutine, public calc_potential_gpw(pot_r, rho_g, poisson_env, pot_g, potential_parameter, dvg, no_transfer)
Calculates potential from a given density in g-space.
Definition: mp2_eri_gpw.F:1110
subroutine, public mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, natom, potential_parameter, sab_orb_sub, L_local_col, kind_of)
Integrates the potential of an RI function.
Definition: mp2_eri_gpw.F:170
subroutine, public mp2_eri_3c_integrate_gpw(mo_coeff, psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env_sub, external_vector, poisson_env, rho_r, pot_g, potential_parameter, mat_munu, qs_env, task_list_sub)
...
Definition: mp2_eri_gpw.F:110
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
Define the data structure for the particle information.
methods of pw_env that have dependence on qs_env
subroutine, public pw_env_rebuild(pw_env, qs_env, external_para_env)
rebuilds the pw_env data (necessary if cell or cutoffs change)
subroutine, public pw_env_create(pw_env)
creates a pw_env, if qs_env is given calls pw_env_rebuild
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_release(pw_env, para_env)
releases the given pw_env (see doc/ReferenceCounting.html)
Definition: pw_env_types.F:176
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
subroutine, public pw_log_deriv_trunc(pw, rcutoff)
Multiply all data points with the logarithmic derivative of the complementary cosine This function is...
Definition: pw_methods.F:10065
subroutine, public pw_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
Definition: pw_methods.F:9806
subroutine, public pw_log_deriv_compl_gauss(pw, omega)
Multiply all data points with the logarithmic derivative of the complementary Gaussian damping factor...
Definition: pw_methods.F:9911
subroutine, public pw_log_deriv_gauss(pw, omega)
Multiply all data points with the logarithmic derivative of a Gaussian.
Definition: pw_methods.F:9837
subroutine, public pw_gauss_damp_mix(pw, omega, scale_coul, scale_long)
Multiply all data points with a Gaussian damping factor and mixes it with the original function Neede...
Definition: pw_methods.F:9958
subroutine, public pw_truncated(pw, rcutoff)
Multiply all data points with a complementary cosine Needed for truncated Coulomb potential V(\vec r)...
Definition: pw_methods.F:10030
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition: pw_methods.F:10106
subroutine, public pw_compl_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
Definition: pw_methods.F:9873
subroutine, public pw_log_deriv_mix_cl(pw, omega, scale_coul, scale_long)
Multiply all data points with the logarithmic derivative of the mixed longrange/Coulomb potential Nee...
Definition: pw_methods.F:9990
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
subroutine, public collocate_single_gaussian(rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, required_function, basis_type)
maps a single gaussian on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Define the neighbor list data types and the corresponding functionality.
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
Transfers densities from PW to RS grids and potentials from PW to RS.
subroutine, public potential_pw2rs(rs_v, v_rspace, pw_env)
transfers a potential from a pw_grid to a vector of realspace multigrids
generate the tasks lists used by collocate and integrate routines
subroutine, public generate_qs_task_list(ks_env, task_list, reorder_rs_grid_ranks, skip_load_balance_distributed, soft_valid, basis_type, pw_env_external, sab_orb_external)
...
types for task lists
subroutine, public deallocate_task_list(task_list)
deallocates the components and the object itself
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type