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 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 :
57 : CONTAINS
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 0 : SUBROUTINE allocate_Bra(bra, max_npgf, max_ncell)
77 : TYPE(bra_t), POINTER :: bra
78 : INTEGER :: max_npgf, max_ncell
79 :
80 0 : IF (ALLOCATED(bra%pgf_idx)) DEALLOCATE (bra%pgf_idx) ! TODO add if n_prm changed
81 0 : IF (ALLOCATED(bra%pgf_scr)) DEALLOCATE (bra%pgf_scr) ! TODO add if n_prm changed
82 0 : IF (ALLOCATED(bra%cell_idx)) DEALLOCATE (bra%cell_idx) ! TODO add if n_prm changed
83 0 : ALLOCATE (bra%pgf_idx(3, max_npgf*max_npgf*max_ncell))
84 0 : ALLOCATE (bra%pgf_scr(5, max_npgf*max_npgf*max_ncell))
85 0 : ALLOCATE (bra%cell_idx(3, max_ncell))
86 0 : 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 0 : 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 0 : element_idx = 0
135 0 : element_off = 0
136 0 : cell_cnt = 0
137 0 : 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 0 : IF (do_periodic) THEN
142 0 : im_B = rb + neighbor_cells(cell_idx)%cell_r(:)
143 : ELSE
144 0 : im_B = rb
145 : END IF
146 0 : AB = ra - im_B
147 0 : 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 0 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
150 :
151 0 : element_cnt = 0
152 0 : DO ipgf = 1, npgfa
153 0 : DO jpgf = 1, npgfb
154 :
155 : ! Second screening on this actual AB pair and the most diffuse CD from before
156 0 : pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
157 0 : IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
158 : ! This pair passed screening. Add it to the bra.
159 0 : element_idx = element_idx + 1
160 0 : bra%pgf_idx(:, element_idx) = [ipgf, jpgf, cell_idx]
161 0 : bra%pgf_scr(1, element_idx) = pgf_max
162 0 : 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 0 : IF (element_cnt == 0) CYCLE
168 0 : cell_cnt = cell_cnt + 1
169 0 : bra%cell_idx(:, cell_cnt) = [cell_idx, element_cnt, element_off]
170 0 : element_off = element_off + element_cnt
171 :
172 : END DO ! cell
173 0 : bra%cell_cnt = cell_cnt
174 0 : nelements = element_idx
175 0 : bra%pgf_cnt = nelements
176 :
177 0 : 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 90774491 : 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 90774491 : DIMENSION(:) :: tmp_product_list
216 :
217 90774491 : nimages1 = list1%nimages
218 90774491 : nimages2 = list2%nimages
219 90774491 : nproducts = 0
220 90774491 : Zeta1 = list1%zetapzetb
221 90774491 : Eta = list2%zetapzetb
222 90774491 : EtaInv = list2%ZetaInv
223 90774491 : Zeta_C = list2%zeta
224 90774491 : Zeta_D = list2%zetb
225 90774491 : temp_CC = 0.0_dp
226 90774491 : temp_DD = 0.0_dp
227 193329177 : DO i = 1, nimages1
228 410218744 : P = list1%image_list(i)%P
229 102554686 : R1 = list1%image_list(i)%R
230 102554686 : S1234a = list1%image_list(i)%S1234
231 102554686 : pgf_max_1 = list1%image_list(i)%pgf_max
232 410218744 : ra = list1%image_list(i)%ra
233 410218744 : rb = list1%image_list(i)%rb
234 346138182 : DO j = 1, nimages2
235 152809005 : pgf_max_2 = list2%image_list(j)%pgf_max
236 152809005 : IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE
237 409691144 : Q = list2%image_list(j)%P
238 102422786 : R2 = list2%image_list(j)%R
239 102422786 : S1234b = list2%image_list(j)%S1234
240 409691144 : rc = list2%image_list(j)%ra
241 409691144 : rd = list2%image_list(j)%rb
242 :
243 102422786 : ZetapEtaInv = Zeta1 + Eta
244 102422786 : ZetapEtaInv = 1.0_dp/ZetapEtaInv
245 102422786 : Rho = Zeta1*Eta*ZetapEtaInv
246 102422786 : RhoInv = 1.0_dp/Rho
247 102422786 : S1234 = EXP(S1234a + S1234b)
248 102422786 : IF (do_periodic) THEN
249 199723748 : temp = P - Q
250 49930937 : PQ = pbc(temp, cell)
251 199723748 : shift = -PQ + temp
252 199723748 : temp_CC = rc + shift
253 199723748 : temp_DD = rd + shift
254 : END IF
255 :
256 3167693800 : DO k = 1, SIZE(neighbor_cells)
257 2962716328 : IF (do_periodic) THEN
258 11640897916 : C11 = temp_CC + neighbor_cells(k)%cell_r(:)
259 11640897916 : tmp_D = temp_DD + neighbor_cells(k)%cell_r(:)
260 : ELSE
261 52491849 : C11 = rc
262 52491849 : tmp_D = rd
263 : END IF
264 11850865312 : Q = (Zeta_C*C11 + Zeta_D*tmp_D)*EtaInv
265 2962716328 : 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 2962716328 : potential_parameter%potential_type == do_potential_short .OR. &
268 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
269 2891759312 : IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) CYCLE
270 : END IF
271 305788079 : IF (potential_parameter%potential_type == do_potential_TShPSC) THEN
272 1160520 : IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius*2.0_dp)**2) CYCLE
273 : END IF
274 305788079 : 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 305788079 : IF (nproducts > SIZE(product_list)) THEN
279 1529 : !$OMP ATOMIC READ
280 : tmp_i4 = pgf_product_list_size
281 1529 : tmp_i4 = MAX(pgf_product_list_size, (3*nproducts + 1)/2)
282 1529 : !$OMP ATOMIC WRITE
283 : pgf_product_list_size = tmp_i4
284 5082462 : ALLOCATE (tmp_product_list(SIZE(product_list)))
285 5006012 : tmp_product_list(:) = product_list
286 1529 : DEALLOCATE (product_list)
287 7587417 : ALLOCATE (product_list(tmp_i4))
288 5006012 : product_list(1:SIZE(tmp_product_list)) = tmp_product_list
289 1529 : DEALLOCATE (tmp_product_list)
290 : END IF
291 :
292 305788079 : T = Rho*rpq2
293 528132763 : SELECT CASE (potential_parameter%potential_type)
294 : CASE (do_potential_truncated)
295 222344684 : R = potential_parameter%cutoff_radius*SQRT(Rho)
296 222344684 : CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
297 222344684 : IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
298 222344684 : factor = 2.0_dp*Pi*RhoInv
299 : CASE (do_potential_TShPSC)
300 1160520 : R = potential_parameter%cutoff_radius*SQRT(Rho)
301 25531440 : product_list(nproducts)%Fm = 0.0_dp
302 1160520 : CALL trunc_CS_poly_n20(product_list(nproducts)%Fm(1), R, T, m_max)
303 1160520 : factor = 2.0_dp*Pi*RhoInv
304 : CASE (do_potential_coulomb)
305 67791700 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
306 67791700 : factor = 2.0_dp*Pi*RhoInv
307 : CASE (do_potential_short)
308 6542109 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
309 6542109 : omega2 = potential_parameter%omega**2
310 6542109 : omega_corr2 = omega2/(omega2 + Rho)
311 6542109 : omega_corr = SQRT(omega_corr2)
312 6542109 : T = T*omega_corr2
313 6542109 : CALL fgamma(m_max, T, Fm)
314 6542109 : tmp = -omega_corr
315 29786308 : DO l = 1, m_max + 1
316 23244199 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp
317 29786308 : tmp = tmp*omega_corr2
318 : END DO
319 6542109 : factor = 2.0_dp*Pi*RhoInv
320 : CASE (do_potential_long)
321 1092506 : omega2 = potential_parameter%omega**2
322 1092506 : omega_corr2 = omega2/(omega2 + Rho)
323 1092506 : omega_corr = SQRT(omega_corr2)
324 1092506 : T = T*omega_corr2
325 1092506 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
326 1092506 : tmp = omega_corr
327 3846061 : DO l = 1, m_max + 1
328 2753555 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*tmp
329 3846061 : tmp = tmp*omega_corr2
330 : END DO
331 1092506 : factor = 2.0_dp*Pi*RhoInv
332 : CASE (do_potential_mix_cl)
333 830676 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
334 830676 : omega2 = potential_parameter%omega**2
335 830676 : omega_corr2 = omega2/(omega2 + Rho)
336 830676 : omega_corr = SQRT(omega_corr2)
337 830676 : T = T*omega_corr2
338 830676 : CALL fgamma(m_max, T, Fm)
339 830676 : tmp = omega_corr
340 3047709 : DO l = 1, m_max + 1
341 : product_list(nproducts)%Fm(l) = &
342 : product_list(nproducts)%Fm(l)*potential_parameter%scale_coulomb &
343 2217033 : + Fm(l)*tmp*potential_parameter%scale_longrange
344 3047709 : tmp = tmp*omega_corr2
345 : END DO
346 830676 : factor = 2.0_dp*Pi*RhoInv
347 : CASE (do_potential_mix_cl_trunc)
348 :
349 : ! truncated
350 5944270 : R = potential_parameter%cutoff_radius*SQRT(rho)
351 5944270 : CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
352 5944270 : IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
353 :
354 : ! Coulomb
355 5944270 : CALL fgamma(m_max, T, Fm)
356 :
357 23027917 : 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 23027917 : Fm(l)*potential_parameter%scale_longrange
361 : END DO
362 :
363 : ! longrange
364 5944270 : omega2 = potential_parameter%omega**2
365 5944270 : omega_corr2 = omega2/(omega2 + Rho)
366 5944270 : omega_corr = SQRT(omega_corr2)
367 5944270 : T = T*omega_corr2
368 5944270 : CALL fgamma(m_max, T, Fm)
369 5944270 : tmp = omega_corr
370 23027917 : DO l = 1, m_max + 1
371 17083647 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
372 23027917 : tmp = tmp*omega_corr2
373 : END DO
374 5944270 : factor = 2.0_dp*Pi*RhoInv
375 :
376 : CASE (do_potential_gaussian)
377 0 : omega2 = potential_parameter%omega**2
378 0 : T = -omega2*T/(Rho + omega2)
379 0 : tmp = 1.0_dp
380 0 : DO l = 1, m_max + 1
381 0 : product_list(nproducts)%Fm(l) = EXP(T)*tmp
382 0 : tmp = tmp*omega2/(Rho + omega2)
383 : END DO
384 0 : factor = (Pi/(Rho + omega2))**(1.5_dp)
385 : CASE (do_potential_mix_lg)
386 29282 : omega2 = potential_parameter%omega**2
387 29282 : omega_corr2 = omega2/(omega2 + Rho)
388 29282 : omega_corr = SQRT(omega_corr2)
389 29282 : T = T*omega_corr2
390 29282 : CALL fgamma(m_max, T, Fm)
391 29282 : tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange
392 122452 : DO l = 1, m_max + 1
393 93170 : Fm(l) = Fm(l)*tmp
394 122452 : tmp = tmp*omega_corr2
395 : END DO
396 : T = Rho*rpq2
397 29282 : T = -omega2*T/(Rho + omega2)
398 29282 : tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
399 122452 : DO l = 1, m_max + 1
400 93170 : product_list(nproducts)%Fm(l) = EXP(T)*tmp + Fm(l)
401 122452 : tmp = tmp*omega2/(Rho + omega2)
402 : END DO
403 52332 : factor = 1.0_dp
404 : CASE (do_potential_id)
405 52332 : num = list1%zeta*list1%zetb
406 52332 : den = list1%zeta + list1%zetb
407 209328 : ssss = -num/den*SUM((ra - rb)**2)
408 :
409 52332 : num = den*Zeta_C
410 52332 : den = den + Zeta_C
411 209328 : ssss = ssss - num/den*SUM((P - rc)**2)
412 :
413 209328 : G(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + Zeta_C*rc(:))/den
414 52332 : num = den*Zeta_D
415 52332 : den = den + Zeta_D
416 209328 : ssss = ssss - num/den*SUM((G - rd)**2)
417 :
418 1151304 : product_list(nproducts)%Fm(:) = EXP(ssss)
419 52332 : factor = 1.0_dp
420 305840411 : IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
421 : END SELECT
422 :
423 305788079 : tmp = (Pi*ZetapEtaInv)**3
424 305788079 : factor = factor*S1234*SQRT(tmp)
425 :
426 1095499048 : DO l = 1, m_max + 1
427 1095499048 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*factor
428 : END DO
429 :
430 1223152316 : W = (Zeta1*P + Eta*Q)*ZetapEtaInv
431 1223152316 : product_list(nproducts)%ra = ra
432 1223152316 : product_list(nproducts)%rb = rb
433 1223152316 : product_list(nproducts)%rc = C11
434 1223152316 : product_list(nproducts)%rd = tmp_D
435 305788079 : product_list(nproducts)%ZetapEtaInv = ZetapEtaInv
436 305788079 : product_list(nproducts)%Rho = Rho
437 305788079 : product_list(nproducts)%RhoInv = RhoInv
438 1223152316 : product_list(nproducts)%P = P
439 1223152316 : product_list(nproducts)%Q = Q
440 1223152316 : product_list(nproducts)%W = W
441 1223152316 : product_list(nproducts)%AB = ra - rb
442 4032889570 : product_list(nproducts)%CD = C11 - tmp_D
443 : END DO
444 : END DO
445 : END DO
446 :
447 90774491 : 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 20682340 : 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 20682340 : 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 76673721 : nimages = 0
491 : ! ** inner loop may never be reached
492 20682340 : nelements = npgfa*npgfb
493 156680124 : DO i = 1, SIZE(neighbor_cells)
494 135997784 : IF (do_periodic) THEN
495 493315976 : im_B = rb + neighbor_cells(i)%cell_r(:)
496 : ELSE
497 12668790 : im_B = rb
498 : END IF
499 543991136 : AB = ra - im_B
500 135997784 : rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
501 135997784 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
502 : element_counter = 0
503 68534547 : DO ipgf = 1, npgfa
504 136295035 : DO jpgf = 1, npgfb
505 67760488 : element_counter = element_counter + 1
506 67760488 : pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
507 67760488 : IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
508 : CYCLE
509 : END IF
510 56958581 : nimages(element_counter) = nimages(element_counter) + 1
511 56958581 : list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
512 227834324 : list(element_counter)%image_list(nimages(element_counter))%ra = ra
513 227834324 : list(element_counter)%image_list(nimages(element_counter))%rb = im_B
514 56958581 : list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
515 :
516 56958581 : Zeta_A = zeta(ipgf)
517 56958581 : Zeta_B = zetb(jpgf)
518 56958581 : Zeta1 = Zeta_A + Zeta_B
519 56958581 : ZetaInv = 1.0_dp/Zeta1
520 :
521 56958581 : IF (nimages(element_counter) == 1) THEN
522 49552051 : list(element_counter)%ipgf = ipgf
523 49552051 : list(element_counter)%jpgf = jpgf
524 49552051 : list(element_counter)%zetaInv = ZetaInv
525 49552051 : list(element_counter)%zetapzetb = Zeta1
526 49552051 : list(element_counter)%zeta = Zeta_A
527 49552051 : list(element_counter)%zetb = Zeta_B
528 : END IF
529 :
530 56958581 : list(element_counter)%image_list(nimages(element_counter))%S1234 = (-Zeta_A*Zeta_B*ZetaInv*rab2)
531 227834324 : 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 56958581 : MAX(0.0_dp, R_pgf(jpgf, ipgf)%x(1)*rab2 + R_pgf(jpgf, ipgf)%x(2))
534 227834324 : list(element_counter)%image_list(nimages(element_counter))%ra = ra
535 227834324 : list(element_counter)%image_list(nimages(element_counter))%rb = im_B
536 56958581 : list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
537 280720941 : list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
538 : END DO
539 : END DO
540 156680124 : nelements = MAX(nelements, element_counter)
541 : END DO
542 76673721 : DO j = 1, nelements
543 76673721 : list(j)%nimages = nimages(j)
544 : END DO
545 : ! ** Remove unused elements
546 :
547 : element_counter = 0
548 76673721 : DO j = 1, nelements
549 55991381 : IF (list(j)%nimages == 0) CYCLE
550 49552051 : element_counter = element_counter + 1
551 49552051 : list(element_counter)%nimages = list(j)%nimages
552 49552051 : list(element_counter)%zetapzetb = list(j)%zetapzetb
553 49552051 : list(element_counter)%ZetaInv = list(j)%ZetaInv
554 49552051 : list(element_counter)%zeta = list(j)%zeta
555 49552051 : list(element_counter)%zetb = list(j)%zetb
556 49552051 : list(element_counter)%ipgf = list(j)%ipgf
557 49552051 : list(element_counter)%jpgf = list(j)%jpgf
558 127192972 : DO i = 1, list(j)%nimages
559 112949962 : list(element_counter)%image_list(i) = list(j)%image_list(i)
560 : END DO
561 : END DO
562 :
563 20682340 : nelements = element_counter
564 :
565 20682340 : 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 54605694 : SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
589 4079466 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
590 4079466 : 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 4079466 : n_element = 0
618 4079466 : nset_ij = 0
619 :
620 8158932 : DO iatom = i_start, i_end
621 12238398 : DO jatom = j_start, j_end
622 4079466 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
623 :
624 3918648 : ikind = kind_of(iatom)
625 3918648 : nseta = basis_parameter(ikind)%nset
626 15674592 : ra = particle_set(iatom)%r(:)
627 :
628 3918648 : IF (jatom < iatom) CYCLE
629 3918648 : jkind = kind_of(jatom)
630 3918648 : nsetb = basis_parameter(jkind)%nset
631 15674592 : rb = particle_set(jatom)%r(:)
632 :
633 3918648 : IF (do_periodic) THEN
634 5252008 : temp = rb - ra
635 1313002 : pbc_B = pbc(temp, cell)
636 5252008 : B11 = ra + pbc_B
637 1313002 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
638 : ELSE
639 2605646 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
640 2605646 : B11 = rb ! ra - rb
641 : END IF
642 3918648 : 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 3908888 : n_element = n_element + 1
646 11726664 : list%elements(n_element)%pair = [iatom, jatom]
647 11726664 : list%elements(n_element)%kind_pair = [ikind, jkind]
648 15635552 : list%elements(n_element)%r1 = ra
649 15635552 : list%elements(n_element)%r2 = B11
650 3908888 : list%elements(n_element)%dist2 = rab2
651 : ! build a list of guaranteed overlapping sets
652 3908888 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
653 12987578 : DO iset = 1, nseta
654 38202240 : DO jset = 1, nsetb
655 25214662 : 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 23841832 : nset_ij = nset_ij + 1
658 81977016 : set_list(nset_ij)%pair = [iset, jset]
659 : END DO
660 : END DO
661 8158932 : list%elements(n_element)%set_bounds(2) = nset_ij
662 : END DO
663 : END DO
664 :
665 4079466 : list%n_element = n_element
666 :
667 4079466 : 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 5934 : SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
684 5934 : 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 90298 : atomic_pair_list = .FALSE.
708 :
709 23996 : DO iblock = 1, SIZE(blocks)
710 18062 : iatom_start = blocks(iblock)%istart
711 18062 : iatom_end = blocks(iblock)%iend
712 90298 : DO jblock = 1, SIZE(blocks)
713 66302 : jatom_start = blocks(jblock)%istart
714 66302 : jatom_end = blocks(jblock)%iend
715 :
716 150666 : DO iatom = iatom_start, iatom_end
717 66302 : ikind = kind_of(iatom)
718 66302 : nseta = basis_parameter(ikind)%nset
719 265208 : ra = particle_set(iatom)%r(:)
720 198906 : DO jatom = jatom_start, jatom_end
721 66302 : IF (jatom < iatom) CYCLE
722 42182 : jkind = kind_of(jatom)
723 42182 : nsetb = basis_parameter(jkind)%nset
724 168728 : rb = particle_set(jatom)%r(:)
725 :
726 42182 : IF (do_periodic) THEN
727 50432 : temp = rb - ra
728 12608 : pbc_B = pbc(temp, cell)
729 50432 : B11 = ra + pbc_B
730 12608 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
731 : ELSE
732 29574 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
733 29574 : B11 = rb ! ra - rb
734 : END IF
735 42182 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
736 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) CYCLE
737 :
738 41724 : atomic_pair_list(jatom, iatom) = .TRUE.
739 132604 : atomic_pair_list(iatom, jatom) = .TRUE.
740 : END DO
741 : END DO
742 : END DO
743 : END DO
744 :
745 5934 : 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 595568 : SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
770 2994 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
771 2994 : 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 2994 : n_element = 0
801 2994 : nset_ij = 0
802 :
803 2994 : my_skip_atom_symmetry = .FALSE.
804 2994 : IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry
805 :
806 6076 : DO iatom = i_start, i_end
807 9454 : DO jatom = j_start, j_end
808 3378 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
809 :
810 3378 : ikind = kind_of(iatom)
811 3378 : nseta = basis_parameter(ikind)%nset
812 13512 : ra = particle_set(iatom)%r(:)
813 :
814 3378 : IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) CYCLE
815 3230 : jkind = kind_of(jatom)
816 3230 : nsetb = basis_parameter(jkind)%nset
817 12920 : rb = particle_set(jatom)%r(:)
818 :
819 3230 : IF (do_periodic) THEN
820 0 : temp = rb - ra
821 0 : pbc_B = pbc(temp, cell)
822 0 : B11 = ra + pbc_B
823 0 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
824 : ELSE
825 3230 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
826 3230 : B11 = rb ! ra - rb
827 : END IF
828 3230 : 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 3230 : n_element = n_element + 1
832 9690 : list%elements(n_element)%pair = [iatom, jatom]
833 9690 : list%elements(n_element)%kind_pair = [ikind, jkind]
834 12920 : list%elements(n_element)%r1 = ra
835 12920 : list%elements(n_element)%r2 = B11
836 3230 : list%elements(n_element)%dist2 = rab2
837 : ! build a list of guaranteed overlapping sets
838 3230 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
839 13988 : DO iset = 1, nseta
840 52586 : DO jset = 1, nsetb
841 38598 : 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 38526 : nset_ij = nset_ij + 1
844 126408 : set_list(nset_ij)%pair = [iset, jset]
845 : END DO
846 : END DO
847 6460 : list%elements(n_element)%set_bounds(2) = nset_ij
848 : END DO
849 : END DO
850 :
851 2994 : list%n_element = n_element
852 :
853 2994 : END SUBROUTINE build_pair_list_mp2
854 :
855 0 : END MODULE hfx_pair_list_methods
|
