Line data Source code
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 : ! **************************************************************************************************
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_fm_types, ONLY: cp_fm_type
21 : USE dbcsr_api, ONLY: dbcsr_p_type,&
22 : dbcsr_set
23 : USE gaussian_gridlevels, ONLY: gaussian_gridlevel
24 : USE input_constants, ONLY: do_potential_coulomb,&
25 : do_potential_id,&
26 : do_potential_long,&
27 : do_potential_mix_cl,&
28 : do_potential_short,&
29 : do_potential_truncated
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,&
37 : pw_env_rebuild
38 : USE pw_env_types, ONLY: pw_env_get,&
39 : pw_env_release,&
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, &
43 : pw_log_deriv_compl_gauss, pw_log_deriv_gauss, pw_log_deriv_mix_cl, pw_log_deriv_trunc, &
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
50 : USE qs_collocate_density, ONLY: calculate_rho_elec,&
51 : calculate_wavefunction,&
52 : collocate_single_gaussian
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,&
59 : get_qs_kind_set,&
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
63 : USE realspace_grid_types, ONLY: map_gaussian_here,&
64 : realspace_grid_desc_p_type,&
65 : realspace_grid_type
66 : USE rs_pw_interface, ONLY: potential_pw2rs
67 : USE task_list_methods, ONLY: generate_qs_task_list
68 : USE task_list_types, ONLY: allocate_task_list,&
69 : deallocate_task_list,&
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 :
81 : PUBLIC :: mp2_eri_2c_integrate_gpw, mp2_eri_3c_integrate_gpw, calc_potential_gpw, cleanup_gpw, prepare_gpw, &
82 : integrate_potential_forces_2c, integrate_potential_forces_3c_1c, integrate_potential_forces_3c_2c, &
83 : virial_gpw_potential
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 27358 : SUBROUTINE mp2_eri_3c_integrate_gpw(mo_coeff, psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
108 13679 : 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 13679 : 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 13679 : external_vector=external_vector)
143 :
144 13679 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
145 :
146 : ! and finally (K|mu nu)
147 13679 : 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 13679 : pw_env_external=pw_env_sub, task_list_external=task_list_sub)
151 :
152 13679 : CALL timestop(handle)
153 :
154 13679 : 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 368 : SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, &
169 368 : 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 368 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: offset_2d
185 368 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
186 368 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
187 : LOGICAL :: map_it_here, use_subpatch
188 : REAL(KIND=dp) :: cutoff_old, radius, relative_cutoff_old
189 368 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old
190 368 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab
191 : REAL(KIND=dp), DIMENSION(3) :: ra
192 368 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
193 368 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, rpgfa, sphi_a, zeta
194 368 : 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 368 : 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 368 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
205 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
206 368 : POINTER :: rs_descs
207 368 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
208 : TYPE(task_list_type), POINTER :: task_list_sub
209 :
210 368 : 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 368 : auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
214 :
215 368 : 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 1378448 : L_local_col = 0.0_dp
218 :
219 368 : i_counter = 0
220 16557 : DO LLL = my_group_L_start, my_group_L_end
221 16189 : 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 16189 : required_function=LLL, basis_type="RI_AUX")
227 :
228 16189 : CALL timeset(routineN//"_pot_lm", handle2)
229 :
230 16189 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
231 :
232 16189 : NULLIFY (rs_v)
233 16189 : NULLIFY (rs_descs)
234 16189 : CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
235 16189 : CALL potential_pw2rs(rs_v, rho_r, pw_env_sub)
236 :
237 16189 : CALL timestop(handle2)
238 :
239 16189 : offset = 0
240 : ! Prepare offset ahead of OMP parallel loop
241 16189 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxnset=max_nseta, basis_type="RI_AUX")
242 64756 : ALLOCATE (offset_2d(max_nseta, natom))
243 :
244 66582 : DO iatom = 1, natom
245 50393 : ikind = kind_of(iatom)
246 50393 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
247 50393 : nseta = basis_set_a%nset
248 50393 : nsgfa => basis_set_a%nsgf_set
249 540937 : DO iset = 1, nseta
250 474355 : offset = offset + nsgfa(iset)
251 524748 : 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 16189 : !$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 48935 : DEALLOCATE (offset_2d)
339 :
340 : END DO
341 :
342 2756528 : 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 368 : task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
346 :
347 368 : CALL timestop(handle)
348 :
349 736 : 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 32376 : 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 10792 : 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 10792 : CALL timeset(routineN, handle)
416 :
417 : ! calculate potential associated to the single aux function
418 10792 : CALL timeset(routineN//"_wf_pot", handle2)
419 : ! pseudo psi_L
420 10792 : CALL pw_zero(rho_r)
421 10792 : 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 10792 : pw_env_sub, required_function=LLL, basis_type='RI_AUX')
425 10792 : IF (use_virial) THEN
426 2145 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter, dvg)
427 : ELSE
428 8647 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
429 : END IF
430 10792 : CALL timestop(handle2)
431 :
432 10792 : IF (use_virial) THEN
433 : ! make a copy of the density in G space
434 2145 : 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 2145 : CALL pw_zero(psi_L)
441 2145 : 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 188146 : external_vector=0.5_dp*factor*G_PQ_local)
446 : ! transfer to reciprocal space and calculate potential
447 2145 : 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 2145 : 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 926972 : -0.25_dp*factor*G_PQ_local, cell, force, use_virial, h_stress, para_env_sub, dft_control)
456 :
457 10792 : CALL timestop(handle)
458 10792 : 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 10294 : 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 10294 : CALL timeset(routineN, handle)
485 :
486 : ! integrate the potential of the single gaussian and update
487 : ! 3-center forces
488 10294 : 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 10294 : task_list_external=task_list_sub)
493 :
494 10294 : CALL timestop(handle)
495 10294 : 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 10294 : 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 10294 : 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 10294 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
559 10294 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
560 : LOGICAL :: map_it_here, skip_shell, use_subpatch
561 : REAL(KIND=dp) :: radius
562 10294 : 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 10294 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
566 10294 : 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 10294 : POINTER :: rs_descs
570 10294 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
571 :
572 10294 : CALL timeset(routineN, handle)
573 :
574 : ! put the gamma density on grid
575 10294 : CALL timeset(routineN//"_Gpot", handle2)
576 10294 : CALL pw_zero(rho_r)
577 10294 : 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 10294 : ks_env=ks_env)
584 : ! calculate associated hartree potential
585 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
586 10294 : CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_parameter)
587 10294 : CALL timestop(handle2)
588 :
589 10294 : 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 10294 : NULLIFY (rs_v)
593 10294 : NULLIFY (rs_descs)
594 10294 : CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
595 10294 : CALL potential_pw2rs(rs_v, rho_r, pw_env_sub)
596 :
597 10294 : offset = 0
598 42402 : DO iatom = 1, SIZE(kind_of)
599 32108 : ikind = kind_of(iatom)
600 32108 : atom_a = atom_of_kind(iatom)
601 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
602 32108 : basis_type="RI_AUX")
603 :
604 32108 : first_sgfa => basis_set_a%first_sgf
605 32108 : la_max => basis_set_a%lmax
606 32108 : la_min => basis_set_a%lmin
607 32108 : npgfa => basis_set_a%npgf
608 32108 : nseta = basis_set_a%nset
609 32108 : nsgfa => basis_set_a%nsgf_set
610 32108 : rpgfa => basis_set_a%pgf_radius
611 32108 : set_radius_a => basis_set_a%set_radius
612 32108 : sphi_a => basis_set_a%sphi
613 32108 : zeta => basis_set_a%zet
614 :
615 32108 : ra(:) = pbc(particle_set(iatom)%r, cell)
616 :
617 32108 : force_a(:) = 0.0_dp
618 32108 : force_b(:) = 0.0_dp
619 32108 : IF (use_virial) THEN
620 6506 : my_virial_a = 0.0_dp
621 6506 : my_virial_b = 0.0_dp
622 : END IF
623 :
624 337559 : DO iset = 1, nseta
625 305451 : ncoa = npgfa(iset)*ncoset(la_max(iset))
626 305451 : sgfa = first_sgfa(1, iset)
627 :
628 916353 : ALLOCATE (I_tmp2(ncoa, 1))
629 2165580 : I_tmp2 = 0.0_dp
630 916353 : ALLOCATE (I_ab(nsgfa(iset), 1))
631 1485748 : I_ab = 0.0_dp
632 916353 : ALLOCATE (pab(ncoa, 1))
633 2165580 : pab = 0.0_dp
634 :
635 305451 : skip_shell = .TRUE.
636 1180297 : DO iorb = 1, nsgfa(iset)
637 1180297 : IF (iorb + offset == LLL) THEN
638 10294 : I_ab(iorb, 1) = 1.0_dp
639 10294 : skip_shell = .FALSE.
640 : END IF
641 : END DO
642 :
643 305451 : IF (skip_shell) THEN
644 295157 : offset = offset + nsgfa(iset)
645 295157 : DEALLOCATE (I_tmp2)
646 295157 : DEALLOCATE (I_ab)
647 295157 : DEALLOCATE (pab)
648 295157 : 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 10294 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
655 10294 : DEALLOCATE (I_ab)
656 :
657 31296 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, MINVAL(zeta(:, iset)))
658 10294 : map_it_here = .FALSE.
659 41176 : use_subpatch = .NOT. ALL(rs_v(igrid_level)%desc%perd == 1)
660 :
661 10294 : IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
662 19916 : DO ipgf = 1, npgfa(iset)
663 10120 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
664 10120 : na2 = ipgf*ncoset(la_max(iset))
665 10120 : 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 10120 : 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 19916 : subpatch_pattern=0)
689 : END DO
690 : END IF
691 :
692 10294 : DEALLOCATE (I_tmp2)
693 10294 : DEALLOCATE (pab)
694 :
695 42402 : offset = offset + nsgfa(iset)
696 :
697 : END DO
698 :
699 : force(ikind)%rho_elec(:, atom_a) = &
700 128432 : force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b(:)
701 42402 : IF (use_virial) THEN
702 84578 : h_stress = h_stress + my_virial_a + my_virial_b
703 : END IF
704 : END DO
705 :
706 10294 : CALL timestop(handle)
707 :
708 20588 : 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 4290 : 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 4290 : CALL timeset(routineN, handle)
737 4290 : e_hartree = 0.0_dp
738 4290 : e_hartree = pw_integral_ab(rho_g_copy, pot_g)
739 :
740 17160 : DO alpha = 1, 3
741 12870 : comp = 0
742 12870 : comp(alpha) = 1
743 12870 : CALL pw_copy(pot_g, rho_g)
744 12870 : CALL pw_derive(rho_g, comp)
745 12870 : CALL factor_virial_gpw(rho_g, potential_parameter)
746 12870 : h_stress(alpha, alpha) = h_stress(alpha, alpha) - e_hartree/REAL(para_env_sub%num_pe, dp)
747 42900 : DO beta = alpha, 3
748 : h_stress(alpha, beta) = h_stress(alpha, beta) &
749 25740 : - 2.0_dp*pw_integral_ab(rho_g, dvg(beta))/fourpi/REAL(para_env_sub%num_pe, dp)
750 38610 : h_stress(beta, alpha) = h_stress(alpha, beta)
751 : END DO
752 : END DO
753 :
754 4290 : CALL timestop(handle)
755 4290 : 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 12870 : 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 12870 : SELECT CASE (potential_parameter%potential_type)
768 : CASE (do_potential_coulomb)
769 1329 : RETURN
770 : CASE (do_potential_long)
771 1329 : CALL pw_log_deriv_gauss(pw, potential_parameter%omega)
772 : CASE (do_potential_short)
773 0 : 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 249 : potential_parameter%scale_coulomb, potential_parameter%scale_longrange)
777 : CASE (do_potential_truncated)
778 0 : CALL pw_log_deriv_trunc(pw, potential_parameter%cutoff_radius)
779 : CASE (do_potential_id)
780 0 : CALL pw_zero(pw)
781 : CASE DEFAULT
782 12870 : 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 10792 : SUBROUTINE integrate_potential(pw_env_sub, pot_r, kind_of, atom_of_kind, particle_set, qs_kind_set, &
804 10792 : 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 10792 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
825 10792 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
826 : LOGICAL :: use_subpatch
827 : REAL(KIND=dp) :: radius
828 10792 : 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 10792 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
832 10792 : 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 10792 : POINTER :: rs_descs
836 10792 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
837 :
838 10792 : CALL timeset(routineN, handle)
839 :
840 10792 : NULLIFY (rs_v)
841 10792 : NULLIFY (rs_descs)
842 10792 : CALL pw_env_get(pw_env_sub, rs_descs=rs_descs, rs_grids=rs_v)
843 10792 : CALL potential_pw2rs(rs_v, pot_r, pw_env_sub)
844 :
845 10792 : offset = 0
846 44394 : DO iatom = 1, SIZE(kind_of)
847 33602 : ikind = kind_of(iatom)
848 33602 : atom_a = atom_of_kind(iatom)
849 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
850 33602 : basis_type="RI_AUX")
851 :
852 33602 : first_sgfa => basis_set_a%first_sgf
853 33602 : la_max => basis_set_a%lmax
854 33602 : la_min => basis_set_a%lmin
855 33602 : npgfa => basis_set_a%npgf
856 33602 : nseta = basis_set_a%nset
857 33602 : nsgfa => basis_set_a%nsgf_set
858 33602 : rpgfa => basis_set_a%pgf_radius
859 33602 : set_radius_a => basis_set_a%set_radius
860 33602 : sphi_a => basis_set_a%sphi
861 33602 : zeta => basis_set_a%zet
862 :
863 33602 : ra(:) = pbc(particle_set(iatom)%r, cell)
864 :
865 33602 : force_a(:) = 0.0_dp
866 33602 : force_b(:) = 0.0_dp
867 33602 : IF (use_virial) THEN
868 7004 : my_virial_a = 0.0_dp
869 7004 : my_virial_b = 0.0_dp
870 : END IF
871 :
872 353495 : DO iset = 1, nseta
873 319893 : ncoa = npgfa(iset)*ncoset(la_max(iset))
874 319893 : sgfa = first_sgfa(1, iset)
875 :
876 959679 : ALLOCATE (I_tmp2(ncoa, 1))
877 2268168 : I_tmp2 = 0.0_dp
878 959679 : ALLOCATE (I_ab(nsgfa(iset), 1))
879 1555966 : I_ab = 0.0_dp
880 959679 : ALLOCATE (pab(ncoa, 1))
881 2268168 : pab = 0.0_dp
882 :
883 1236073 : 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 319893 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
889 :
890 1555966 : I_ab = 0.0_dp
891 :
892 961335 : igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, MINVAL(zeta(:, iset)))
893 1279572 : use_subpatch = .NOT. ALL(rs_v(igrid_level)%desc%perd == 1)
894 :
895 319893 : IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, offset, para_env_sub%num_pe, para_env_sub%mepos)) THEN
896 611730 : DO ipgf = 1, npgfa(iset)
897 306279 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
898 306279 : na2 = ipgf*ncoset(la_max(iset))
899 306279 : 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 306279 : 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 611730 : subpatch_pattern=0)
923 :
924 : END DO
925 :
926 : END IF
927 :
928 319893 : DEALLOCATE (I_tmp2)
929 319893 : DEALLOCATE (I_ab)
930 319893 : DEALLOCATE (pab)
931 :
932 353495 : offset = offset + nsgfa(iset)
933 :
934 : END DO
935 :
936 : force(ikind)%rho_elec(:, atom_a) = &
937 134408 : force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b
938 44394 : IF (use_virial) THEN
939 91052 : h_stress = h_stress + my_virial_a + my_virial_b
940 : END IF
941 : END DO
942 :
943 10792 : CALL timestop(handle)
944 :
945 21584 : 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 950 : 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 950 : CALL timeset(routineN, handle)
991 :
992 950 : 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 950 : progression_factor = dft_control%qs_control%progression_factor
996 950 : n_multigrid = SIZE(dft_control%qs_control%e_cutoff)
997 2850 : ALLOCATE (e_cutoff_old(n_multigrid))
998 4750 : e_cutoff_old(:) = dft_control%qs_control%e_cutoff
999 950 : cutoff_old = dft_control%qs_control%cutoff
1000 :
1001 950 : dft_control%qs_control%cutoff = qs_env%mp2_env%mp2_gpw%cutoff*0.5_dp
1002 950 : dft_control%qs_control%e_cutoff(1) = dft_control%qs_control%cutoff
1003 3800 : 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 3800 : /progression_factor
1006 : END DO
1007 :
1008 950 : relative_cutoff_old = dft_control%qs_control%relative_cutoff
1009 950 : dft_control%qs_control%relative_cutoff = qs_env%mp2_env%mp2_gpw%relative_cutoff*0.5_dp
1010 :
1011 : ! a pw_env
1012 950 : NULLIFY (pw_env_sub)
1013 950 : CALL pw_env_create(pw_env_sub)
1014 950 : 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 950 : poisson_env=poisson_env)
1018 : ! hack hack hack XXXXXXXXXXXXXXX
1019 :
1020 : ! now we need a task list, hard code skip_load_balance_distributed
1021 950 : NULLIFY (task_list_sub)
1022 950 : skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
1023 950 : 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 950 : pw_env_external=pw_env_sub, sab_orb_external=sab_orb_sub)
1028 :
1029 : ! get some of the grids ready
1030 950 : CALL auxbas_pw_pool%create_pw(rho_r)
1031 950 : CALL auxbas_pw_pool%create_pw(rho_g)
1032 950 : CALL auxbas_pw_pool%create_pw(pot_g)
1033 950 : 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 950 : CALL pw_zero(rho_r)
1037 950 : CALL pw_transfer(rho_r, rho_g)
1038 :
1039 950 : CALL timestop(handle)
1040 :
1041 950 : 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 950 : 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 950 : CALL timeset(routineN, handle)
1078 :
1079 : ! and now free the whole lot
1080 950 : CALL auxbas_pw_pool%give_back_pw(rho_r)
1081 950 : CALL auxbas_pw_pool%give_back_pw(rho_g)
1082 950 : CALL auxbas_pw_pool%give_back_pw(pot_g)
1083 950 : CALL auxbas_pw_pool%give_back_pw(psi_L)
1084 :
1085 950 : CALL deallocate_task_list(task_list_sub)
1086 950 : CALL pw_env_release(pw_env_sub, para_env=para_env_sub)
1087 :
1088 950 : CALL get_qs_env(qs_env, dft_control=dft_control)
1089 :
1090 : ! restore the initial value of the cutoff
1091 4750 : dft_control%qs_control%e_cutoff = e_cutoff_old
1092 950 : dft_control%qs_control%cutoff = cutoff_old
1093 950 : dft_control%qs_control%relative_cutoff = relative_cutoff_old
1094 :
1095 950 : CALL timestop(handle)
1096 :
1097 950 : 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 54858 : 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 54858 : CALL timeset(routineN, handle)
1125 :
1126 54858 : my_potential_type = do_potential_coulomb
1127 54858 : IF (PRESENT(potential_parameter)) THEN
1128 54858 : my_potential_type = potential_parameter%potential_type
1129 : END IF
1130 :
1131 54858 : my_transfer = .TRUE.
1132 54858 : 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 54858 : IF (my_potential_type /= do_potential_id) THEN
1137 53862 : IF (PRESENT(dvg)) THEN
1138 2311 : CALL pw_poisson_solve(poisson_env, rho_g, vhartree=pot_g, dvhartree=dvg)
1139 : ELSE
1140 51551 : CALL pw_poisson_solve(poisson_env, rho_g, vhartree=pot_g)
1141 : END IF
1142 53862 : IF (my_potential_type == do_potential_long) CALL pw_gauss_damp(pot_g, potential_parameter%omega)
1143 53862 : IF (my_potential_type == do_potential_short) CALL pw_compl_gauss_damp(pot_g, potential_parameter%omega)
1144 53862 : IF (my_potential_type == do_potential_truncated) CALL pw_truncated(pot_g, potential_parameter%cutoff_radius)
1145 53862 : IF (my_potential_type == do_potential_mix_cl) CALL pw_gauss_damp_mix(pot_g, potential_parameter%omega, &
1146 : potential_parameter%scale_coulomb, &
1147 332 : potential_parameter%scale_longrange)
1148 53862 : 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 996 : CALL pw_copy(rho_g, pot_g)
1152 996 : IF (PRESENT(dvg)) THEN
1153 0 : DO idir = 1, 3
1154 0 : CALL pw_copy(pot_g, dvg(idir))
1155 0 : comp = 0
1156 0 : comp(idir) = 1
1157 0 : CALL pw_derive(dvg(idir), comp)
1158 : END DO
1159 : END IF
1160 996 : IF (my_transfer) CALL pw_transfer(rho_g, pot_r)
1161 : END IF
1162 :
1163 54858 : IF (my_transfer) CALL pw_scale(pot_r, pot_r%pw_grid%dvol)
1164 54858 : CALL timestop(handle)
1165 :
1166 54858 : END SUBROUTINE calc_potential_gpw
1167 :
1168 : END MODULE mp2_eri_gpw
|
