Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines 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 :
45 : ! an initial estimate for the size of the product list
46 : INTEGER, SAVE :: pgf_product_list_size = 128
47 :
48 : !***
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief ...
54 : !> \param list1 ...
55 : !> \param list2 ...
56 : !> \param product_list ...
57 : !> \param nproducts ...
58 : !> \param log10_pmax ...
59 : !> \param log10_eps_schwarz ...
60 : !> \param neighbor_cells ...
61 : !> \param cell ...
62 : !> \param potential_parameter ...
63 : !> \param m_max ...
64 : !> \param do_periodic ...
65 : ! **************************************************************************************************
66 74816522 : SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
67 : log10_pmax, log10_eps_schwarz, neighbor_cells, &
68 : cell, potential_parameter, m_max, do_periodic)
69 :
70 : TYPE(hfx_pgf_list) :: list1, list2
71 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
72 : DIMENSION(:), INTENT(INOUT) :: product_list
73 : INTEGER, INTENT(OUT) :: nproducts
74 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
75 : TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
76 : TYPE(cell_type), POINTER :: cell
77 : TYPE(hfx_potential_type) :: potential_parameter
78 : INTEGER, INTENT(IN) :: m_max
79 : LOGICAL, INTENT(IN) :: do_periodic
80 :
81 : INTEGER :: i, j, k, l, nimages1, nimages2, tmp_i4
82 : LOGICAL :: use_gamma
83 : REAL(dp) :: C11(3), den, Eta, EtaInv, factor, Fm(prim_data_f_size), G(3), num, omega2, &
84 : omega_corr, omega_corr2, P(3), pgf_max_1, pgf_max_2, PQ(3), Q(3), R, R1, R2, ra(3), &
85 : rb(3), rc(3), rd(3), Rho, RhoInv, rpq2, S1234, S1234a, S1234b, shift(3), ssss, T, &
86 : temp(3), temp_CC(3), temp_DD(3), tmp, tmp_D(3), W(3), Zeta1, Zeta_C, Zeta_D, ZetapEtaInv
87 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
88 74816522 : DIMENSION(:) :: tmp_product_list
89 :
90 74816522 : nimages1 = list1%nimages
91 74816522 : nimages2 = list2%nimages
92 74816522 : nproducts = 0
93 74816522 : Zeta1 = list1%zetapzetb
94 74816522 : Eta = list2%zetapzetb
95 74816522 : EtaInv = list2%ZetaInv
96 74816522 : Zeta_C = list2%zeta
97 74816522 : Zeta_D = list2%zetb
98 74816522 : temp_CC = 0.0_dp
99 74816522 : temp_DD = 0.0_dp
100 159599649 : DO i = 1, nimages1
101 339132508 : P = list1%image_list(i)%P
102 84783127 : R1 = list1%image_list(i)%R
103 84783127 : S1234a = list1%image_list(i)%S1234
104 84783127 : pgf_max_1 = list1%image_list(i)%pgf_max
105 339132508 : ra = list1%image_list(i)%ra
106 339132508 : rb = list1%image_list(i)%rb
107 283362132 : DO j = 1, nimages2
108 123762483 : pgf_max_2 = list2%image_list(j)%pgf_max
109 123762483 : IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE
110 339901672 : Q = list2%image_list(j)%P
111 84975418 : R2 = list2%image_list(j)%R
112 84975418 : S1234b = list2%image_list(j)%S1234
113 339901672 : rc = list2%image_list(j)%ra
114 339901672 : rd = list2%image_list(j)%rb
115 :
116 84975418 : ZetapEtaInv = Zeta1 + Eta
117 84975418 : ZetapEtaInv = 1.0_dp/ZetapEtaInv
118 84975418 : Rho = Zeta1*Eta*ZetapEtaInv
119 84975418 : RhoInv = 1.0_dp/Rho
120 84975418 : S1234 = EXP(S1234a + S1234b)
121 84975418 : IF (do_periodic) THEN
122 152664568 : temp = P - Q
123 38166142 : PQ = pbc(temp, cell)
124 152664568 : shift = -PQ + temp
125 152664568 : temp_CC = rc + shift
126 152664568 : temp_DD = rd + shift
127 : END IF
128 :
129 2350508837 : DO k = 1, SIZE(neighbor_cells)
130 2180750292 : IF (do_periodic) THEN
131 8535764064 : C11 = temp_CC + neighbor_cells(k)%cell_r(:)
132 8535764064 : tmp_D = temp_DD + neighbor_cells(k)%cell_r(:)
133 : ELSE
134 46809276 : C11 = rc
135 46809276 : tmp_D = rd
136 : END IF
137 8723001168 : Q = (Zeta_C*C11 + Zeta_D*tmp_D)*EtaInv
138 2180750292 : rpq2 = (P(1) - Q(1))**2 + (P(2) - Q(2))**2 + (P(3) - Q(3))**2
139 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
140 2180750292 : potential_parameter%potential_type == do_potential_short .OR. &
141 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
142 2122995357 : IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) CYCLE
143 : END IF
144 185254862 : IF (potential_parameter%potential_type == do_potential_TShPSC) THEN
145 1160608 : IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius*2.0_dp)**2) CYCLE
146 : END IF
147 185254862 : nproducts = nproducts + 1
148 :
149 : ! allocate size as needed,
150 : ! updating the global size estimate to make this a rare event in longer simulations
151 185254862 : IF (nproducts > SIZE(product_list)) THEN
152 1311 : !$OMP ATOMIC READ
153 : tmp_i4 = pgf_product_list_size
154 1311 : tmp_i4 = MAX(pgf_product_list_size, (3*nproducts + 1)/2)
155 1311 : !$OMP ATOMIC WRITE
156 : pgf_product_list_size = tmp_i4
157 2702601 : ALLOCATE (tmp_product_list(SIZE(product_list)))
158 2637051 : tmp_product_list(:) = product_list
159 1311 : DEALLOCATE (product_list)
160 4022802 : ALLOCATE (product_list(tmp_i4))
161 2637051 : product_list(1:SIZE(tmp_product_list)) = tmp_product_list
162 1311 : DEALLOCATE (tmp_product_list)
163 : END IF
164 :
165 185254862 : T = Rho*rpq2
166 300276358 : SELECT CASE (potential_parameter%potential_type)
167 : CASE (do_potential_truncated)
168 115021496 : R = potential_parameter%cutoff_radius*SQRT(Rho)
169 115021496 : CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
170 115021496 : IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
171 115021496 : factor = 2.0_dp*Pi*RhoInv
172 : CASE (do_potential_TShPSC)
173 1160608 : R = potential_parameter%cutoff_radius*SQRT(Rho)
174 25533376 : product_list(nproducts)%Fm = 0.0_dp
175 1160608 : CALL trunc_CS_poly_n20(product_list(nproducts)%Fm(1), R, T, m_max)
176 1160608 : factor = 2.0_dp*Pi*RhoInv
177 : CASE (do_potential_coulomb)
178 54594650 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
179 54594650 : factor = 2.0_dp*Pi*RhoInv
180 : CASE (do_potential_short)
181 6538556 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
182 6538556 : omega2 = potential_parameter%omega**2
183 6538556 : omega_corr2 = omega2/(omega2 + Rho)
184 6538556 : omega_corr = SQRT(omega_corr2)
185 6538556 : T = T*omega_corr2
186 6538556 : CALL fgamma(m_max, T, Fm)
187 6538556 : tmp = -omega_corr
188 29777583 : DO l = 1, m_max + 1
189 23239027 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp
190 29777583 : tmp = tmp*omega_corr2
191 : END DO
192 6538556 : factor = 2.0_dp*Pi*RhoInv
193 : CASE (do_potential_long)
194 1087576 : omega2 = potential_parameter%omega**2
195 1087576 : omega_corr2 = omega2/(omega2 + Rho)
196 1087576 : omega_corr = SQRT(omega_corr2)
197 1087576 : T = T*omega_corr2
198 1087576 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
199 1087576 : tmp = omega_corr
200 3835774 : DO l = 1, m_max + 1
201 2748198 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*tmp
202 3835774 : tmp = tmp*omega_corr2
203 : END DO
204 1087576 : factor = 2.0_dp*Pi*RhoInv
205 : CASE (do_potential_mix_cl)
206 830487 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
207 830487 : omega2 = potential_parameter%omega**2
208 830487 : omega_corr2 = omega2/(omega2 + Rho)
209 830487 : omega_corr = SQRT(omega_corr2)
210 830487 : T = T*omega_corr2
211 830487 : CALL fgamma(m_max, T, Fm)
212 830487 : tmp = omega_corr
213 3047189 : DO l = 1, m_max + 1
214 : product_list(nproducts)%Fm(l) = &
215 : product_list(nproducts)%Fm(l)*potential_parameter%scale_coulomb &
216 2216702 : + Fm(l)*tmp*potential_parameter%scale_longrange
217 3047189 : tmp = tmp*omega_corr2
218 : END DO
219 830487 : factor = 2.0_dp*Pi*RhoInv
220 : CASE (do_potential_mix_cl_trunc)
221 :
222 : ! truncated
223 5939875 : R = potential_parameter%cutoff_radius*SQRT(rho)
224 5939875 : CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
225 5939875 : IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
226 :
227 : ! Coulomb
228 5939875 : CALL fgamma(m_max, T, Fm)
229 :
230 23020373 : DO l = 1, m_max + 1
231 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)* &
232 : (potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
233 23020373 : Fm(l)*potential_parameter%scale_longrange
234 : END DO
235 :
236 : ! longrange
237 5939875 : omega2 = potential_parameter%omega**2
238 5939875 : omega_corr2 = omega2/(omega2 + Rho)
239 5939875 : omega_corr = SQRT(omega_corr2)
240 5939875 : T = T*omega_corr2
241 5939875 : CALL fgamma(m_max, T, Fm)
242 5939875 : tmp = omega_corr
243 23020373 : DO l = 1, m_max + 1
244 17080498 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
245 23020373 : tmp = tmp*omega_corr2
246 : END DO
247 5939875 : factor = 2.0_dp*Pi*RhoInv
248 :
249 : CASE (do_potential_gaussian)
250 0 : omega2 = potential_parameter%omega**2
251 0 : T = -omega2*T/(Rho + omega2)
252 0 : tmp = 1.0_dp
253 0 : DO l = 1, m_max + 1
254 0 : product_list(nproducts)%Fm(l) = EXP(T)*tmp
255 0 : tmp = tmp*omega2/(Rho + omega2)
256 : END DO
257 0 : factor = (Pi/(Rho + omega2))**(1.5_dp)
258 : CASE (do_potential_mix_lg)
259 29282 : omega2 = potential_parameter%omega**2
260 29282 : omega_corr2 = omega2/(omega2 + Rho)
261 29282 : omega_corr = SQRT(omega_corr2)
262 29282 : T = T*omega_corr2
263 29282 : CALL fgamma(m_max, T, Fm)
264 29282 : tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange
265 122452 : DO l = 1, m_max + 1
266 93170 : Fm(l) = Fm(l)*tmp
267 122452 : tmp = tmp*omega_corr2
268 : END DO
269 : T = Rho*rpq2
270 29282 : T = -omega2*T/(Rho + omega2)
271 29282 : tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
272 122452 : DO l = 1, m_max + 1
273 93170 : product_list(nproducts)%Fm(l) = EXP(T)*tmp + Fm(l)
274 122452 : tmp = tmp*omega2/(Rho + omega2)
275 : END DO
276 52332 : factor = 1.0_dp
277 : CASE (do_potential_id)
278 52332 : num = list1%zeta*list1%zetb
279 52332 : den = list1%zeta + list1%zetb
280 209328 : ssss = -num/den*SUM((ra - rb)**2)
281 :
282 52332 : num = den*Zeta_C
283 52332 : den = den + Zeta_C
284 209328 : ssss = ssss - num/den*SUM((P - rc)**2)
285 :
286 209328 : G(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + Zeta_C*rc(:))/den
287 52332 : num = den*Zeta_D
288 52332 : den = den + Zeta_D
289 209328 : ssss = ssss - num/den*SUM((G - rd)**2)
290 :
291 1151304 : product_list(nproducts)%Fm(:) = EXP(ssss)
292 52332 : factor = 1.0_dp
293 185307194 : IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
294 : END SELECT
295 :
296 185254862 : tmp = (Pi*ZetapEtaInv)**3
297 185254862 : factor = factor*S1234*SQRT(tmp)
298 :
299 699086109 : DO l = 1, m_max + 1
300 699086109 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*factor
301 : END DO
302 :
303 741019448 : W = (Zeta1*P + Eta*Q)*ZetapEtaInv
304 741019448 : product_list(nproducts)%ra = ra
305 741019448 : product_list(nproducts)%rb = rb
306 741019448 : product_list(nproducts)%rc = C11
307 741019448 : product_list(nproducts)%rd = tmp_D
308 185254862 : product_list(nproducts)%ZetapEtaInv = ZetapEtaInv
309 185254862 : product_list(nproducts)%Rho = Rho
310 185254862 : product_list(nproducts)%RhoInv = RhoInv
311 741019448 : product_list(nproducts)%P = P
312 741019448 : product_list(nproducts)%Q = Q
313 741019448 : product_list(nproducts)%W = W
314 741019448 : product_list(nproducts)%AB = ra - rb
315 2860277361 : product_list(nproducts)%CD = C11 - tmp_D
316 : END DO
317 : END DO
318 : END DO
319 :
320 74816522 : END SUBROUTINE build_pgf_product_list
321 :
322 : ! **************************************************************************************************
323 : !> \brief ...
324 : !> \param npgfa ...
325 : !> \param npgfb ...
326 : !> \param list ...
327 : !> \param zeta ...
328 : !> \param zetb ...
329 : !> \param screen1 ...
330 : !> \param screen2 ...
331 : !> \param pgf ...
332 : !> \param R_pgf ...
333 : !> \param log10_pmax ...
334 : !> \param log10_eps_schwarz ...
335 : !> \param ra ...
336 : !> \param rb ...
337 : !> \param nelements ...
338 : !> \param neighbor_cells ...
339 : !> \param nimages ...
340 : !> \param do_periodic ...
341 : ! **************************************************************************************************
342 17254192 : SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, &
343 : log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
344 17254192 : neighbor_cells, nimages, do_periodic)
345 : INTEGER, INTENT(IN) :: npgfa, npgfb
346 : TYPE(hfx_pgf_list), DIMENSION(npgfa*npgfb) :: list
347 : REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
348 : REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
349 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2)
350 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
351 : POINTER :: pgf, R_pgf
352 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), &
353 : rb(3)
354 : INTEGER, INTENT(OUT) :: nelements
355 : TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
356 : INTEGER :: nimages(npgfa*npgfb)
357 : LOGICAL, INTENT(IN) :: do_periodic
358 :
359 : INTEGER :: element_counter, i, ipgf, j, jpgf
360 : REAL(dp) :: AB(3), im_B(3), pgf_max, rab2, Zeta1, &
361 : Zeta_A, Zeta_B, ZetaInv
362 :
363 62017986 : nimages = 0
364 : ! ** inner loop may never be reached
365 17254192 : nelements = npgfa*npgfb
366 137146092 : DO i = 1, SIZE(neighbor_cells)
367 119891900 : IF (do_periodic) THEN
368 433813008 : im_B = rb + neighbor_cells(i)%cell_r(:)
369 : ELSE
370 11438648 : im_B = rb
371 : END IF
372 479567600 : AB = ra - im_B
373 119891900 : rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
374 119891900 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
375 : element_counter = 0
376 56657680 : DO ipgf = 1, npgfa
377 111549948 : DO jpgf = 1, npgfb
378 54892268 : element_counter = element_counter + 1
379 54892268 : pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
380 54892268 : IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
381 : CYCLE
382 : END IF
383 46320113 : nimages(element_counter) = nimages(element_counter) + 1
384 46320113 : list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
385 185280452 : list(element_counter)%image_list(nimages(element_counter))%ra = ra
386 185280452 : list(element_counter)%image_list(nimages(element_counter))%rb = im_B
387 46320113 : list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
388 :
389 46320113 : Zeta_A = zeta(ipgf)
390 46320113 : Zeta_B = zetb(jpgf)
391 46320113 : Zeta1 = Zeta_A + Zeta_B
392 46320113 : ZetaInv = 1.0_dp/Zeta1
393 :
394 46320113 : IF (nimages(element_counter) == 1) THEN
395 40187382 : list(element_counter)%ipgf = ipgf
396 40187382 : list(element_counter)%jpgf = jpgf
397 40187382 : list(element_counter)%zetaInv = ZetaInv
398 40187382 : list(element_counter)%zetapzetb = Zeta1
399 40187382 : list(element_counter)%zeta = Zeta_A
400 40187382 : list(element_counter)%zetb = Zeta_B
401 : END IF
402 :
403 46320113 : list(element_counter)%image_list(nimages(element_counter))%S1234 = (-Zeta_A*Zeta_B*ZetaInv*rab2)
404 185280452 : list(element_counter)%image_list(nimages(element_counter))%P = (Zeta_A*ra + Zeta_B*im_B)*ZetaInv
405 : list(element_counter)%image_list(nimages(element_counter))%R = &
406 46320113 : MAX(0.0_dp, R_pgf(jpgf, ipgf)%x(1)*rab2 + R_pgf(jpgf, ipgf)%x(2))
407 185280452 : list(element_counter)%image_list(nimages(element_counter))%ra = ra
408 185280452 : list(element_counter)%image_list(nimages(element_counter))%rb = im_B
409 46320113 : list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
410 228422139 : list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
411 : END DO
412 : END DO
413 137146092 : nelements = MAX(nelements, element_counter)
414 : END DO
415 62017986 : DO j = 1, nelements
416 62017986 : list(j)%nimages = nimages(j)
417 : END DO
418 : ! ** Remove unused elements
419 :
420 : element_counter = 0
421 62017986 : DO j = 1, nelements
422 44763794 : IF (list(j)%nimages == 0) CYCLE
423 40187382 : element_counter = element_counter + 1
424 40187382 : list(element_counter)%nimages = list(j)%nimages
425 40187382 : list(element_counter)%zetapzetb = list(j)%zetapzetb
426 40187382 : list(element_counter)%ZetaInv = list(j)%ZetaInv
427 40187382 : list(element_counter)%zeta = list(j)%zeta
428 40187382 : list(element_counter)%zetb = list(j)%zetb
429 40187382 : list(element_counter)%ipgf = list(j)%ipgf
430 40187382 : list(element_counter)%jpgf = list(j)%jpgf
431 103761687 : DO i = 1, list(j)%nimages
432 91083907 : list(element_counter)%image_list(i) = list(j)%image_list(i)
433 : END DO
434 : END DO
435 :
436 17254192 : nelements = element_counter
437 :
438 17254192 : END SUBROUTINE build_pair_list_pgf
439 :
440 : ! **************************************************************************************************
441 : !> \brief ...
442 : !> \param natom ...
443 : !> \param list ...
444 : !> \param set_list ...
445 : !> \param i_start ...
446 : !> \param i_end ...
447 : !> \param j_start ...
448 : !> \param j_end ...
449 : !> \param kind_of ...
450 : !> \param basis_parameter ...
451 : !> \param particle_set ...
452 : !> \param do_periodic ...
453 : !> \param coeffs_set ...
454 : !> \param coeffs_kind ...
455 : !> \param coeffs_kind_max0 ...
456 : !> \param log10_eps_schwarz ...
457 : !> \param cell ...
458 : !> \param pmax_blocks ...
459 : !> \param atomic_pair_list ...
460 : ! **************************************************************************************************
461 3026779942 : SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
462 2868410 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
463 2868410 : pmax_blocks, atomic_pair_list)
464 :
465 : INTEGER, INTENT(IN) :: natom
466 : TYPE(pair_list_type), INTENT(OUT) :: list
467 : TYPE(pair_set_list_type), DIMENSION(:), &
468 : INTENT(OUT) :: set_list
469 : INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
470 : kind_of(*)
471 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
472 : POINTER :: basis_parameter
473 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
474 : POINTER :: particle_set
475 : LOGICAL, INTENT(IN) :: do_periodic
476 : TYPE(hfx_screen_coeff_type), &
477 : DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
478 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
479 : INTENT(IN) :: coeffs_kind
480 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
481 : TYPE(cell_type), POINTER :: cell
482 : REAL(dp), INTENT(IN) :: pmax_blocks
483 : LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
484 :
485 : INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
486 : n_element, nset_ij, nseta, nsetb
487 : REAL(KIND=dp) :: rab2
488 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
489 :
490 2868410 : n_element = 0
491 2868410 : nset_ij = 0
492 :
493 5736820 : DO iatom = i_start, i_end
494 8605230 : DO jatom = j_start, j_end
495 2868410 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
496 :
497 2707592 : ikind = kind_of(iatom)
498 2707592 : nseta = basis_parameter(ikind)%nset
499 10830368 : ra = particle_set(iatom)%r(:)
500 :
501 2707592 : IF (jatom < iatom) CYCLE
502 2707592 : jkind = kind_of(jatom)
503 2707592 : nsetb = basis_parameter(jkind)%nset
504 10830368 : rb = particle_set(jatom)%r(:)
505 :
506 2707592 : IF (do_periodic) THEN
507 4507696 : temp = rb - ra
508 1126924 : pbc_B = pbc(temp, cell)
509 4507696 : B11 = ra + pbc_B
510 1126924 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
511 : ELSE
512 1580668 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
513 1580668 : B11 = rb ! ra - rb
514 : END IF
515 2707592 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
516 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
517 :
518 2697832 : n_element = n_element + 1
519 8093496 : list%elements(n_element)%pair = (/iatom, jatom/)
520 8093496 : list%elements(n_element)%kind_pair = (/ikind, jkind/)
521 10791328 : list%elements(n_element)%r1 = ra
522 10791328 : list%elements(n_element)%r2 = B11
523 2697832 : list%elements(n_element)%dist2 = rab2
524 : ! build a list of guaranteed overlapping sets
525 2697832 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
526 9755836 : DO iset = 1, nseta
527 30506544 : DO jset = 1, nsetb
528 20750708 : IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
529 : coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
530 19433094 : nset_ij = nset_ij + 1
531 66674900 : set_list(nset_ij)%pair = (/iset, jset/)
532 : END DO
533 : END DO
534 5736820 : list%elements(n_element)%set_bounds(2) = nset_ij
535 : END DO
536 : END DO
537 :
538 2868410 : list%n_element = n_element
539 :
540 2868410 : END SUBROUTINE build_pair_list
541 :
542 : ! **************************************************************************************************
543 : !> \brief ...
544 : !> \param natom ...
545 : !> \param atomic_pair_list ...
546 : !> \param kind_of ...
547 : !> \param basis_parameter ...
548 : !> \param particle_set ...
549 : !> \param do_periodic ...
550 : !> \param coeffs_kind ...
551 : !> \param coeffs_kind_max0 ...
552 : !> \param log10_eps_schwarz ...
553 : !> \param cell ...
554 : !> \param blocks ...
555 : ! **************************************************************************************************
556 5174 : SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
557 5174 : do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
558 : blocks)
559 : INTEGER, INTENT(IN) :: natom
560 : LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list
561 : INTEGER, INTENT(IN) :: kind_of(*)
562 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
563 : POINTER :: basis_parameter
564 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
565 : POINTER :: particle_set
566 : LOGICAL, INTENT(IN) :: do_periodic
567 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
568 : INTENT(IN) :: coeffs_kind
569 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
570 : TYPE(cell_type), POINTER :: cell
571 : TYPE(hfx_block_range_type), DIMENSION(:), &
572 : INTENT(IN), POINTER :: blocks
573 :
574 : INTEGER :: iatom, iatom_end, iatom_start, iblock, &
575 : ikind, jatom, jatom_end, jatom_start, &
576 : jblock, jkind, nseta, nsetb
577 : REAL(KIND=dp) :: rab2
578 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
579 :
580 77238 : atomic_pair_list = .FALSE.
581 :
582 20834 : DO iblock = 1, SIZE(blocks)
583 15660 : iatom_start = blocks(iblock)%istart
584 15660 : iatom_end = blocks(iblock)%iend
585 77238 : DO jblock = 1, SIZE(blocks)
586 56404 : jatom_start = blocks(jblock)%istart
587 56404 : jatom_end = blocks(jblock)%iend
588 :
589 128468 : DO iatom = iatom_start, iatom_end
590 56404 : ikind = kind_of(iatom)
591 56404 : nseta = basis_parameter(ikind)%nset
592 225616 : ra = particle_set(iatom)%r(:)
593 169212 : DO jatom = jatom_start, jatom_end
594 56404 : IF (jatom < iatom) CYCLE
595 36032 : jkind = kind_of(jatom)
596 36032 : nsetb = basis_parameter(jkind)%nset
597 144128 : rb = particle_set(jatom)%r(:)
598 :
599 36032 : IF (do_periodic) THEN
600 47280 : temp = rb - ra
601 11820 : pbc_B = pbc(temp, cell)
602 47280 : B11 = ra + pbc_B
603 11820 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
604 : ELSE
605 24212 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
606 24212 : B11 = rb ! ra - rb
607 : END IF
608 36032 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
609 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) CYCLE
610 :
611 35574 : atomic_pair_list(jatom, iatom) = .TRUE.
612 112808 : atomic_pair_list(iatom, jatom) = .TRUE.
613 : END DO
614 : END DO
615 : END DO
616 : END DO
617 :
618 5174 : END SUBROUTINE build_atomic_pair_list
619 :
620 : ! **************************************************************************************************
621 : !> \brief ...
622 : !> \param natom ...
623 : !> \param list ...
624 : !> \param set_list ...
625 : !> \param i_start ...
626 : !> \param i_end ...
627 : !> \param j_start ...
628 : !> \param j_end ...
629 : !> \param kind_of ...
630 : !> \param basis_parameter ...
631 : !> \param particle_set ...
632 : !> \param do_periodic ...
633 : !> \param coeffs_set ...
634 : !> \param coeffs_kind ...
635 : !> \param coeffs_kind_max0 ...
636 : !> \param log10_eps_schwarz ...
637 : !> \param cell ...
638 : !> \param pmax_blocks ...
639 : !> \param atomic_pair_list ...
640 : !> \param skip_atom_symmetry ...
641 : ! **************************************************************************************************
642 595568 : SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
643 2994 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
644 2994 : pmax_blocks, atomic_pair_list, skip_atom_symmetry)
645 :
646 : INTEGER, INTENT(IN) :: natom
647 : TYPE(pair_list_type_mp2) :: list
648 : TYPE(pair_set_list_type), DIMENSION(:), &
649 : INTENT(OUT) :: set_list
650 : INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
651 : kind_of(*)
652 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
653 : POINTER :: basis_parameter
654 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
655 : POINTER :: particle_set
656 : LOGICAL, INTENT(IN) :: do_periodic
657 : TYPE(hfx_screen_coeff_type), &
658 : DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
659 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
660 : INTENT(IN) :: coeffs_kind
661 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
662 : TYPE(cell_type), POINTER :: cell
663 : REAL(dp), INTENT(IN) :: pmax_blocks
664 : LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
665 : LOGICAL, INTENT(IN), OPTIONAL :: skip_atom_symmetry
666 :
667 : INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
668 : n_element, nset_ij, nseta, nsetb
669 : LOGICAL :: my_skip_atom_symmetry
670 : REAL(KIND=dp) :: rab2
671 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
672 :
673 2994 : n_element = 0
674 2994 : nset_ij = 0
675 :
676 2994 : my_skip_atom_symmetry = .FALSE.
677 2994 : IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry
678 :
679 6076 : DO iatom = i_start, i_end
680 9454 : DO jatom = j_start, j_end
681 3378 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
682 :
683 3378 : ikind = kind_of(iatom)
684 3378 : nseta = basis_parameter(ikind)%nset
685 13512 : ra = particle_set(iatom)%r(:)
686 :
687 3378 : IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) CYCLE
688 3230 : jkind = kind_of(jatom)
689 3230 : nsetb = basis_parameter(jkind)%nset
690 12920 : rb = particle_set(jatom)%r(:)
691 :
692 3230 : IF (do_periodic) THEN
693 0 : temp = rb - ra
694 0 : pbc_B = pbc(temp, cell)
695 0 : B11 = ra + pbc_B
696 0 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
697 : ELSE
698 3230 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
699 3230 : B11 = rb ! ra - rb
700 : END IF
701 3230 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
702 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
703 :
704 3230 : n_element = n_element + 1
705 9690 : list%elements(n_element)%pair = (/iatom, jatom/)
706 9690 : list%elements(n_element)%kind_pair = (/ikind, jkind/)
707 12920 : list%elements(n_element)%r1 = ra
708 12920 : list%elements(n_element)%r2 = B11
709 3230 : list%elements(n_element)%dist2 = rab2
710 : ! build a list of guaranteed overlapping sets
711 3230 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
712 13988 : DO iset = 1, nseta
713 52586 : DO jset = 1, nsetb
714 38598 : IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
715 : coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
716 38526 : nset_ij = nset_ij + 1
717 126408 : set_list(nset_ij)%pair = (/iset, jset/)
718 : END DO
719 : END DO
720 6460 : list%elements(n_element)%set_bounds(2) = nset_ij
721 : END DO
722 : END DO
723 :
724 2994 : list%n_element = n_element
725 :
726 2994 : END SUBROUTINE build_pair_list_mp2
727 :
728 : END MODULE hfx_pair_list_methods
|
