Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : ! **************************************************************************************************
13 : MODULE mp2_eri_gpw
14 : USE ao_util, ONLY: exp_radius_very_extended
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_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_id,&
25 : do_potential_long,&
26 : do_potential_mix_cl,&
27 : do_potential_short,&
28 : do_potential_truncated
29 : USE kinds, ONLY: dp
30 : USE libint_2c_3c, ONLY: libint_potential_type
31 : USE mathconstants, ONLY: fourpi
32 : USE message_passing, ONLY: mp_para_env_type
33 : USE orbital_pointers, ONLY: ncoset
34 : USE particle_types, ONLY: particle_type
35 : USE pw_env_methods, ONLY: pw_env_create,&
36 : pw_env_rebuild
37 : USE pw_env_types, ONLY: pw_env_get,&
38 : pw_env_release,&
39 : pw_env_type
40 : USE pw_methods, ONLY: &
41 : pw_compl_gauss_damp, pw_copy, pw_derive, pw_gauss_damp, pw_gauss_damp_mix, pw_integral_ab, &
42 : pw_log_deriv_compl_gauss, pw_log_deriv_gauss, pw_log_deriv_mix_cl, pw_log_deriv_trunc, &
43 : pw_scale, pw_transfer, pw_truncated, pw_zero
44 : USE pw_poisson_methods, ONLY: pw_poisson_solve
45 : USE pw_poisson_types, ONLY: pw_poisson_type
46 : USE pw_pool_types, ONLY: pw_pool_type
47 : USE pw_types, ONLY: pw_c1d_gs_type,&
48 : pw_r3d_rs_type
49 : USE qs_collocate_density, ONLY: calculate_rho_elec,&
50 : collocate_function,&
51 : collocate_single_gaussian
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_force_types, ONLY: qs_force_type
55 : USE qs_integrate_potential, ONLY: integrate_pgf_product,&
56 : integrate_v_rspace
57 : USE qs_kind_types, ONLY: get_qs_kind,&
58 : get_qs_kind_set,&
59 : qs_kind_type
60 : USE qs_ks_types, ONLY: qs_ks_env_type
61 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
62 : USE realspace_grid_types, ONLY: map_gaussian_here,&
63 : realspace_grid_desc_p_type,&
64 : realspace_grid_type
65 : USE rs_pw_interface, ONLY: potential_pw2rs
66 : USE task_list_methods, ONLY: generate_qs_task_list
67 : USE task_list_types, ONLY: allocate_task_list,&
68 : deallocate_task_list,&
69 : task_list_type
70 :
71 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 :
76 : PRIVATE
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_eri_gpw'
79 :
80 : PUBLIC :: mp2_eri_2c_integrate_gpw, mp2_eri_3c_integrate_gpw, calc_potential_gpw, cleanup_gpw, prepare_gpw, &
81 : integrate_potential_forces_2c, integrate_potential_forces_3c_1c, integrate_potential_forces_3c_2c, &
82 : virial_gpw_potential
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief ...
88 : !> \param psi_L ...
89 : !> \param rho_g ...
90 : !> \param atomic_kind_set ...
91 : !> \param qs_kind_set ...
92 : !> \param cell ...
93 : !> \param dft_control ...
94 : !> \param particle_set ...
95 : !> \param pw_env_sub ...
96 : !> \param external_vector ...
97 : !> \param poisson_env ...
98 : !> \param rho_r ...
99 : !> \param pot_g ...
100 : !> \param potential_parameter ...
101 : !> \param mat_munu ...
102 : !> \param qs_env ...
103 : !> \param task_list_sub ...
104 : ! **************************************************************************************************
105 29518 : SUBROUTINE mp2_eri_3c_integrate_gpw(psi_L, rho_g, atomic_kind_set, qs_kind_set, &
106 : cell, dft_control, particle_set, &
107 14759 : pw_env_sub, external_vector, poisson_env, rho_r, pot_g, &
108 : potential_parameter, mat_munu, qs_env, task_list_sub)
109 :
110 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: psi_L
111 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
112 : TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
113 : POINTER :: atomic_kind_set
114 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
115 : POINTER :: qs_kind_set
116 : TYPE(cell_type), INTENT(IN), POINTER :: cell
117 : TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
118 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
119 : POINTER :: particle_set
120 : TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
121 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: external_vector
122 : TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
123 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
124 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
125 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
126 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu
127 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
128 : TYPE(task_list_type), INTENT(IN), POINTER :: task_list_sub
129 :
130 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_eri_3c_integrate_gpw'
131 :
132 : INTEGER :: handle
133 :
134 14759 : CALL timeset(routineN, handle)
135 :
136 : ! pseudo psi_L
137 : CALL collocate_function(external_vector, psi_L, rho_g, atomic_kind_set, &
138 : qs_kind_set, cell, particle_set, pw_env_sub, &
139 14759 : dft_control%qs_control%eps_rho_rspace, basis_type="RI_AUX")
140 :
141 14759 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
142 :
143 : ! and finally (K|mu nu)
144 14759 : CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
145 : CALL integrate_v_rspace(rho_r, hmat=mat_munu, qs_env=qs_env, &
146 : calculate_forces=.FALSE., compute_tau=.FALSE., gapw=.FALSE., &
147 14759 : pw_env_external=pw_env_sub, task_list_external=task_list_sub)
148 :
149 14759 : CALL timestop(handle)
150 :
151 14759 : END SUBROUTINE mp2_eri_3c_integrate_gpw
152 :
153 : ! **************************************************************************************************
154 : !> \brief Integrates the potential of an RI function
155 : !> \param qs_env ...
156 : !> \param para_env_sub ...
157 : !> \param my_group_L_start ...
158 : !> \param my_group_L_end ...
159 : !> \param natom ...
160 : !> \param potential_parameter ...
161 : !> \param sab_orb_sub ...
162 : !> \param L_local_col ...
163 : !> \param kind_of ...
164 : ! **************************************************************************************************
165 398 : SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, &
166 398 : natom, potential_parameter, sab_orb_sub, L_local_col, kind_of)
167 :
168 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
169 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
170 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, natom
171 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
172 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
173 : INTENT(IN), POINTER :: sab_orb_sub
174 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: L_local_col
175 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
176 :
177 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_eri_2c_integrate_gpw'
178 :
179 : INTEGER :: dir, handle, handle2, i_counter, iatom, igrid_level, ikind, ipgf, iset, lb(3), &
180 : LLL, location(3), max_nseta, na1, na2, ncoa, nseta, offset, sgfa, tp(3), ub(3)
181 398 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: offset_2d
182 398 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
183 398 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
184 : LOGICAL :: map_it_here, use_subpatch
185 : REAL(KIND=dp) :: cutoff_old, radius, relative_cutoff_old
186 398 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old
187 398 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab
188 : REAL(KIND=dp), DIMENSION(3) :: ra
189 398 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
190 398 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, rpgfa, sphi_a, zeta
191 398 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
192 : TYPE(cell_type), POINTER :: cell
193 : TYPE(dft_control_type), POINTER :: dft_control
194 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
195 398 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
196 : TYPE(pw_c1d_gs_type) :: pot_g, rho_g
197 : TYPE(pw_env_type), POINTER :: pw_env_sub
198 : TYPE(pw_poisson_type), POINTER :: poisson_env
199 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
200 : TYPE(pw_r3d_rs_type) :: psi_L, rho_r
201 398 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
202 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
203 398 : POINTER :: rs_descs
204 398 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
205 : TYPE(task_list_type), POINTER :: task_list_sub
206 :
207 398 : 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 398 : auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
211 :
212 398 : CALL get_qs_env(qs_env, cell=cell, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, particle_set=particle_set)
213 :
214 1468442 : L_local_col = 0.0_dp
215 :
216 398 : i_counter = 0
217 17667 : DO LLL = my_group_L_start, my_group_L_end
218 17269 : i_counter = i_counter + 1
219 :
220 : ! pseudo psi_L
221 : CALL collocate_single_gaussian(psi_L, rho_g, atomic_kind_set, &
222 : qs_kind_set, cell, dft_control, particle_set, pw_env_sub, &
223 17269 : required_function=LLL, basis_type="RI_AUX")
224 :
225 17269 : CALL timeset(routineN//"_pot_lm", handle2)
226 :
227 17269 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
228 :
229 17269 : NULLIFY (rs_v)
230 17269 : NULLIFY (rs_descs)
231 17269 : CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
232 17269 : CALL potential_pw2rs(rs_v, rho_r, pw_env_sub)
233 :
234 17269 : CALL timestop(handle2)
235 :
236 17269 : offset = 0
237 : ! Prepare offset ahead of OMP parallel loop
238 17269 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxnset=max_nseta, basis_type="RI_AUX")
239 69076 : ALLOCATE (offset_2d(max_nseta, natom))
240 :
241 71012 : DO iatom = 1, natom
242 53743 : ikind = kind_of(iatom)
243 53743 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
244 53743 : nseta = basis_set_a%nset
245 53743 : nsgfa => basis_set_a%nsgf_set
246 576731 : DO iset = 1, nseta
247 505719 : offset = offset + nsgfa(iset)
248 559462 : offset_2d(iset, iatom) = offset
249 : END DO
250 : END DO
251 :
252 : ! integrate the little bastards
253 : !$OMP PARALLEL DO DEFAULT(NONE) &
254 : !$OMP SHARED(natom, particle_set, cell, pw_env_sub, rs_v, offset_2d, &
255 : !$OMP qs_kind_set, ncoset, para_env_sub, dft_control, i_counter, &
256 : !$OMP kind_of, l_local_col) &
257 : !$OMP PRIVATE(iatom, ikind, basis_set_a, first_sgfa, la_max, la_min, npgfa, &
258 : !$OMP nseta, nsgfa, rpgfa, set_radius_a, sphi_a, zeta, &
259 : !$OMP ra, iset, ncoa, I_tmp2, I_ab, igrid_level, &
260 : !$OMP map_it_here, dir, tp, lb, ub, location, ipgf, &
261 17269 : !$OMP sgfa, na1, na2, radius, offset, use_subpatch)
262 : DO iatom = 1, natom
263 : ikind = kind_of(iatom)
264 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
265 :
266 : first_sgfa => basis_set_a%first_sgf
267 : la_max => basis_set_a%lmax
268 : la_min => basis_set_a%lmin
269 : npgfa => basis_set_a%npgf
270 : nseta = basis_set_a%nset
271 : nsgfa => basis_set_a%nsgf_set
272 : rpgfa => basis_set_a%pgf_radius
273 : set_radius_a => basis_set_a%set_radius
274 : sphi_a => basis_set_a%sphi
275 : zeta => basis_set_a%zet
276 :
277 : ra(:) = pbc(particle_set(iatom)%r, cell)
278 :
279 : DO iset = 1, nseta
280 : ncoa = npgfa(iset)*ncoset(la_max(iset))
281 : sgfa = first_sgfa(1, iset)
282 :
283 : ALLOCATE (I_tmp2(ncoa, 1))
284 : I_tmp2 = 0.0_dp
285 : ALLOCATE (I_ab(nsgfa(iset), 1))
286 : I_ab = 0.0_dp
287 :
288 : offset = offset_2d(iset, iatom)
289 :
290 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, MINVAL(zeta(:, iset)))
291 : use_subpatch = .NOT. ALL(rs_v(igrid_level)%desc%perd == 1)
292 :
293 : IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, &
294 : offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
295 : DO ipgf = 1, npgfa(iset)
296 : sgfa = first_sgfa(1, iset)
297 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
298 : na2 = ipgf*ncoset(la_max(iset))
299 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))
300 :
301 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
302 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
303 : zetp=zeta(ipgf, iset), &
304 : eps=dft_control%qs_control%eps_gvg_rspace, &
305 : prefactor=1.0_dp, cutoff=1.0_dp)
306 :
307 : CALL integrate_pgf_product( &
308 : la_max=la_max(iset), zeta=zeta(ipgf, iset), la_min=la_min(iset), &
309 : lb_max=0, zetb=0.0_dp, lb_min=0, &
310 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
311 : rsgrid=rs_v(igrid_level), &
312 : hab=I_tmp2, &
313 : o1=na1 - 1, &
314 : o2=0, &
315 : radius=radius, &
316 : calculate_forces=.FALSE., &
317 : use_subpatch=use_subpatch, &
318 : subpatch_pattern=0)
319 : END DO
320 :
321 : CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, &
322 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
323 : I_tmp2(1, 1), SIZE(I_tmp2, 1), &
324 : 1.0_dp, I_ab(1, 1), SIZE(I_ab, 1))
325 :
326 : L_local_col(offset - nsgfa(iset) + 1:offset, i_counter) = I_ab(1:nsgfa(iset), 1)
327 : END IF
328 :
329 : DEALLOCATE (I_tmp2)
330 : DEALLOCATE (I_ab)
331 :
332 : END DO
333 : END DO
334 : !$OMP END PARALLEL DO
335 52205 : DEALLOCATE (offset_2d)
336 :
337 : END DO
338 :
339 2936486 : CALL para_env_sub%sum(L_local_col)
340 :
341 : CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
342 398 : task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
343 :
344 398 : CALL timestop(handle)
345 :
346 796 : END SUBROUTINE
347 :
348 : ! **************************************************************************************************
349 : !> \brief Integrates the potential of a RI function obtaining the forces and stress tensor
350 : !> \param rho_r ...
351 : !> \param LLL ...
352 : !> \param rho_g ...
353 : !> \param atomic_kind_set ...
354 : !> \param qs_kind_set ...
355 : !> \param particle_set ...
356 : !> \param cell ...
357 : !> \param pw_env_sub ...
358 : !> \param poisson_env ...
359 : !> \param pot_g ...
360 : !> \param potential_parameter ...
361 : !> \param use_virial ...
362 : !> \param rho_g_copy ...
363 : !> \param dvg ...
364 : !> \param kind_of ...
365 : !> \param atom_of_kind ...
366 : !> \param G_PQ_local ...
367 : !> \param force ...
368 : !> \param h_stress ...
369 : !> \param para_env_sub ...
370 : !> \param dft_control ...
371 : !> \param psi_L ...
372 : !> \param factor ...
373 : ! **************************************************************************************************
374 33912 : SUBROUTINE integrate_potential_forces_2c(rho_r, LLL, rho_g, atomic_kind_set, &
375 : qs_kind_set, particle_set, cell, pw_env_sub, poisson_env, pot_g, &
376 : potential_parameter, use_virial, rho_g_copy, dvg, &
377 11304 : kind_of, atom_of_kind, G_PQ_local, force, h_stress, para_env_sub, &
378 : dft_control, psi_L, factor)
379 :
380 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
381 : INTEGER, INTENT(IN) :: LLL
382 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
383 : TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
384 : POINTER :: atomic_kind_set
385 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
386 : POINTER :: qs_kind_set
387 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
388 : POINTER :: particle_set
389 : TYPE(cell_type), INTENT(IN), POINTER :: cell
390 : TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
391 : TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
392 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
393 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
394 : LOGICAL, INTENT(IN) :: use_virial
395 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g_copy, dvg(3)
396 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of, atom_of_kind
397 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: G_PQ_local
398 : TYPE(qs_force_type), DIMENSION(:), INTENT(IN), &
399 : POINTER :: force
400 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
401 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
402 : TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
403 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: psi_L
404 : REAL(KIND=dp), INTENT(IN) :: factor
405 :
406 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential_forces_2c'
407 :
408 : INTEGER :: handle, handle2
409 :
410 11304 : CALL timeset(routineN, handle)
411 :
412 : ! calculate potential associated to the single aux function
413 11304 : CALL timeset(routineN//"_wf_pot", handle2)
414 : ! pseudo psi_L
415 11304 : CALL pw_zero(rho_r)
416 11304 : CALL pw_zero(rho_g)
417 : CALL collocate_single_gaussian(rho_r, rho_g, atomic_kind_set, &
418 : qs_kind_set, cell, dft_control, particle_set, &
419 11304 : pw_env_sub, required_function=LLL, basis_type='RI_AUX')
420 11304 : IF (use_virial) THEN
421 2325 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter, dvg)
422 : ELSE
423 8979 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
424 : END IF
425 11304 : CALL timestop(handle2)
426 :
427 11304 : IF (use_virial) THEN
428 : ! make a copy of the density in G space
429 2325 : CALL pw_copy(rho_g, rho_g_copy)
430 :
431 : ! add the volume contribution to the virial due to
432 : ! the (P|Q) integrals, first we put the full gamme_PQ
433 : ! pseudo wave-function into grid in order to calculate the
434 : ! hartree potential derivatives
435 2325 : CALL pw_zero(psi_L)
436 2325 : CALL pw_zero(rho_g)
437 : CALL collocate_function(0.5_dp*factor*G_PQ_local, psi_L, rho_g, atomic_kind_set, &
438 : qs_kind_set, cell, particle_set, pw_env_sub, &
439 : dft_control%qs_control%eps_rho_rspace, &
440 204624 : basis_type="RI_AUX")
441 : ! transfer to reciprocal space and calculate potential
442 2325 : CALL calc_potential_gpw(psi_L, rho_g, poisson_env, pot_g, potential_parameter, no_transfer=.TRUE.)
443 : ! update virial with volume term (first calculate hartree like energy (diagonal part of the virial))
444 2325 : CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
445 : END IF
446 :
447 : ! integrate the potential of the single gaussian and update
448 : ! 2-center forces with Gamma_PQ
449 : CALL integrate_potential(pw_env_sub, rho_r, kind_of, atom_of_kind, particle_set, qs_kind_set, &
450 971338 : -0.25_dp*factor*G_PQ_local, cell, force, use_virial, h_stress, para_env_sub, dft_control)
451 :
452 11304 : CALL timestop(handle)
453 11304 : END SUBROUTINE integrate_potential_forces_2c
454 :
455 : ! **************************************************************************************************
456 : !> \brief Takes the precomputed potential of an RI wave-function and determines matrix element and
457 : !> gradients with product of Gaussians
458 : !> \param mat_munu ...
459 : !> \param rho_r ...
460 : !> \param matrix_P_munu ...
461 : !> \param qs_env ...
462 : !> \param pw_env_sub ...
463 : !> \param task_list_sub ...
464 : ! **************************************************************************************************
465 10806 : SUBROUTINE integrate_potential_forces_3c_1c(mat_munu, rho_r, matrix_P_munu, qs_env, pw_env_sub, &
466 : task_list_sub)
467 :
468 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu
469 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_r
470 : TYPE(dbcsr_p_type), INTENT(IN) :: matrix_P_munu
471 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
472 : TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
473 : TYPE(task_list_type), INTENT(INOUT), POINTER :: task_list_sub
474 :
475 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential_forces_3c_1c'
476 :
477 : INTEGER :: handle
478 :
479 10806 : CALL timeset(routineN, handle)
480 :
481 : ! integrate the potential of the single gaussian and update
482 : ! 3-center forces
483 10806 : CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
484 : CALL integrate_v_rspace(rho_r, hmat=mat_munu, pmat=matrix_P_munu, &
485 : qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.FALSE., gapw=.FALSE., &
486 : pw_env_external=pw_env_sub, &
487 10806 : task_list_external=task_list_sub)
488 :
489 10806 : CALL timestop(handle)
490 10806 : END SUBROUTINE integrate_potential_forces_3c_1c
491 :
492 : ! **************************************************************************************************
493 : !> \brief Integrates potential of two Gaussians to a potential
494 : !> \param matrix_P_munu ...
495 : !> \param rho_r ...
496 : !> \param rho_g ...
497 : !> \param task_list_sub ...
498 : !> \param pw_env_sub ...
499 : !> \param potential_parameter ...
500 : !> \param ks_env ...
501 : !> \param poisson_env ...
502 : !> \param pot_g ...
503 : !> \param use_virial ...
504 : !> \param rho_g_copy ...
505 : !> \param dvg ...
506 : !> \param h_stress ...
507 : !> \param para_env_sub ...
508 : !> \param kind_of ...
509 : !> \param atom_of_kind ...
510 : !> \param qs_kind_set ...
511 : !> \param particle_set ...
512 : !> \param cell ...
513 : !> \param LLL ...
514 : !> \param force ...
515 : !> \param dft_control ...
516 : ! **************************************************************************************************
517 10806 : SUBROUTINE integrate_potential_forces_3c_2c(matrix_P_munu, rho_r, rho_g, task_list_sub, pw_env_sub, &
518 : potential_parameter, &
519 : ks_env, poisson_env, pot_g, use_virial, rho_g_copy, dvg, &
520 10806 : h_stress, para_env_sub, kind_of, atom_of_kind, &
521 : qs_kind_set, particle_set, cell, LLL, force, dft_control)
522 :
523 : TYPE(dbcsr_p_type), INTENT(IN) :: matrix_P_munu
524 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
525 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
526 : TYPE(task_list_type), INTENT(IN), POINTER :: task_list_sub
527 : TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
528 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
529 : TYPE(qs_ks_env_type), INTENT(IN), POINTER :: ks_env
530 : TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
531 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
532 : LOGICAL, INTENT(IN) :: use_virial
533 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g_copy
534 : TYPE(pw_c1d_gs_type), INTENT(IN) :: dvg(3)
535 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
536 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
537 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of, atom_of_kind
538 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
539 : POINTER :: qs_kind_set
540 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
541 : POINTER :: particle_set
542 : TYPE(cell_type), INTENT(IN), POINTER :: cell
543 : INTEGER, INTENT(IN) :: LLL
544 : TYPE(qs_force_type), DIMENSION(:), INTENT(IN), &
545 : POINTER :: force
546 : TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
547 :
548 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential_forces_3c_2c'
549 :
550 : INTEGER :: atom_a, handle, handle2, iatom, &
551 : igrid_level, ikind, iorb, ipgf, iset, &
552 : na1, na2, ncoa, nseta, offset, sgfa
553 10806 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
554 10806 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
555 : LOGICAL :: map_it_here, skip_shell, use_subpatch
556 : REAL(KIND=dp) :: radius
557 10806 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab
558 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
559 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
560 10806 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
561 10806 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, pab, rpgfa, sphi_a, zeta
562 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
563 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
564 10806 : POINTER :: rs_descs
565 10806 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
566 :
567 10806 : CALL timeset(routineN, handle)
568 :
569 : ! put the gamma density on grid
570 10806 : CALL timeset(routineN//"_Gpot", handle2)
571 10806 : CALL pw_zero(rho_r)
572 10806 : CALL pw_zero(rho_g)
573 : CALL calculate_rho_elec(matrix_p=matrix_P_munu%matrix, &
574 : rho=rho_r, &
575 : rho_gspace=rho_g, &
576 : task_list_external=task_list_sub, &
577 : pw_env_external=pw_env_sub, &
578 10806 : ks_env=ks_env)
579 : ! calculate associated hartree potential
580 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581 10806 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
582 10806 : CALL timestop(handle2)
583 :
584 10806 : IF (use_virial) CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
585 :
586 : ! integrate potential with auxiliary basis function derivatives
587 10806 : NULLIFY (rs_v)
588 10806 : NULLIFY (rs_descs)
589 10806 : CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
590 10806 : CALL potential_pw2rs(rs_v, rho_r, pw_env_sub)
591 :
592 10806 : offset = 0
593 44547 : DO iatom = 1, SIZE(kind_of)
594 33741 : ikind = kind_of(iatom)
595 33741 : atom_a = atom_of_kind(iatom)
596 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
597 33741 : basis_type="RI_AUX")
598 :
599 33741 : first_sgfa => basis_set_a%first_sgf
600 33741 : la_max => basis_set_a%lmax
601 33741 : la_min => basis_set_a%lmin
602 33741 : npgfa => basis_set_a%npgf
603 33741 : nseta = basis_set_a%nset
604 33741 : nsgfa => basis_set_a%nsgf_set
605 33741 : rpgfa => basis_set_a%pgf_radius
606 33741 : set_radius_a => basis_set_a%set_radius
607 33741 : sphi_a => basis_set_a%sphi
608 33741 : zeta => basis_set_a%zet
609 :
610 33741 : ra(:) = pbc(particle_set(iatom)%r, cell)
611 :
612 33741 : force_a(:) = 0.0_dp
613 33741 : force_b(:) = 0.0_dp
614 33741 : IF (use_virial) THEN
615 7143 : my_virial_a = 0.0_dp
616 7143 : my_virial_b = 0.0_dp
617 : END IF
618 :
619 354622 : DO iset = 1, nseta
620 320881 : ncoa = npgfa(iset)*ncoset(la_max(iset))
621 320881 : sgfa = first_sgfa(1, iset)
622 :
623 962643 : ALLOCATE (I_tmp2(ncoa, 1))
624 2274253 : I_tmp2 = 0.0_dp
625 962643 : ALLOCATE (I_ab(nsgfa(iset), 1))
626 1560462 : I_ab = 0.0_dp
627 641762 : ALLOCATE (pab(ncoa, 1))
628 2274253 : pab = 0.0_dp
629 :
630 320881 : skip_shell = .TRUE.
631 1239581 : DO iorb = 1, nsgfa(iset)
632 1239581 : IF (iorb + offset == LLL) THEN
633 10806 : I_ab(iorb, 1) = 1.0_dp
634 10806 : skip_shell = .FALSE.
635 : END IF
636 : END DO
637 :
638 320881 : IF (skip_shell) THEN
639 310075 : offset = offset + nsgfa(iset)
640 310075 : DEALLOCATE (I_tmp2)
641 310075 : DEALLOCATE (I_ab)
642 310075 : DEALLOCATE (pab)
643 310075 : CYCLE
644 : END IF
645 :
646 : CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
647 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
648 : I_ab(1, 1), SIZE(I_ab, 1), &
649 10806 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
650 10806 : DEALLOCATE (I_ab)
651 :
652 32832 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, MINVAL(zeta(:, iset)))
653 10806 : map_it_here = .FALSE.
654 43224 : use_subpatch = .NOT. ALL(rs_v(igrid_level)%desc%perd == 1)
655 :
656 10806 : IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
657 20940 : DO ipgf = 1, npgfa(iset)
658 10632 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
659 10632 : na2 = ipgf*ncoset(la_max(iset))
660 10632 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))
661 :
662 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
663 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
664 : zetp=zeta(ipgf, iset), &
665 : eps=dft_control%qs_control%eps_gvg_rspace, &
666 10632 : prefactor=1.0_dp, cutoff=1.0_dp)
667 :
668 : CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), &
669 : lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, &
670 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
671 : rsgrid=rs_v(igrid_level), &
672 : hab=I_tmp2, &
673 : pab=pab, &
674 : o1=na1 - 1, &
675 : o2=0, &
676 : radius=radius, &
677 : calculate_forces=.TRUE., &
678 : force_a=force_a, force_b=force_b, &
679 : use_virial=use_virial, &
680 : my_virial_a=my_virial_a, &
681 : my_virial_b=my_virial_b, &
682 : use_subpatch=use_subpatch, &
683 20940 : subpatch_pattern=0)
684 : END DO
685 : END IF
686 :
687 10806 : DEALLOCATE (I_tmp2)
688 10806 : DEALLOCATE (pab)
689 :
690 44547 : offset = offset + nsgfa(iset)
691 :
692 : END DO
693 :
694 : force(ikind)%rho_elec(:, atom_a) = &
695 134964 : force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b(:)
696 44547 : IF (use_virial) THEN
697 92859 : h_stress = h_stress + my_virial_a + my_virial_b
698 : END IF
699 : END DO
700 :
701 10806 : CALL timestop(handle)
702 :
703 21612 : END SUBROUTINE integrate_potential_forces_3c_2c
704 :
705 : ! **************************************************************************************************
706 : !> \brief Calculates stress tensor contribution from the operator
707 : !> \param rho_g_copy ...
708 : !> \param pot_g ...
709 : !> \param rho_g ...
710 : !> \param dvg ...
711 : !> \param h_stress ...
712 : !> \param potential_parameter ...
713 : !> \param para_env_sub ...
714 : ! **************************************************************************************************
715 4650 : SUBROUTINE virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
716 :
717 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_g_copy, pot_g
718 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g
719 : TYPE(pw_c1d_gs_type), INTENT(IN) :: dvg(3)
720 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
721 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
722 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
723 :
724 : CHARACTER(LEN=*), PARAMETER :: routineN = 'virial_gpw_potential'
725 :
726 : INTEGER :: alpha, beta, handle
727 : INTEGER, DIMENSION(3) :: comp
728 : REAL(KIND=dp) :: e_hartree
729 :
730 : ! add the volume contribution
731 4650 : CALL timeset(routineN, handle)
732 4650 : e_hartree = 0.0_dp
733 4650 : e_hartree = pw_integral_ab(rho_g_copy, pot_g)
734 :
735 18600 : DO alpha = 1, 3
736 13950 : comp = 0
737 13950 : comp(alpha) = 1
738 13950 : CALL pw_copy(pot_g, rho_g)
739 13950 : CALL pw_derive(rho_g, comp)
740 13950 : CALL factor_virial_gpw(rho_g, potential_parameter)
741 13950 : h_stress(alpha, alpha) = h_stress(alpha, alpha) - e_hartree/REAL(para_env_sub%num_pe, dp)
742 46500 : DO beta = alpha, 3
743 : h_stress(alpha, beta) = h_stress(alpha, beta) &
744 27900 : - 2.0_dp*pw_integral_ab(rho_g, dvg(beta))/fourpi/REAL(para_env_sub%num_pe, dp)
745 41850 : h_stress(beta, alpha) = h_stress(alpha, beta)
746 : END DO
747 : END DO
748 :
749 4650 : CALL timestop(handle)
750 4650 : END SUBROUTINE virial_gpw_potential
751 :
752 : ! **************************************************************************************************
753 : !> \brief Multiplies a potential in g space with the factor ln(V(g)/Vc(g))' with Vc(g) being the
754 : !> Fourier-transformed of the Coulomg potential
755 : !> \param pw ...
756 : !> \param potential_parameter parameters of potential V(g)
757 : ! **************************************************************************************************
758 13950 : SUBROUTINE factor_virial_gpw(pw, potential_parameter)
759 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
760 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
761 :
762 13950 : SELECT CASE (potential_parameter%potential_type)
763 : CASE (do_potential_coulomb)
764 1329 : RETURN
765 : CASE (do_potential_long)
766 1329 : CALL pw_log_deriv_gauss(pw, potential_parameter%omega)
767 : CASE (do_potential_short)
768 0 : CALL pw_log_deriv_compl_gauss(pw, potential_parameter%omega)
769 : CASE (do_potential_mix_cl)
770 : CALL pw_log_deriv_mix_cl(pw, potential_parameter%omega, &
771 249 : potential_parameter%scale_coulomb, potential_parameter%scale_longrange)
772 : CASE (do_potential_truncated)
773 0 : CALL pw_log_deriv_trunc(pw, potential_parameter%cutoff_radius)
774 : CASE (do_potential_id)
775 0 : CALL pw_zero(pw)
776 : CASE DEFAULT
777 13950 : CPABORT("Unknown potential type")
778 : END SELECT
779 :
780 : END SUBROUTINE factor_virial_gpw
781 :
782 : ! **************************************************************************************************
783 : !> \brief Integrate potential of RI function with RI function
784 : !> \param pw_env_sub ...
785 : !> \param pot_r ...
786 : !> \param kind_of ...
787 : !> \param atom_of_kind ...
788 : !> \param particle_set ...
789 : !> \param qs_kind_set ...
790 : !> \param G_PQ_local ...
791 : !> \param cell ...
792 : !> \param force ...
793 : !> \param use_virial ...
794 : !> \param h_stress ...
795 : !> \param para_env_sub ...
796 : !> \param dft_control ...
797 : ! **************************************************************************************************
798 11304 : SUBROUTINE integrate_potential(pw_env_sub, pot_r, kind_of, atom_of_kind, particle_set, qs_kind_set, &
799 11304 : G_PQ_local, cell, force, use_virial, h_stress, para_env_sub, dft_control)
800 :
801 : TYPE(pw_env_type), INTENT(IN), POINTER :: pw_env_sub
802 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pot_r
803 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of, atom_of_kind
804 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
805 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
806 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: G_PQ_local
807 : TYPE(cell_type), INTENT(IN), POINTER :: cell
808 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
809 : LOGICAL, INTENT(IN) :: use_virial
810 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
811 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
812 : TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
813 :
814 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential'
815 :
816 : INTEGER :: atom_a, handle, iatom, igrid_level, &
817 : ikind, ipgf, iset, na1, na2, ncoa, &
818 : nseta, offset, sgfa
819 11304 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
820 11304 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
821 : LOGICAL :: use_subpatch
822 : REAL(KIND=dp) :: radius
823 11304 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab
824 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
825 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
826 11304 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
827 11304 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, pab, rpgfa, sphi_a, zeta
828 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
829 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
830 11304 : POINTER :: rs_descs
831 11304 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
832 :
833 11304 : CALL timeset(routineN, handle)
834 :
835 11304 : NULLIFY (rs_v)
836 11304 : NULLIFY (rs_descs)
837 11304 : CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
838 11304 : CALL potential_pw2rs(rs_v, pot_r, pw_env_sub)
839 :
840 11304 : offset = 0
841 46539 : DO iatom = 1, SIZE(kind_of)
842 35235 : ikind = kind_of(iatom)
843 35235 : atom_a = atom_of_kind(iatom)
844 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
845 35235 : basis_type="RI_AUX")
846 :
847 35235 : first_sgfa => basis_set_a%first_sgf
848 35235 : la_max => basis_set_a%lmax
849 35235 : la_min => basis_set_a%lmin
850 35235 : npgfa => basis_set_a%npgf
851 35235 : nseta = basis_set_a%nset
852 35235 : nsgfa => basis_set_a%nsgf_set
853 35235 : rpgfa => basis_set_a%pgf_radius
854 35235 : set_radius_a => basis_set_a%set_radius
855 35235 : sphi_a => basis_set_a%sphi
856 35235 : zeta => basis_set_a%zet
857 :
858 35235 : ra(:) = pbc(particle_set(iatom)%r, cell)
859 :
860 35235 : force_a(:) = 0.0_dp
861 35235 : force_b(:) = 0.0_dp
862 35235 : IF (use_virial) THEN
863 7641 : my_virial_a = 0.0_dp
864 7641 : my_virial_b = 0.0_dp
865 : END IF
866 :
867 370558 : DO iset = 1, nseta
868 335323 : ncoa = npgfa(iset)*ncoset(la_max(iset))
869 335323 : sgfa = first_sgfa(1, iset)
870 :
871 1005969 : ALLOCATE (I_tmp2(ncoa, 1))
872 2376841 : I_tmp2 = 0.0_dp
873 1005969 : ALLOCATE (I_ab(nsgfa(iset), 1))
874 1630680 : I_ab = 0.0_dp
875 670646 : ALLOCATE (pab(ncoa, 1))
876 2376841 : pab = 0.0_dp
877 :
878 1295357 : I_ab(1:nsgfa(iset), 1) = -4.0_dp*G_PQ_local(offset + 1:offset + nsgfa(iset))
879 :
880 : CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
881 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
882 : I_ab(1, 1), SIZE(I_ab, 1), &
883 335323 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
884 :
885 1630680 : I_ab = 0.0_dp
886 :
887 1007625 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, MINVAL(zeta(:, iset)))
888 1341292 : use_subpatch = .NOT. ALL(rs_v(igrid_level)%desc%perd == 1)
889 :
890 335323 : IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
891 642590 : DO ipgf = 1, npgfa(iset)
892 321709 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
893 321709 : na2 = ipgf*ncoset(la_max(iset))
894 321709 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))
895 :
896 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
897 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
898 : zetp=zeta(ipgf, iset), &
899 : eps=dft_control%qs_control%eps_gvg_rspace, &
900 321709 : prefactor=1.0_dp, cutoff=1.0_dp)
901 :
902 : CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), &
903 : lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, &
904 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
905 : rsgrid=rs_v(igrid_level), &
906 : hab=I_tmp2, &
907 : pab=pab, &
908 : o1=na1 - 1, &
909 : o2=0, &
910 : radius=radius, &
911 : calculate_forces=.TRUE., &
912 : force_a=force_a, force_b=force_b, &
913 : use_virial=use_virial, &
914 : my_virial_a=my_virial_a, &
915 : my_virial_b=my_virial_b, &
916 : use_subpatch=use_subpatch, &
917 642590 : subpatch_pattern=0)
918 :
919 : END DO
920 :
921 : END IF
922 :
923 335323 : DEALLOCATE (I_tmp2)
924 335323 : DEALLOCATE (I_ab)
925 335323 : DEALLOCATE (pab)
926 :
927 370558 : offset = offset + nsgfa(iset)
928 :
929 : END DO
930 :
931 : force(ikind)%rho_elec(:, atom_a) = &
932 140940 : force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b
933 46539 : IF (use_virial) THEN
934 99333 : h_stress = h_stress + my_virial_a + my_virial_b
935 : END IF
936 : END DO
937 :
938 11304 : CALL timestop(handle)
939 :
940 22608 : END SUBROUTINE
941 :
942 : ! **************************************************************************************************
943 : !> \brief Prepares GPW calculation for RI-MP2/RI-RPA
944 : !> \param qs_env ...
945 : !> \param dft_control ...
946 : !> \param e_cutoff_old ...
947 : !> \param cutoff_old ...
948 : !> \param relative_cutoff_old ...
949 : !> \param para_env_sub ...
950 : !> \param pw_env_sub ...
951 : !> \param auxbas_pw_pool ...
952 : !> \param poisson_env ...
953 : !> \param task_list_sub ...
954 : !> \param rho_r ...
955 : !> \param rho_g ...
956 : !> \param pot_g ...
957 : !> \param psi_L ...
958 : !> \param sab_orb_sub ...
959 : ! **************************************************************************************************
960 1022 : SUBROUTINE prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
961 : auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
962 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
963 : TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
964 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
965 : INTENT(OUT) :: e_cutoff_old
966 : REAL(KIND=dp), INTENT(OUT) :: cutoff_old, relative_cutoff_old
967 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
968 : TYPE(pw_env_type), POINTER :: pw_env_sub
969 : TYPE(pw_pool_type), INTENT(IN), POINTER :: auxbas_pw_pool
970 : TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
971 : TYPE(task_list_type), POINTER :: task_list_sub
972 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: rho_r
973 : TYPE(pw_c1d_gs_type), INTENT(OUT) :: rho_g, pot_g
974 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: psi_L
975 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
976 : INTENT(IN), POINTER :: sab_orb_sub
977 :
978 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_gpw'
979 :
980 : INTEGER :: handle, i_multigrid, n_multigrid
981 : LOGICAL :: skip_load_balance_distributed
982 : REAL(KIND=dp) :: progression_factor
983 : TYPE(qs_ks_env_type), POINTER :: ks_env
984 :
985 1022 : CALL timeset(routineN, handle)
986 :
987 1022 : CALL get_qs_env(qs_env, dft_control=dft_control, ks_env=ks_env)
988 :
989 : ! hack hack hack XXXXXXXXXXXXXXX rebuilds the pw_en with the new cutoffs
990 1022 : progression_factor = dft_control%qs_control%progression_factor
991 1022 : n_multigrid = SIZE(dft_control%qs_control%e_cutoff)
992 3066 : ALLOCATE (e_cutoff_old(n_multigrid))
993 5110 : e_cutoff_old(:) = dft_control%qs_control%e_cutoff
994 1022 : cutoff_old = dft_control%qs_control%cutoff
995 :
996 1022 : dft_control%qs_control%cutoff = qs_env%mp2_env%mp2_gpw%cutoff*0.5_dp
997 1022 : dft_control%qs_control%e_cutoff(1) = dft_control%qs_control%cutoff
998 4088 : DO i_multigrid = 2, n_multigrid
999 : dft_control%qs_control%e_cutoff(i_multigrid) = dft_control%qs_control%e_cutoff(i_multigrid - 1) &
1000 4088 : /progression_factor
1001 : END DO
1002 :
1003 1022 : relative_cutoff_old = dft_control%qs_control%relative_cutoff
1004 1022 : dft_control%qs_control%relative_cutoff = qs_env%mp2_env%mp2_gpw%relative_cutoff*0.5_dp
1005 :
1006 : ! a pw_env
1007 1022 : NULLIFY (pw_env_sub)
1008 1022 : CALL pw_env_create(pw_env_sub)
1009 1022 : CALL pw_env_rebuild(pw_env_sub, qs_env, para_env_sub)
1010 :
1011 : CALL pw_env_get(pw_env_sub, auxbas_pw_pool=auxbas_pw_pool, &
1012 1022 : poisson_env=poisson_env)
1013 : ! hack hack hack XXXXXXXXXXXXXXX
1014 :
1015 : ! now we need a task list, hard code skip_load_balance_distributed
1016 1022 : NULLIFY (task_list_sub)
1017 1022 : skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
1018 1022 : CALL allocate_task_list(task_list_sub)
1019 : CALL generate_qs_task_list(ks_env, task_list_sub, &
1020 : reorder_rs_grid_ranks=.TRUE., soft_valid=.FALSE., &
1021 : skip_load_balance_distributed=skip_load_balance_distributed, &
1022 1022 : pw_env_external=pw_env_sub, sab_orb_external=sab_orb_sub)
1023 :
1024 : ! get some of the grids ready
1025 1022 : CALL auxbas_pw_pool%create_pw(rho_r)
1026 1022 : CALL auxbas_pw_pool%create_pw(rho_g)
1027 1022 : CALL auxbas_pw_pool%create_pw(pot_g)
1028 1022 : CALL auxbas_pw_pool%create_pw(psi_L)
1029 :
1030 : ! run the FFT once, to set up buffers and to take into account the memory
1031 1022 : CALL pw_zero(rho_r)
1032 1022 : CALL pw_transfer(rho_r, rho_g)
1033 :
1034 1022 : CALL timestop(handle)
1035 :
1036 1022 : END SUBROUTINE prepare_gpw
1037 :
1038 : ! **************************************************************************************************
1039 : !> \brief Cleanup GPW integration for RI-MP2/RI-RPA
1040 : !> \param qs_env ...
1041 : !> \param e_cutoff_old ...
1042 : !> \param cutoff_old ...
1043 : !> \param relative_cutoff_old ...
1044 : !> \param para_env_sub ...
1045 : !> \param pw_env_sub ...
1046 : !> \param task_list_sub ...
1047 : !> \param auxbas_pw_pool ...
1048 : !> \param rho_r ...
1049 : !> \param rho_g ...
1050 : !> \param pot_g ...
1051 : !> \param psi_L ...
1052 : ! **************************************************************************************************
1053 1022 : SUBROUTINE cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
1054 : task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
1055 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1056 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1057 : INTENT(IN) :: e_cutoff_old
1058 : REAL(KIND=dp), INTENT(IN) :: cutoff_old, relative_cutoff_old
1059 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
1060 : TYPE(pw_env_type), POINTER :: pw_env_sub
1061 : TYPE(task_list_type), POINTER :: task_list_sub
1062 : TYPE(pw_pool_type), INTENT(IN), POINTER :: auxbas_pw_pool
1063 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_r
1064 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_g, pot_g
1065 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: psi_L
1066 :
1067 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cleanup_gpw'
1068 :
1069 : INTEGER :: handle
1070 : TYPE(dft_control_type), POINTER :: dft_control
1071 :
1072 1022 : CALL timeset(routineN, handle)
1073 :
1074 : ! and now free the whole lot
1075 1022 : CALL auxbas_pw_pool%give_back_pw(rho_r)
1076 1022 : CALL auxbas_pw_pool%give_back_pw(rho_g)
1077 1022 : CALL auxbas_pw_pool%give_back_pw(pot_g)
1078 1022 : CALL auxbas_pw_pool%give_back_pw(psi_L)
1079 :
1080 1022 : CALL deallocate_task_list(task_list_sub)
1081 1022 : CALL pw_env_release(pw_env_sub, para_env=para_env_sub)
1082 :
1083 1022 : CALL get_qs_env(qs_env, dft_control=dft_control)
1084 :
1085 : ! restore the initial value of the cutoff
1086 5110 : dft_control%qs_control%e_cutoff = e_cutoff_old
1087 1022 : dft_control%qs_control%cutoff = cutoff_old
1088 1022 : dft_control%qs_control%relative_cutoff = relative_cutoff_old
1089 :
1090 1022 : CALL timestop(handle)
1091 :
1092 1022 : END SUBROUTINE cleanup_gpw
1093 :
1094 : ! **************************************************************************************************
1095 : !> \brief Calculates potential from a given density in g-space
1096 : !> \param pot_r on output: potential in r space
1097 : !> \param rho_g on input: rho in g space
1098 : !> \param poisson_env ...
1099 : !> \param pot_g on output: potential in g space
1100 : !> \param potential_parameter Potential parameters, if not provided, assume Coulomb potential
1101 : !> \param dvg in output: first derivatives of the corresponding Coulomb potential
1102 : !> \param no_transfer whether NOT to transform potential from g-space to r-space (default: do it)
1103 : ! **************************************************************************************************
1104 58222 : SUBROUTINE calc_potential_gpw(pot_r, rho_g, poisson_env, pot_g, potential_parameter, dvg, no_transfer)
1105 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pot_r
1106 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_g
1107 : TYPE(pw_poisson_type), INTENT(IN), POINTER :: poisson_env
1108 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pot_g
1109 : TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
1110 : TYPE(pw_c1d_gs_type), DIMENSION(3), &
1111 : INTENT(INOUT), OPTIONAL :: dvg
1112 : LOGICAL, INTENT(IN), OPTIONAL :: no_transfer
1113 :
1114 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_potential_gpw'
1115 :
1116 : INTEGER :: comp(3), handle, idir, my_potential_type
1117 : LOGICAL :: my_transfer
1118 :
1119 58222 : CALL timeset(routineN, handle)
1120 :
1121 58222 : my_potential_type = do_potential_coulomb
1122 58222 : IF (PRESENT(potential_parameter)) THEN
1123 58222 : my_potential_type = potential_parameter%potential_type
1124 : END IF
1125 :
1126 58222 : my_transfer = .TRUE.
1127 58222 : IF (PRESENT(no_transfer)) my_transfer = .NOT. no_transfer
1128 :
1129 : ! in case we do Coulomb metric for RI, we need the Coulomb operator, but for RI with the
1130 : ! overlap metric, we do not need the Coulomb operator
1131 58222 : IF (my_potential_type /= do_potential_id) THEN
1132 57226 : IF (PRESENT(dvg)) THEN
1133 2491 : CALL pw_poisson_solve(poisson_env, rho_g, vhartree=pot_g, dvhartree=dvg)
1134 : ELSE
1135 54735 : CALL pw_poisson_solve(poisson_env, rho_g, vhartree=pot_g)
1136 : END IF
1137 57226 : IF (my_potential_type == do_potential_long) CALL pw_gauss_damp(pot_g, potential_parameter%omega)
1138 57226 : IF (my_potential_type == do_potential_short) CALL pw_compl_gauss_damp(pot_g, potential_parameter%omega)
1139 57226 : IF (my_potential_type == do_potential_truncated) CALL pw_truncated(pot_g, potential_parameter%cutoff_radius)
1140 57226 : IF (my_potential_type == do_potential_mix_cl) CALL pw_gauss_damp_mix(pot_g, potential_parameter%omega, &
1141 : potential_parameter%scale_coulomb, &
1142 332 : potential_parameter%scale_longrange)
1143 57226 : IF (my_transfer) CALL pw_transfer(pot_g, pot_r)
1144 : ELSE
1145 : ! If we use an overlap metric, make sure to use the correct potential=density on output
1146 996 : CALL pw_copy(rho_g, pot_g)
1147 996 : IF (PRESENT(dvg)) THEN
1148 0 : DO idir = 1, 3
1149 0 : CALL pw_copy(pot_g, dvg(idir))
1150 0 : comp = 0
1151 0 : comp(idir) = 1
1152 0 : CALL pw_derive(dvg(idir), comp)
1153 : END DO
1154 : END IF
1155 996 : IF (my_transfer) CALL pw_transfer(rho_g, pot_r)
1156 : END IF
1157 :
1158 55731 : IF (my_transfer) CALL pw_scale(pot_r, pot_r%pw_grid%dvol)
1159 58222 : CALL timestop(handle)
1160 :
1161 58222 : END SUBROUTINE calc_potential_gpw
1162 :
1163 : END MODULE mp2_eri_gpw
|
