Line data Source code
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 : ! **************************************************************************************************
15 : MODULE hfx_pair_list_methods
16 : USE cell_types, ONLY: cell_type,&
17 : pbc
18 : USE four_center_integral_methods, ONLY: min_image_gap_sq
19 : USE hfx_types, ONLY: hfx_basis_type,&
20 : hfx_block_range_type,&
21 : hfx_cell_type,&
22 : hfx_screen_coeff_type,&
23 : pair_list_type,&
24 : pair_set_list_type
25 : USE kinds, ONLY: dp
26 : USE mp2_types, ONLY: pair_list_type_mp2
27 : USE particle_types, ONLY: particle_type
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 : PRIVATE
32 :
33 : PUBLIC :: build_pair_list, &
34 : build_pair_list_mp2, &
35 : build_atomic_pair_list, &
36 : build_pair_list_pbc_pgf, &
37 : bra_t, allocate_bra, min_image_gap_sq
38 :
39 : ! Store screened primitive gaussian function pairs between two shells (A and B)
40 : TYPE :: Bra_t
41 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pgf_scr
42 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: pgf_idx, cell_idx
43 : INTEGER :: cell_cnt = 0, pgf_cnt = 0
44 : END TYPE Bra_t
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Allocates all arrays within a Bra_t structure.
50 : !>
51 : !> This subroutine ensures that a given `Bra_t` pointer has allocated memory
52 : !> for storing primitive Gaussian pair indices, screening data, and cell indices.
53 : !> Existing allocations are safely deallocated first.
54 : !>
55 : !> \param[in,out] bra Pointer to the Bra_t structure to allocate.
56 : !> \param[in] max_npgf Maximum number of primitives per shell.
57 : !> \param[in] max_ncell Maximum number of neighbor cells.
58 : !>
59 : !> \note The arrays are allocated as follows:
60 : !> - `pgf_idx(3, max_npgf^2 * max_ncell)`
61 : !> - `pgf_scr(5, max_npgf^2 * max_ncell)`
62 : !> - `cell_idx(3, max_ncell)`
63 : !>
64 : ! **************************************************************************************************
65 0 : SUBROUTINE allocate_Bra(bra, max_npgf, max_ncell)
66 : TYPE(bra_t), POINTER :: bra
67 : INTEGER :: max_npgf, max_ncell
68 :
69 0 : IF (ALLOCATED(bra%pgf_idx)) DEALLOCATE (bra%pgf_idx) ! TODO add if n_prm changed
70 0 : IF (ALLOCATED(bra%pgf_scr)) DEALLOCATE (bra%pgf_scr) ! TODO add if n_prm changed
71 0 : IF (ALLOCATED(bra%cell_idx)) DEALLOCATE (bra%cell_idx) ! TODO add if n_prm changed
72 0 : ALLOCATE (bra%pgf_idx(3, max_npgf*max_npgf*max_ncell))
73 0 : ALLOCATE (bra%pgf_scr(5, max_npgf*max_npgf*max_ncell))
74 0 : ALLOCATE (bra%cell_idx(3, max_ncell))
75 0 : END SUBROUTINE allocate_Bra
76 :
77 : ! **************************************************************************************************
78 : !> \brief Builds a screened list of primitives from centers A and B, intersecting with another shell
79 : !>
80 : !> This subroutine populates the `Bra_t` structure with indices and screening data
81 : !> for primitive Gaussian pairs (ipgf,jpgf,n3) that pass the Schwarz screening criterion
82 : !>
83 : !> \param[in] npgfa Number of primitives in shell A
84 : !> \param[in] npgfb Number of primitives in shell B
85 : !> \param[out] bra Pointer to Bra_t structure to populate
86 : !> \param[in] screen1 Screening coeffiecients for the AB pair based on most diffuse pgf
87 : !> \param[in] screen2 Screening coeffiecients for the CD pair based on most diffuse pgf
88 : !> \param[in] pgf Screening coeeficents for primitive gaussians
89 : !> \param[in] log10_pmax Log10 of maximum integral prefactor
90 : !> \param[in] log10_eps_schwarz Log10 of Schwarz screening threshold
91 : !> \param[in] ra Cartesian coordinates of center A
92 : !> \param[in] rb Cartesian coordinates of center B
93 : !> \param[out] nelements Total number of valid primitives found
94 : !> \param[in] neighbor_cells Array of lattice vectors
95 : !> \param[in] do_periodic Logical flag for periodic boundary conditions on B
96 : !>
97 : !> \note Loops over all cells, shifting the B shell position accordingly
98 : !> Only pairs satisfying the screening condition are stored.
99 : !>
100 : !> \note Each valid cell contributes a block of pairs stored contiguously in
101 : !> bra%pgf_idx, with cell metadata stored in bra%cell_idx
102 : !> \note screening primitive pairs (ipgf, jpgf) using both coarse and fine screening thresholds
103 : ! **************************************************************************************************
104 0 : SUBROUTINE build_pair_list_pbc_pgf(npgfa, npgfb, bra, screen1, screen2, pgf, &
105 : log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
106 : neighbor_cells, do_periodic)
107 :
108 : INTEGER, INTENT(IN) :: npgfa, npgfb
109 : TYPE(bra_t), POINTER :: bra
110 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2)
111 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
112 : POINTER :: pgf
113 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), &
114 : rb(3)
115 : INTEGER, INTENT(OUT) :: nelements
116 : TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
117 : LOGICAL, INTENT(IN) :: do_periodic
118 :
119 : INTEGER :: cell_cnt, cell_idx, element_cnt, &
120 : element_idx, element_off, ipgf, jpgf
121 : REAL(dp) :: AB(3), im_B(3), pgf_max, rab2
122 :
123 0 : element_idx = 0
124 0 : element_off = 0
125 0 : cell_cnt = 0
126 0 : DO cell_idx = 1, SIZE(neighbor_cells)
127 : ! move B to this cell while keeping A fixed
128 : ! NOTE rb has already been moved by build_pair_list
129 : ! so that when cell_r is (000) AB is in the pbc cell
130 0 : IF (do_periodic) THEN
131 0 : im_B = rb + neighbor_cells(cell_idx)%cell_r(:)
132 : ELSE
133 0 : im_B = rb
134 : END IF
135 0 : AB = ra - im_B
136 0 : rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
137 : ! First screening on the most diffuse AB gaussians along with the most diffuse CD gaussian
138 0 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
139 :
140 0 : element_cnt = 0
141 0 : DO ipgf = 1, npgfa
142 0 : DO jpgf = 1, npgfb
143 :
144 : ! Second screening on this actual AB pair and the most diffuse CD from before
145 0 : pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
146 0 : IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
147 : ! This pair passed screening. Add it to the bra.
148 0 : element_idx = element_idx + 1
149 0 : bra%pgf_idx(:, element_idx) = [ipgf, jpgf, cell_idx]
150 0 : bra%pgf_scr(1, element_idx) = pgf_max
151 0 : element_cnt = element_cnt + 1
152 : END DO
153 : END DO
154 :
155 : ! If this cell produced any pair, add this cell to the bra
156 0 : IF (element_cnt == 0) CYCLE
157 0 : cell_cnt = cell_cnt + 1
158 0 : bra%cell_idx(:, cell_cnt) = [cell_idx, element_cnt, element_off]
159 0 : element_off = element_off + element_cnt
160 :
161 : END DO ! cell
162 0 : bra%cell_cnt = cell_cnt
163 0 : nelements = element_idx
164 0 : bra%pgf_cnt = nelements
165 :
166 0 : END SUBROUTINE build_pair_list_pbc_pgf
167 :
168 : ! **************************************************************************************************
169 : !> \brief ...
170 : !> \param natom ...
171 : !> \param list ...
172 : !> \param set_list ...
173 : !> \param i_start ...
174 : !> \param i_end ...
175 : !> \param j_start ...
176 : !> \param j_end ...
177 : !> \param kind_of ...
178 : !> \param basis_parameter ...
179 : !> \param particle_set ...
180 : !> \param do_periodic ...
181 : !> \param coeffs_set ...
182 : !> \param coeffs_kind ...
183 : !> \param coeffs_kind_max0 ...
184 : !> \param log10_eps_schwarz ...
185 : !> \param cell ...
186 : !> \param pmax_blocks ...
187 : !> \param atomic_pair_list ...
188 : ! **************************************************************************************************
189 52847672 : SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
190 3822772 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
191 3822772 : pmax_blocks, atomic_pair_list)
192 :
193 : INTEGER, INTENT(IN) :: natom
194 : TYPE(pair_list_type), INTENT(INOUT) :: list
195 : TYPE(pair_set_list_type), DIMENSION(:), &
196 : INTENT(OUT) :: set_list
197 : INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
198 : kind_of(*)
199 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
200 : POINTER :: basis_parameter
201 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
202 : POINTER :: particle_set
203 : LOGICAL, INTENT(IN) :: do_periodic
204 : TYPE(hfx_screen_coeff_type), &
205 : DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
206 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
207 : INTENT(IN) :: coeffs_kind
208 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
209 : TYPE(cell_type), POINTER :: cell
210 : REAL(dp), INTENT(IN) :: pmax_blocks
211 : LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
212 :
213 : INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
214 : n_element, nset_ij, nseta, nsetb
215 : REAL(KIND=dp) :: rab2
216 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
217 :
218 3822772 : n_element = 0
219 3822772 : nset_ij = 0
220 :
221 7645544 : DO iatom = i_start, i_end
222 11468316 : DO jatom = j_start, j_end
223 3822772 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
224 :
225 3689122 : ikind = kind_of(iatom)
226 3689122 : nseta = basis_parameter(ikind)%nset
227 14756488 : ra = particle_set(iatom)%r(:)
228 :
229 3689122 : IF (jatom < iatom) CYCLE
230 3689122 : jkind = kind_of(jatom)
231 3689122 : nsetb = basis_parameter(jkind)%nset
232 14756488 : rb = particle_set(jatom)%r(:)
233 :
234 3689122 : IF (do_periodic) THEN
235 4695160 : temp = rb - ra
236 1173790 : pbc_B = pbc(temp, cell)
237 4695160 : B11 = ra + pbc_B
238 1173790 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
239 : ELSE
240 2515332 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
241 2515332 : B11 = rb ! ra - rb
242 : END IF
243 3689122 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
244 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
245 :
246 3682634 : n_element = n_element + 1
247 11047902 : list%elements(n_element)%pair = [iatom, jatom]
248 11047902 : list%elements(n_element)%kind_pair = [ikind, jkind]
249 14730536 : list%elements(n_element)%r1 = ra
250 14730536 : list%elements(n_element)%r2 = B11
251 3682634 : list%elements(n_element)%dist2 = rab2
252 : ! build a list of guaranteed overlapping sets
253 3682634 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
254 12429552 : DO iset = 1, nseta
255 36933438 : DO jset = 1, nsetb
256 24503886 : IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
257 : coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
258 23300342 : nset_ij = nset_ij + 1
259 79851488 : set_list(nset_ij)%pair = [iset, jset]
260 : END DO
261 : END DO
262 7645544 : list%elements(n_element)%set_bounds(2) = nset_ij
263 : END DO
264 : END DO
265 :
266 3822772 : list%n_element = n_element
267 :
268 3822772 : END SUBROUTINE build_pair_list
269 :
270 : ! **************************************************************************************************
271 : !> \brief ...
272 : !> \param natom ...
273 : !> \param atomic_pair_list ...
274 : !> \param kind_of ...
275 : !> \param basis_parameter ...
276 : !> \param particle_set ...
277 : !> \param do_periodic ...
278 : !> \param coeffs_kind ...
279 : !> \param coeffs_kind_max0 ...
280 : !> \param log10_eps_schwarz ...
281 : !> \param cell ...
282 : !> \param blocks ...
283 : ! **************************************************************************************************
284 6238 : SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
285 6238 : do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
286 : blocks)
287 : INTEGER, INTENT(IN) :: natom
288 : LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list
289 : INTEGER, INTENT(IN) :: kind_of(*)
290 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
291 : POINTER :: basis_parameter
292 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
293 : POINTER :: particle_set
294 : LOGICAL, INTENT(IN) :: do_periodic
295 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
296 : INTENT(IN) :: coeffs_kind
297 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
298 : TYPE(cell_type), POINTER :: cell
299 : TYPE(hfx_block_range_type), DIMENSION(:), &
300 : INTENT(IN), POINTER :: blocks
301 :
302 : INTEGER :: iatom, iatom_end, iatom_start, iblock, &
303 : ikind, jatom, jatom_end, jatom_start, &
304 : jblock, jkind, nseta, nsetb
305 : REAL(KIND=dp) :: rab2
306 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
307 :
308 93586 : atomic_pair_list = .FALSE.
309 :
310 25096 : DO iblock = 1, SIZE(blocks)
311 18858 : iatom_start = blocks(iblock)%istart
312 18858 : iatom_end = blocks(iblock)%iend
313 93586 : DO jblock = 1, SIZE(blocks)
314 68490 : jatom_start = blocks(jblock)%istart
315 68490 : jatom_end = blocks(jblock)%iend
316 :
317 155838 : DO iatom = iatom_start, iatom_end
318 68490 : ikind = kind_of(iatom)
319 68490 : nseta = basis_parameter(ikind)%nset
320 273960 : ra = particle_set(iatom)%r(:)
321 205470 : DO jatom = jatom_start, jatom_end
322 68490 : IF (jatom < iatom) CYCLE
323 43674 : jkind = kind_of(jatom)
324 43674 : nsetb = basis_parameter(jkind)%nset
325 174696 : rb = particle_set(jatom)%r(:)
326 :
327 43674 : IF (do_periodic) THEN
328 50504 : temp = rb - ra
329 12626 : pbc_B = pbc(temp, cell)
330 50504 : B11 = ra + pbc_B
331 12626 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
332 : ELSE
333 31048 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
334 31048 : B11 = rb ! ra - rb
335 : END IF
336 43674 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
337 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) CYCLE
338 :
339 43218 : atomic_pair_list(jatom, iatom) = .TRUE.
340 136980 : atomic_pair_list(iatom, jatom) = .TRUE.
341 : END DO
342 : END DO
343 : END DO
344 : END DO
345 :
346 6238 : END SUBROUTINE build_atomic_pair_list
347 :
348 : ! **************************************************************************************************
349 : !> \brief ...
350 : !> \param natom ...
351 : !> \param list ...
352 : !> \param set_list ...
353 : !> \param i_start ...
354 : !> \param i_end ...
355 : !> \param j_start ...
356 : !> \param j_end ...
357 : !> \param kind_of ...
358 : !> \param basis_parameter ...
359 : !> \param particle_set ...
360 : !> \param do_periodic ...
361 : !> \param coeffs_set ...
362 : !> \param coeffs_kind ...
363 : !> \param coeffs_kind_max0 ...
364 : !> \param log10_eps_schwarz ...
365 : !> \param cell ...
366 : !> \param pmax_blocks ...
367 : !> \param atomic_pair_list ...
368 : !> \param skip_atom_symmetry ...
369 : ! **************************************************************************************************
370 595568 : SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
371 2994 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
372 2994 : pmax_blocks, atomic_pair_list, skip_atom_symmetry)
373 :
374 : INTEGER, INTENT(IN) :: natom
375 : TYPE(pair_list_type_mp2), INTENT(INOUT) :: list
376 : TYPE(pair_set_list_type), DIMENSION(:), &
377 : INTENT(OUT) :: set_list
378 : INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
379 : kind_of(*)
380 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
381 : POINTER :: basis_parameter
382 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
383 : POINTER :: particle_set
384 : LOGICAL, INTENT(IN) :: do_periodic
385 : TYPE(hfx_screen_coeff_type), &
386 : DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
387 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
388 : INTENT(IN) :: coeffs_kind
389 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
390 : TYPE(cell_type), POINTER :: cell
391 : REAL(dp), INTENT(IN) :: pmax_blocks
392 : LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
393 : LOGICAL, INTENT(IN), OPTIONAL :: skip_atom_symmetry
394 :
395 : INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
396 : n_element, nset_ij, nseta, nsetb
397 : LOGICAL :: my_skip_atom_symmetry
398 : REAL(KIND=dp) :: rab2
399 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
400 :
401 2994 : n_element = 0
402 2994 : nset_ij = 0
403 :
404 2994 : my_skip_atom_symmetry = .FALSE.
405 2994 : IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry
406 :
407 6076 : DO iatom = i_start, i_end
408 9454 : DO jatom = j_start, j_end
409 3378 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
410 :
411 3378 : ikind = kind_of(iatom)
412 3378 : nseta = basis_parameter(ikind)%nset
413 13512 : ra = particle_set(iatom)%r(:)
414 :
415 3378 : IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) CYCLE
416 3230 : jkind = kind_of(jatom)
417 3230 : nsetb = basis_parameter(jkind)%nset
418 12920 : rb = particle_set(jatom)%r(:)
419 :
420 3230 : IF (do_periodic) THEN
421 0 : temp = rb - ra
422 0 : pbc_B = pbc(temp, cell)
423 0 : B11 = ra + pbc_B
424 0 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
425 : ELSE
426 3230 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
427 3230 : B11 = rb ! ra - rb
428 : END IF
429 3230 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
430 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
431 :
432 3230 : n_element = n_element + 1
433 9690 : list%elements(n_element)%pair = [iatom, jatom]
434 9690 : list%elements(n_element)%kind_pair = [ikind, jkind]
435 12920 : list%elements(n_element)%r1 = ra
436 12920 : list%elements(n_element)%r2 = B11
437 3230 : list%elements(n_element)%dist2 = rab2
438 : ! build a list of guaranteed overlapping sets
439 3230 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
440 13988 : DO iset = 1, nseta
441 52586 : DO jset = 1, nsetb
442 38598 : IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
443 : coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
444 38526 : nset_ij = nset_ij + 1
445 126408 : set_list(nset_ij)%pair = [iset, jset]
446 : END DO
447 : END DO
448 6460 : list%elements(n_element)%set_bounds(2) = nset_ij
449 : END DO
450 : END DO
451 :
452 2994 : list%n_element = n_element
453 :
454 2994 : END SUBROUTINE build_pair_list_mp2
455 :
456 0 : END MODULE hfx_pair_list_methods
|
