(git:1f9fd2c)
Loading...
Searching...
No Matches
hfx_pair_list_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines for optimizing load balance between processes in HFX calculations
10!> \par History
11!> 04.2008 created [Manuel Guidon]
12!> 11.2019 fixed initial value for potential_id (A. Bussy)
13!> \author Manuel Guidon
14! **************************************************************************************************
15MODULE hfx_pair_list_methods
16 USE cell_types, ONLY: cell_type,&
17 pbc
18 USE gamma, ONLY: fgamma => fgamma_0
19 USE hfx_types, ONLY: &
20 hfx_basis_type, hfx_block_range_type, hfx_cell_type, hfx_pgf_list, hfx_pgf_product_list, &
21 hfx_potential_type, hfx_screen_coeff_type, pair_list_type, pair_set_list_type
22 USE input_constants, ONLY: &
23 do_potential_tshpsc, do_potential_coulomb, do_potential_gaussian, do_potential_id, &
24 do_potential_long, do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, &
25 do_potential_short, do_potential_truncated
26 USE kinds, ONLY: dp
27 USE libint_wrapper, ONLY: prim_data_f_size
28 USE mathconstants, ONLY: pi
29 USE mp2_types, ONLY: pair_list_type_mp2
30 USE particle_types, ONLY: particle_type
31 USE t_c_g0, ONLY: t_c_g0_n
32 USE t_sh_p_s_c, ONLY: trunc_cs_poly_n20
33#include "./base/base_uses.f90"
34
35 IMPLICIT NONE
36 PRIVATE
37
38 PUBLIC :: build_pair_list, &
39 build_pair_list_mp2, &
40 build_pair_list_pgf, &
41 build_pgf_product_list, &
42 build_atomic_pair_list, &
43 pgf_product_list_size, &
44 build_pair_list_pbc_pgf, &
45 bra_t, allocate_bra
46
47 ! an initial estimate for the size of the product list
48 INTEGER, SAVE :: pgf_product_list_size = 128
49
50 ! Store screened primitive gaussian function pairs between two shells (A and B)
51 TYPE :: bra_t
52 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pgf_scr
53 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: pgf_idx, cell_idx
54 INTEGER :: cell_cnt = 0, pgf_cnt = 0
55 END TYPE bra_t
56
57CONTAINS
58
59! **************************************************************************************************
60!> \brief Allocates all arrays within a Bra_t structure.
61!>
62!> This subroutine ensures that a given `Bra_t` pointer has allocated memory
63!> for storing primitive Gaussian pair indices, screening data, and cell indices.
64!> Existing allocations are safely deallocated first.
65!>
66!> \param[in,out] bra Pointer to the Bra_t structure to allocate.
67!> \param[in] max_npgf Maximum number of primitives per shell.
68!> \param[in] max_ncell Maximum number of neighbor cells.
69!>
70!> \note The arrays are allocated as follows:
71!> - `pgf_idx(3, max_npgf^2 * max_ncell)`
72!> - `pgf_scr(5, max_npgf^2 * max_ncell)`
73!> - `cell_idx(3, max_ncell)`
74!>
75! **************************************************************************************************
76 SUBROUTINE allocate_bra(bra, max_npgf, max_ncell)
77 TYPE(bra_t), POINTER :: bra
78 INTEGER :: max_npgf, max_ncell
79
80 IF (ALLOCATED(bra%pgf_idx)) DEALLOCATE (bra%pgf_idx) ! TODO add if n_prm changed
81 IF (ALLOCATED(bra%pgf_scr)) DEALLOCATE (bra%pgf_scr) ! TODO add if n_prm changed
82 IF (ALLOCATED(bra%cell_idx)) DEALLOCATE (bra%cell_idx) ! TODO add if n_prm changed
83 ALLOCATE (bra%pgf_idx(3, max_npgf*max_npgf*max_ncell))
84 ALLOCATE (bra%pgf_scr(5, max_npgf*max_npgf*max_ncell))
85 ALLOCATE (bra%cell_idx(3, max_ncell))
86 END SUBROUTINE allocate_bra
87
88! **************************************************************************************************
89!> \brief Builds a screened list of primitives from centers A and B, intersecting with another shell
90!>
91!> This subroutine populates the `Bra_t` structure with indices and screening data
92!> for primitive Gaussian pairs (ipgf,jpgf,n3) that pass the Schwarz screening criterion
93!>
94!> \param[in] npgfa Number of primitives in shell A
95!> \param[in] npgfb Number of primitives in shell B
96!> \param[out] bra Pointer to Bra_t structure to populate
97!> \param[in] screen1 Screening coeffiecients for the AB pair based on most diffuse pgf
98!> \param[in] screen2 Screening coeffiecients for the CD pair based on most diffuse pgf
99!> \param[in] pgf Screening coeeficents for primitive gaussians
100!> \param[in] log10_pmax Log10 of maximum integral prefactor
101!> \param[in] log10_eps_schwarz Log10 of Schwarz screening threshold
102!> \param[in] ra Cartesian coordinates of center A
103!> \param[in] rb Cartesian coordinates of center B
104!> \param[out] nelements Total number of valid primitives found
105!> \param[in] neighbor_cells Array of lattice vectors
106!> \param[in] do_periodic Logical flag for periodic boundary conditions on B
107!>
108!> \note Loops over all cells, shifting the B shell position accordingly
109!> Only pairs satisfying the screening condition are stored.
110!>
111!> \note Each valid cell contributes a block of pairs stored contiguously in
112!> bra%pgf_idx, with cell metadata stored in bra%cell_idx
113!> \note screening primitive pairs (ipgf, jpgf) using both coarse and fine screening thresholds
114! **************************************************************************************************
115 SUBROUTINE build_pair_list_pbc_pgf(npgfa, npgfb, bra, screen1, screen2, pgf, &
116 log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
117 neighbor_cells, do_periodic)
118
119 INTEGER, INTENT(IN) :: npgfa, npgfb
120 TYPE(bra_t), POINTER :: bra
121 REAL(dp), INTENT(IN) :: screen1(2), screen2(2)
122 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
123 POINTER :: pgf
124 REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), &
125 rb(3)
126 INTEGER, INTENT(OUT) :: nelements
127 TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
128 LOGICAL, INTENT(IN) :: do_periodic
129
130 INTEGER :: cell_cnt, cell_idx, element_cnt, &
131 element_idx, element_off, ipgf, jpgf
132 REAL(dp) :: ab(3), im_b(3), pgf_max, rab2
133
134 element_idx = 0
135 element_off = 0
136 cell_cnt = 0
137 DO cell_idx = 1, SIZE(neighbor_cells)
138 ! move B to this cell while keeping A fixed
139 ! NOTE rb has already been moved by build_pair_list
140 ! so that when cell_r is (000) AB is in the pbc cell
141 IF (do_periodic) THEN
142 im_b = rb + neighbor_cells(cell_idx)%cell_r(:)
143 ELSE
144 im_b = rb
145 END IF
146 ab = ra - im_b
147 rab2 = ab(1)**2 + ab(2)**2 + ab(3)**2
148 ! First screening on the most diffuse AB gaussians along with the most diffuse CD gaussian
149 IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) cycle
150
151 element_cnt = 0
152 DO ipgf = 1, npgfa
153 DO jpgf = 1, npgfb
154
155 ! Second screening on this actual AB pair and the most diffuse CD from before
156 pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
157 IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) cycle
158 ! This pair passed screening. Add it to the bra.
159 element_idx = element_idx + 1
160 bra%pgf_idx(:, element_idx) = [ipgf, jpgf, cell_idx]
161 bra%pgf_scr(1, element_idx) = pgf_max
162 element_cnt = element_cnt + 1
163 END DO
164 END DO
165
166 ! If this cell produced any pair, add this cell to the bra
167 IF (element_cnt == 0) cycle
168 cell_cnt = cell_cnt + 1
169 bra%cell_idx(:, cell_cnt) = [cell_idx, element_cnt, element_off]
170 element_off = element_off + element_cnt
171
172 END DO ! cell
173 bra%cell_cnt = cell_cnt
174 nelements = element_idx
175 bra%pgf_cnt = nelements
176
177 END SUBROUTINE build_pair_list_pbc_pgf
178
179! **************************************************************************************************
180!> \brief ...
181!> \param list1 ...
182!> \param list2 ...
183!> \param product_list ...
184!> \param nproducts ...
185!> \param log10_pmax ...
186!> \param log10_eps_schwarz ...
187!> \param neighbor_cells ...
188!> \param cell ...
189!> \param potential_parameter ...
190!> \param m_max ...
191!> \param do_periodic ...
192! **************************************************************************************************
193 SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
194 log10_pmax, log10_eps_schwarz, neighbor_cells, &
195 cell, potential_parameter, m_max, do_periodic)
196
197 TYPE(hfx_pgf_list) :: list1, list2
198 TYPE(hfx_pgf_product_list), ALLOCATABLE, &
199 DIMENSION(:), INTENT(INOUT) :: product_list
200 INTEGER, INTENT(OUT) :: nproducts
201 REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
202 TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
203 TYPE(cell_type), POINTER :: cell
204 TYPE(hfx_potential_type) :: potential_parameter
205 INTEGER, INTENT(IN) :: m_max
206 LOGICAL, INTENT(IN) :: do_periodic
207
208 INTEGER :: i, j, k, l, nimages1, nimages2, tmp_i4
209 LOGICAL :: use_gamma
210 REAL(dp) :: c11(3), den, eta, etainv, factor, fm(prim_data_f_size), g(3), num, omega2, &
211 omega_corr, omega_corr2, p(3), pgf_max_1, pgf_max_2, pq(3), q(3), r, r1, r2, ra(3), &
212 rb(3), rc(3), rd(3), rho, rhoinv, rpq2, s1234, s1234a, s1234b, shift(3), ssss, t, &
213 temp(3), temp_cc(3), temp_dd(3), tmp, tmp_d(3), w(3), zeta1, zeta_c, zeta_d, zetapetainv
214 TYPE(hfx_pgf_product_list), ALLOCATABLE, &
215 DIMENSION(:) :: tmp_product_list
216
217 nimages1 = list1%nimages
218 nimages2 = list2%nimages
219 nproducts = 0
220 zeta1 = list1%zetapzetb
221 eta = list2%zetapzetb
222 etainv = list2%ZetaInv
223 zeta_c = list2%zeta
224 zeta_d = list2%zetb
225 temp_cc = 0.0_dp
226 temp_dd = 0.0_dp
227 DO i = 1, nimages1
228 p = list1%image_list(i)%P
229 r1 = list1%image_list(i)%R
230 s1234a = list1%image_list(i)%S1234
231 pgf_max_1 = list1%image_list(i)%pgf_max
232 ra = list1%image_list(i)%ra
233 rb = list1%image_list(i)%rb
234 DO j = 1, nimages2
235 pgf_max_2 = list2%image_list(j)%pgf_max
236 IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) cycle
237 q = list2%image_list(j)%P
238 r2 = list2%image_list(j)%R
239 s1234b = list2%image_list(j)%S1234
240 rc = list2%image_list(j)%ra
241 rd = list2%image_list(j)%rb
242
243 zetapetainv = zeta1 + eta
244 zetapetainv = 1.0_dp/zetapetainv
245 rho = zeta1*eta*zetapetainv
246 rhoinv = 1.0_dp/rho
247 s1234 = exp(s1234a + s1234b)
248 IF (do_periodic) THEN
249 temp = p - q
250 pq = pbc(temp, cell)
251 shift = -pq + temp
252 temp_cc = rc + shift
253 temp_dd = rd + shift
254 END IF
255
256 DO k = 1, SIZE(neighbor_cells)
257 IF (do_periodic) THEN
258 c11 = temp_cc + neighbor_cells(k)%cell_r(:)
259 tmp_d = temp_dd + neighbor_cells(k)%cell_r(:)
260 ELSE
261 c11 = rc
262 tmp_d = rd
263 END IF
264 q = (zeta_c*c11 + zeta_d*tmp_d)*etainv
265 rpq2 = (p(1) - q(1))**2 + (p(2) - q(2))**2 + (p(3) - q(3))**2
266 IF (potential_parameter%potential_type == do_potential_truncated .OR. &
267 potential_parameter%potential_type == do_potential_short .OR. &
268 potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
269 IF (rpq2 > (r1 + r2 + potential_parameter%cutoff_radius)**2) cycle
270 END IF
271 IF (potential_parameter%potential_type == do_potential_tshpsc) THEN
272 IF (rpq2 > (r1 + r2 + potential_parameter%cutoff_radius*2.0_dp)**2) cycle
273 END IF
274 nproducts = nproducts + 1
275
276 ! allocate size as needed,
277 ! updating the global size estimate to make this a rare event in longer simulations
278 IF (nproducts > SIZE(product_list)) THEN
279!$OMP ATOMIC READ
280 tmp_i4 = pgf_product_list_size
281 tmp_i4 = max(pgf_product_list_size, (3*nproducts + 1)/2)
282!$OMP ATOMIC WRITE
283 pgf_product_list_size = tmp_i4
284 ALLOCATE (tmp_product_list(SIZE(product_list)))
285 tmp_product_list(:) = product_list
286 DEALLOCATE (product_list)
287 ALLOCATE (product_list(tmp_i4))
288 product_list(1:SIZE(tmp_product_list)) = tmp_product_list
289 DEALLOCATE (tmp_product_list)
290 END IF
291
292 t = rho*rpq2
293 SELECT CASE (potential_parameter%potential_type)
294 CASE (do_potential_truncated)
295 r = potential_parameter%cutoff_radius*sqrt(rho)
296 CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, r, t, m_max)
297 IF (use_gamma) CALL fgamma(m_max, t, product_list(nproducts)%Fm(1))
298 factor = 2.0_dp*pi*rhoinv
299 CASE (do_potential_tshpsc)
300 r = potential_parameter%cutoff_radius*sqrt(rho)
301 product_list(nproducts)%Fm = 0.0_dp
302 CALL trunc_cs_poly_n20(product_list(nproducts)%Fm(1), r, t, m_max)
303 factor = 2.0_dp*pi*rhoinv
304 CASE (do_potential_coulomb)
305 CALL fgamma(m_max, t, product_list(nproducts)%Fm(1))
306 factor = 2.0_dp*pi*rhoinv
307 CASE (do_potential_short)
308 CALL fgamma(m_max, t, product_list(nproducts)%Fm(1))
309 omega2 = potential_parameter%omega**2
310 omega_corr2 = omega2/(omega2 + rho)
311 omega_corr = sqrt(omega_corr2)
312 t = t*omega_corr2
313 CALL fgamma(m_max, t, fm)
314 tmp = -omega_corr
315 DO l = 1, m_max + 1
316 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + fm(l)*tmp
317 tmp = tmp*omega_corr2
318 END DO
319 factor = 2.0_dp*pi*rhoinv
320 CASE (do_potential_long)
321 omega2 = potential_parameter%omega**2
322 omega_corr2 = omega2/(omega2 + rho)
323 omega_corr = sqrt(omega_corr2)
324 t = t*omega_corr2
325 CALL fgamma(m_max, t, product_list(nproducts)%Fm(1))
326 tmp = omega_corr
327 DO l = 1, m_max + 1
328 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*tmp
329 tmp = tmp*omega_corr2
330 END DO
331 factor = 2.0_dp*pi*rhoinv
332 CASE (do_potential_mix_cl)
333 CALL fgamma(m_max, t, product_list(nproducts)%Fm(1))
334 omega2 = potential_parameter%omega**2
335 omega_corr2 = omega2/(omega2 + rho)
336 omega_corr = sqrt(omega_corr2)
337 t = t*omega_corr2
338 CALL fgamma(m_max, t, fm)
339 tmp = omega_corr
340 DO l = 1, m_max + 1
341 product_list(nproducts)%Fm(l) = &
342 product_list(nproducts)%Fm(l)*potential_parameter%scale_coulomb &
343 + fm(l)*tmp*potential_parameter%scale_longrange
344 tmp = tmp*omega_corr2
345 END DO
346 factor = 2.0_dp*pi*rhoinv
347 CASE (do_potential_mix_cl_trunc)
348
349 ! truncated
350 r = potential_parameter%cutoff_radius*sqrt(rho)
351 CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, r, t, m_max)
352 IF (use_gamma) CALL fgamma(m_max, t, product_list(nproducts)%Fm(1))
353
354 ! Coulomb
355 CALL fgamma(m_max, t, fm)
356
357 DO l = 1, m_max + 1
358 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)* &
359 (potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
360 fm(l)*potential_parameter%scale_longrange
361 END DO
362
363 ! longrange
364 omega2 = potential_parameter%omega**2
365 omega_corr2 = omega2/(omega2 + rho)
366 omega_corr = sqrt(omega_corr2)
367 t = t*omega_corr2
368 CALL fgamma(m_max, t, fm)
369 tmp = omega_corr
370 DO l = 1, m_max + 1
371 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + fm(l)*tmp*potential_parameter%scale_longrange
372 tmp = tmp*omega_corr2
373 END DO
374 factor = 2.0_dp*pi*rhoinv
375
376 CASE (do_potential_gaussian)
377 omega2 = potential_parameter%omega**2
378 t = -omega2*t/(rho + omega2)
379 tmp = 1.0_dp
380 DO l = 1, m_max + 1
381 product_list(nproducts)%Fm(l) = exp(t)*tmp
382 tmp = tmp*omega2/(rho + omega2)
383 END DO
384 factor = (pi/(rho + omega2))**(1.5_dp)
385 CASE (do_potential_mix_lg)
386 omega2 = potential_parameter%omega**2
387 omega_corr2 = omega2/(omega2 + rho)
388 omega_corr = sqrt(omega_corr2)
389 t = t*omega_corr2
390 CALL fgamma(m_max, t, fm)
391 tmp = omega_corr*2.0_dp*pi*rhoinv*potential_parameter%scale_longrange
392 DO l = 1, m_max + 1
393 fm(l) = fm(l)*tmp
394 tmp = tmp*omega_corr2
395 END DO
396 t = rho*rpq2
397 t = -omega2*t/(rho + omega2)
398 tmp = (pi/(rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
399 DO l = 1, m_max + 1
400 product_list(nproducts)%Fm(l) = exp(t)*tmp + fm(l)
401 tmp = tmp*omega2/(rho + omega2)
402 END DO
403 factor = 1.0_dp
404 CASE (do_potential_id)
405 num = list1%zeta*list1%zetb
406 den = list1%zeta + list1%zetb
407 ssss = -num/den*sum((ra - rb)**2)
408
409 num = den*zeta_c
410 den = den + zeta_c
411 ssss = ssss - num/den*sum((p - rc)**2)
412
413 g(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + zeta_c*rc(:))/den
414 num = den*zeta_d
415 den = den + zeta_d
416 ssss = ssss - num/den*sum((g - rd)**2)
417
418 product_list(nproducts)%Fm(:) = exp(ssss)
419 factor = 1.0_dp
420 IF (s1234 > epsilon(0.0_dp)) factor = 1.0_dp/s1234
421 END SELECT
422
423 tmp = (pi*zetapetainv)**3
424 factor = factor*s1234*sqrt(tmp)
425
426 DO l = 1, m_max + 1
427 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*factor
428 END DO
429
430 w = (zeta1*p + eta*q)*zetapetainv
431 product_list(nproducts)%ra = ra
432 product_list(nproducts)%rb = rb
433 product_list(nproducts)%rc = c11
434 product_list(nproducts)%rd = tmp_d
435 product_list(nproducts)%ZetapEtaInv = zetapetainv
436 product_list(nproducts)%Rho = rho
437 product_list(nproducts)%RhoInv = rhoinv
438 product_list(nproducts)%P = p
439 product_list(nproducts)%Q = q
440 product_list(nproducts)%W = w
441 product_list(nproducts)%AB = ra - rb
442 product_list(nproducts)%CD = c11 - tmp_d
443 END DO
444 END DO
445 END DO
446
447 END SUBROUTINE build_pgf_product_list
448
449! **************************************************************************************************
450!> \brief ...
451!> \param npgfa ...
452!> \param npgfb ...
453!> \param list ...
454!> \param zeta ...
455!> \param zetb ...
456!> \param screen1 ...
457!> \param screen2 ...
458!> \param pgf ...
459!> \param R_pgf ...
460!> \param log10_pmax ...
461!> \param log10_eps_schwarz ...
462!> \param ra ...
463!> \param rb ...
464!> \param nelements ...
465!> \param neighbor_cells ...
466!> \param nimages ...
467!> \param do_periodic ...
468! **************************************************************************************************
469 SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, &
470 log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
471 neighbor_cells, nimages, do_periodic)
472 INTEGER, INTENT(IN) :: npgfa, npgfb
473 TYPE(hfx_pgf_list), DIMENSION(npgfa*npgfb) :: list
474 REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
475 REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
476 REAL(dp), INTENT(IN) :: screen1(2), screen2(2)
477 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
478 POINTER :: pgf, r_pgf
479 REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), &
480 rb(3)
481 INTEGER, INTENT(OUT) :: nelements
482 TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
483 INTEGER :: nimages(npgfa*npgfb)
484 LOGICAL, INTENT(IN) :: do_periodic
485
486 INTEGER :: element_counter, i, ipgf, j, jpgf
487 REAL(dp) :: ab(3), im_b(3), pgf_max, rab2, zeta1, &
488 zeta_a, zeta_b, zetainv
489
490 nimages = 0
491 ! ** inner loop may never be reached
492 nelements = npgfa*npgfb
493 DO i = 1, SIZE(neighbor_cells)
494 IF (do_periodic) THEN
495 im_b = rb + neighbor_cells(i)%cell_r(:)
496 ELSE
497 im_b = rb
498 END IF
499 ab = ra - im_b
500 rab2 = ab(1)**2 + ab(2)**2 + ab(3)**2
501 IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) cycle
502 element_counter = 0
503 DO ipgf = 1, npgfa
504 DO jpgf = 1, npgfb
505 element_counter = element_counter + 1
506 pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
507 IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
508 cycle
509 END IF
510 nimages(element_counter) = nimages(element_counter) + 1
511 list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
512 list(element_counter)%image_list(nimages(element_counter))%ra = ra
513 list(element_counter)%image_list(nimages(element_counter))%rb = im_b
514 list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
515
516 zeta_a = zeta(ipgf)
517 zeta_b = zetb(jpgf)
518 zeta1 = zeta_a + zeta_b
519 zetainv = 1.0_dp/zeta1
520
521 IF (nimages(element_counter) == 1) THEN
522 list(element_counter)%ipgf = ipgf
523 list(element_counter)%jpgf = jpgf
524 list(element_counter)%zetaInv = zetainv
525 list(element_counter)%zetapzetb = zeta1
526 list(element_counter)%zeta = zeta_a
527 list(element_counter)%zetb = zeta_b
528 END IF
529
530 list(element_counter)%image_list(nimages(element_counter))%S1234 = (-zeta_a*zeta_b*zetainv*rab2)
531 list(element_counter)%image_list(nimages(element_counter))%P = (zeta_a*ra + zeta_b*im_b)*zetainv
532 list(element_counter)%image_list(nimages(element_counter))%R = &
533 max(0.0_dp, r_pgf(jpgf, ipgf)%x(1)*rab2 + r_pgf(jpgf, ipgf)%x(2))
534 list(element_counter)%image_list(nimages(element_counter))%ra = ra
535 list(element_counter)%image_list(nimages(element_counter))%rb = im_b
536 list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
537 list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
538 END DO
539 END DO
540 nelements = max(nelements, element_counter)
541 END DO
542 DO j = 1, nelements
543 list(j)%nimages = nimages(j)
544 END DO
545 ! ** Remove unused elements
546
547 element_counter = 0
548 DO j = 1, nelements
549 IF (list(j)%nimages == 0) cycle
550 element_counter = element_counter + 1
551 list(element_counter)%nimages = list(j)%nimages
552 list(element_counter)%zetapzetb = list(j)%zetapzetb
553 list(element_counter)%ZetaInv = list(j)%ZetaInv
554 list(element_counter)%zeta = list(j)%zeta
555 list(element_counter)%zetb = list(j)%zetb
556 list(element_counter)%ipgf = list(j)%ipgf
557 list(element_counter)%jpgf = list(j)%jpgf
558 DO i = 1, list(j)%nimages
559 list(element_counter)%image_list(i) = list(j)%image_list(i)
560 END DO
561 END DO
562
563 nelements = element_counter
564
565 END SUBROUTINE build_pair_list_pgf
566
567! **************************************************************************************************
568!> \brief ...
569!> \param natom ...
570!> \param list ...
571!> \param set_list ...
572!> \param i_start ...
573!> \param i_end ...
574!> \param j_start ...
575!> \param j_end ...
576!> \param kind_of ...
577!> \param basis_parameter ...
578!> \param particle_set ...
579!> \param do_periodic ...
580!> \param coeffs_set ...
581!> \param coeffs_kind ...
582!> \param coeffs_kind_max0 ...
583!> \param log10_eps_schwarz ...
584!> \param cell ...
585!> \param pmax_blocks ...
586!> \param atomic_pair_list ...
587! **************************************************************************************************
588 SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
589 do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
590 pmax_blocks, atomic_pair_list)
591
592 INTEGER, INTENT(IN) :: natom
593 TYPE(pair_list_type), INTENT(INOUT) :: list
594 TYPE(pair_set_list_type), DIMENSION(:), &
595 INTENT(OUT) :: set_list
596 INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
597 kind_of(*)
598 TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
599 POINTER :: basis_parameter
600 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
601 POINTER :: particle_set
602 LOGICAL, INTENT(IN) :: do_periodic
603 TYPE(hfx_screen_coeff_type), &
604 DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
605 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
606 INTENT(IN) :: coeffs_kind
607 REAL(kind=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
608 TYPE(cell_type), POINTER :: cell
609 REAL(dp), INTENT(IN) :: pmax_blocks
610 LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
611
612 INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
613 n_element, nset_ij, nseta, nsetb
614 REAL(kind=dp) :: rab2
615 REAL(kind=dp), DIMENSION(3) :: b11, pbc_b, ra, rb, temp
616
617 n_element = 0
618 nset_ij = 0
619
620 DO iatom = i_start, i_end
621 DO jatom = j_start, j_end
622 IF (atomic_pair_list(jatom, iatom) .EQV. .false.) cycle
623
624 ikind = kind_of(iatom)
625 nseta = basis_parameter(ikind)%nset
626 ra = particle_set(iatom)%r(:)
627
628 IF (jatom < iatom) cycle
629 jkind = kind_of(jatom)
630 nsetb = basis_parameter(jkind)%nset
631 rb = particle_set(jatom)%r(:)
632
633 IF (do_periodic) THEN
634 temp = rb - ra
635 pbc_b = pbc(temp, cell)
636 b11 = ra + pbc_b
637 rab2 = (ra(1) - b11(1))**2 + (ra(2) - b11(2))**2 + (ra(3) - b11(3))**2
638 ELSE
639 rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
640 b11 = rb ! ra - rb
641 END IF
642 IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
643 coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
644
645 n_element = n_element + 1
646 list%elements(n_element)%pair = [iatom, jatom]
647 list%elements(n_element)%kind_pair = [ikind, jkind]
648 list%elements(n_element)%r1 = ra
649 list%elements(n_element)%r2 = b11
650 list%elements(n_element)%dist2 = rab2
651 ! build a list of guaranteed overlapping sets
652 list%elements(n_element)%set_bounds(1) = nset_ij + 1
653 DO iset = 1, nseta
654 DO jset = 1, nsetb
655 IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
656 coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
657 nset_ij = nset_ij + 1
658 set_list(nset_ij)%pair = [iset, jset]
659 END DO
660 END DO
661 list%elements(n_element)%set_bounds(2) = nset_ij
662 END DO
663 END DO
664
665 list%n_element = n_element
666
667 END SUBROUTINE build_pair_list
668
669! **************************************************************************************************
670!> \brief ...
671!> \param natom ...
672!> \param atomic_pair_list ...
673!> \param kind_of ...
674!> \param basis_parameter ...
675!> \param particle_set ...
676!> \param do_periodic ...
677!> \param coeffs_kind ...
678!> \param coeffs_kind_max0 ...
679!> \param log10_eps_schwarz ...
680!> \param cell ...
681!> \param blocks ...
682! **************************************************************************************************
683 SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
684 do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
685 blocks)
686 INTEGER, INTENT(IN) :: natom
687 LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list
688 INTEGER, INTENT(IN) :: kind_of(*)
689 TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
690 POINTER :: basis_parameter
691 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
692 POINTER :: particle_set
693 LOGICAL, INTENT(IN) :: do_periodic
694 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
695 INTENT(IN) :: coeffs_kind
696 REAL(kind=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
697 TYPE(cell_type), POINTER :: cell
698 TYPE(hfx_block_range_type), DIMENSION(:), &
699 INTENT(IN), POINTER :: blocks
700
701 INTEGER :: iatom, iatom_end, iatom_start, iblock, &
702 ikind, jatom, jatom_end, jatom_start, &
703 jblock, jkind, nseta, nsetb
704 REAL(kind=dp) :: rab2
705 REAL(kind=dp), DIMENSION(3) :: b11, pbc_b, ra, rb, temp
706
707 atomic_pair_list = .false.
708
709 DO iblock = 1, SIZE(blocks)
710 iatom_start = blocks(iblock)%istart
711 iatom_end = blocks(iblock)%iend
712 DO jblock = 1, SIZE(blocks)
713 jatom_start = blocks(jblock)%istart
714 jatom_end = blocks(jblock)%iend
715
716 DO iatom = iatom_start, iatom_end
717 ikind = kind_of(iatom)
718 nseta = basis_parameter(ikind)%nset
719 ra = particle_set(iatom)%r(:)
720 DO jatom = jatom_start, jatom_end
721 IF (jatom < iatom) cycle
722 jkind = kind_of(jatom)
723 nsetb = basis_parameter(jkind)%nset
724 rb = particle_set(jatom)%r(:)
725
726 IF (do_periodic) THEN
727 temp = rb - ra
728 pbc_b = pbc(temp, cell)
729 b11 = ra + pbc_b
730 rab2 = (ra(1) - b11(1))**2 + (ra(2) - b11(2))**2 + (ra(3) - b11(3))**2
731 ELSE
732 rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
733 b11 = rb ! ra - rb
734 END IF
735 IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
736 coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) cycle
737
738 atomic_pair_list(jatom, iatom) = .true.
739 atomic_pair_list(iatom, jatom) = .true.
740 END DO
741 END DO
742 END DO
743 END DO
744
745 END SUBROUTINE build_atomic_pair_list
746
747! **************************************************************************************************
748!> \brief ...
749!> \param natom ...
750!> \param list ...
751!> \param set_list ...
752!> \param i_start ...
753!> \param i_end ...
754!> \param j_start ...
755!> \param j_end ...
756!> \param kind_of ...
757!> \param basis_parameter ...
758!> \param particle_set ...
759!> \param do_periodic ...
760!> \param coeffs_set ...
761!> \param coeffs_kind ...
762!> \param coeffs_kind_max0 ...
763!> \param log10_eps_schwarz ...
764!> \param cell ...
765!> \param pmax_blocks ...
766!> \param atomic_pair_list ...
767!> \param skip_atom_symmetry ...
768! **************************************************************************************************
769 SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
770 do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
771 pmax_blocks, atomic_pair_list, skip_atom_symmetry)
772
773 INTEGER, INTENT(IN) :: natom
774 TYPE(pair_list_type_mp2), INTENT(INOUT) :: list
775 TYPE(pair_set_list_type), DIMENSION(:), &
776 INTENT(OUT) :: set_list
777 INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
778 kind_of(*)
779 TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
780 POINTER :: basis_parameter
781 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
782 POINTER :: particle_set
783 LOGICAL, INTENT(IN) :: do_periodic
784 TYPE(hfx_screen_coeff_type), &
785 DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
786 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
787 INTENT(IN) :: coeffs_kind
788 REAL(kind=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
789 TYPE(cell_type), POINTER :: cell
790 REAL(dp), INTENT(IN) :: pmax_blocks
791 LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
792 LOGICAL, INTENT(IN), OPTIONAL :: skip_atom_symmetry
793
794 INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
795 n_element, nset_ij, nseta, nsetb
796 LOGICAL :: my_skip_atom_symmetry
797 REAL(kind=dp) :: rab2
798 REAL(kind=dp), DIMENSION(3) :: b11, pbc_b, ra, rb, temp
799
800 n_element = 0
801 nset_ij = 0
802
803 my_skip_atom_symmetry = .false.
804 IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry
805
806 DO iatom = i_start, i_end
807 DO jatom = j_start, j_end
808 IF (atomic_pair_list(jatom, iatom) .EQV. .false.) cycle
809
810 ikind = kind_of(iatom)
811 nseta = basis_parameter(ikind)%nset
812 ra = particle_set(iatom)%r(:)
813
814 IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) cycle
815 jkind = kind_of(jatom)
816 nsetb = basis_parameter(jkind)%nset
817 rb = particle_set(jatom)%r(:)
818
819 IF (do_periodic) THEN
820 temp = rb - ra
821 pbc_b = pbc(temp, cell)
822 b11 = ra + pbc_b
823 rab2 = (ra(1) - b11(1))**2 + (ra(2) - b11(2))**2 + (ra(3) - b11(3))**2
824 ELSE
825 rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
826 b11 = rb ! ra - rb
827 END IF
828 IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
829 coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
830
831 n_element = n_element + 1
832 list%elements(n_element)%pair = [iatom, jatom]
833 list%elements(n_element)%kind_pair = [ikind, jkind]
834 list%elements(n_element)%r1 = ra
835 list%elements(n_element)%r2 = b11
836 list%elements(n_element)%dist2 = rab2
837 ! build a list of guaranteed overlapping sets
838 list%elements(n_element)%set_bounds(1) = nset_ij + 1
839 DO iset = 1, nseta
840 DO jset = 1, nsetb
841 IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
842 coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
843 nset_ij = nset_ij + 1
844 set_list(nset_ij)%pair = [iset, jset]
845 END DO
846 END DO
847 list%elements(n_element)%set_bounds(2) = nset_ij
848 END DO
849 END DO
850
851 list%n_element = n_element
852
853 END SUBROUTINE build_pair_list_mp2
854
855END MODULE hfx_pair_list_methods
cell_types::pbc
Definition cell_types.F:103
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
hfx_pair_list_methods
Routines for optimizing load balance between processes in HFX calculations.
Definition hfx_pair_list_methods.F:15
hfx_pair_list_methods::build_pgf_product_list
subroutine, public build_pgf_product_list(list1, list2, product_list, nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, potential_parameter, m_max, do_periodic)
...
Definition hfx_pair_list_methods.F:196
hfx_pair_list_methods::build_atomic_pair_list
subroutine, public build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, blocks)
...
Definition hfx_pair_list_methods.F:686
hfx_pair_list_methods::build_pair_list
subroutine, public build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, pmax_blocks, atomic_pair_list)
...
Definition hfx_pair_list_methods.F:591
hfx_pair_list_methods::build_pair_list_pgf
subroutine, public build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, r_pgf, log10_pmax, log10_eps_schwarz, ra, rb, nelements, neighbor_cells, nimages, do_periodic)
...
Definition hfx_pair_list_methods.F:472
hfx_pair_list_methods::build_pair_list_mp2
subroutine, public build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, pmax_blocks, atomic_pair_list, skip_atom_symmetry)
...
Definition hfx_pair_list_methods.F:772
hfx_pair_list_methods::allocate_bra
subroutine, public allocate_bra(bra, max_npgf, max_ncell)
Allocates all arrays within a Bra_t structure.
Definition hfx_pair_list_methods.F:77
hfx_pair_list_methods::build_pair_list_pbc_pgf
subroutine, public build_pair_list_pbc_pgf(npgfa, npgfb, bra, screen1, screen2, pgf, log10_pmax, log10_eps_schwarz, ra, rb, nelements, neighbor_cells, do_periodic)
Builds a screened list of primitives from centers A and B, intersecting with another shell.
Definition hfx_pair_list_methods.F:118
hfx_pair_list_methods::pgf_product_list_size
integer, save, public pgf_product_list_size
Definition hfx_pair_list_methods.F:48
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:16
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_potential_mix_cl
integer, parameter, public do_potential_mix_cl
Definition input_constants.F:834
input_constants::do_potential_gaussian
integer, parameter, public do_potential_gaussian
Definition input_constants.F:834
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:834
input_constants::do_potential_mix_lg
integer, parameter, public do_potential_mix_lg
Definition input_constants.F:834
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:834
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:834
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:834
input_constants::do_potential_mix_cl_trunc
integer, parameter, public do_potential_mix_cl_trunc
Definition input_constants.F:834
input_constants::do_potential_long
integer, parameter, public do_potential_long
Definition input_constants.F:834
input_constants::do_potential_tshpsc
integer, parameter, public do_potential_tshpsc
Definition input_constants.F:834
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
libint_wrapper::prim_data_f_size
integer, parameter, public prim_data_f_size
Definition libint_wrapper.F:65
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:58
t_c_g0::t_c_g0_n
subroutine, public t_c_g0_n(res, use_gamma, r, t, nderiv)
...
Definition t_c_g0.F:88
t_sh_p_s_c
...
Definition t_sh_p_s_c.F:39
t_sh_p_s_c::trunc_cs_poly_n20
subroutine, public trunc_cs_poly_n20(res, r, t, nderiv)
...
Definition t_sh_p_s_c.F:65
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
hfx_pair_list_methods::bra_t
Definition hfx_pair_list_methods.F:51
hfx_types::hfx_basis_type
Definition hfx_types.F:267
hfx_types::hfx_block_range_type
Definition hfx_types.F:344
hfx_types::hfx_cell_type
Definition hfx_types.F:208
hfx_types::hfx_pgf_list
Definition hfx_types.F:323
hfx_types::hfx_pgf_product_list
Definition hfx_types.F:334
hfx_types::hfx_potential_type
Definition hfx_types.F:141
hfx_types::hfx_screen_coeff_type
Definition hfx_types.F:296
hfx_types::pair_list_type
Definition hfx_types.F:239
hfx_types::pair_set_list_type
Definition hfx_types.F:234
mp2_types::pair_list_type_mp2
Definition mp2_types.F:423
particle_types::particle_type
Definition particle_types.F:35