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, min_image_gap_sq
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 92588479 : 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 : REAL(dp), PARAMETER :: r_cutoff_eps = 1.0E-12_dp
209 :
210 : INTEGER :: i, j, k, k_max, l, nimages1, nimages2, &
211 : tmp_i4
212 : LOGICAL :: has_r_cutoff, use_gamma
213 : REAL(dp) :: C11(3), den, Eta, EtaInv, factor, Fm(prim_data_f_size), G(3), min_img_delta, &
214 : num, omega2, omega_corr, omega_corr2, P(3), pgf_max_1, pgf_max_2, PQ(3), pq2, Q(3), R, &
215 : R1, R2, r_cutoff_sq, r_cutoff_sq_guard, ra(3), rb(3), rc(3), rcut_eff, rd(3), Rho, &
216 : RhoInv, rpq2, S1234, S1234a, S1234b, shift(3), ssss, T, temp(3), temp_CC(3), temp_DD(3), &
217 : tmp, tmp_D(3), W(3), Zeta1, Zeta_C, Zeta_D, ZetapEtaInv
218 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
219 92588479 : DIMENSION(:) :: tmp_product_list
220 :
221 92588479 : nimages1 = list1%nimages
222 92588479 : nimages2 = list2%nimages
223 92588479 : nproducts = 0
224 92588479 : Zeta1 = list1%zetapzetb
225 92588479 : Eta = list2%zetapzetb
226 92588479 : EtaInv = list2%ZetaInv
227 92588479 : Zeta_C = list2%zeta
228 92588479 : Zeta_D = list2%zetb
229 92588479 : temp_CC = 0.0_dp
230 92588479 : temp_DD = 0.0_dp
231 92588479 : r_cutoff_sq = 0.0_dp
232 92588479 : r_cutoff_sq_guard = 0.0_dp
233 104475808 : SELECT CASE (potential_parameter%potential_type)
234 : CASE (do_potential_truncated, do_potential_short, do_potential_mix_cl_trunc)
235 11887329 : has_r_cutoff = .TRUE.
236 11887329 : rcut_eff = potential_parameter%cutoff_radius
237 : CASE (do_potential_TShPSC)
238 1253728 : has_r_cutoff = .TRUE.
239 1253728 : rcut_eff = potential_parameter%cutoff_radius*2.0_dp
240 : CASE DEFAULT
241 : has_r_cutoff = .FALSE.
242 92588479 : rcut_eff = 0.0_dp
243 : END SELECT
244 196953673 : DO i = 1, nimages1
245 417460776 : P = list1%image_list(i)%P
246 104365194 : R1 = list1%image_list(i)%R
247 104365194 : S1234a = list1%image_list(i)%S1234
248 104365194 : pgf_max_1 = list1%image_list(i)%pgf_max
249 417460776 : ra = list1%image_list(i)%ra
250 417460776 : rb = list1%image_list(i)%rb
251 351645024 : DO j = 1, nimages2
252 154691351 : pgf_max_2 = list2%image_list(j)%pgf_max
253 154691351 : IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE
254 416817092 : Q = list2%image_list(j)%P
255 104204273 : R2 = list2%image_list(j)%R
256 104204273 : S1234b = list2%image_list(j)%S1234
257 416817092 : rc = list2%image_list(j)%ra
258 416817092 : rd = list2%image_list(j)%rb
259 :
260 104204273 : IF (has_r_cutoff) THEN
261 32612694 : r_cutoff_sq = (R1 + R2 + rcut_eff)**2
262 32612694 : r_cutoff_sq_guard = r_cutoff_sq*(1.0_dp + r_cutoff_eps)
263 : END IF
264 104204273 : k_max = SIZE(neighbor_cells)
265 104204273 : IF (do_periodic) THEN
266 201056652 : temp = P - Q
267 50264163 : PQ = pbc(temp, cell)
268 200902918 : IF (has_r_cutoff .AND. cell%orthorhombic .AND. ALL(cell%perd == 1) .AND. k_max >= 1) THEN
269 102865072 : IF (ALL(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
270 25716268 : pq2 = PQ(1)**2 + PQ(2)**2 + PQ(3)**2
271 25716268 : IF (pq2 > r_cutoff_sq_guard) CYCLE
272 25549582 : min_img_delta = min_image_gap_sq(PQ, cell%hmat)
273 25549582 : IF (min_img_delta >= 0.0_dp .AND. pq2 + min_img_delta > r_cutoff_sq_guard) k_max = 1
274 : END IF
275 : END IF
276 200389908 : shift = -PQ + temp
277 200389908 : temp_CC = rc + shift
278 200389908 : temp_DD = rd + shift
279 : END IF
280 :
281 104037587 : ZetapEtaInv = Zeta1 + Eta
282 104037587 : ZetapEtaInv = 1.0_dp/ZetapEtaInv
283 104037587 : Rho = Zeta1*Eta*ZetapEtaInv
284 104037587 : RhoInv = 1.0_dp/Rho
285 104037587 : S1234 = EXP(S1234a + S1234b)
286 :
287 2627621190 : DO k = 1, k_max
288 2419218409 : IF (do_periodic) THEN
289 2365278299 : IF (has_r_cutoff .AND. k_max > 1) THEN
290 2334829497 : IF ((PQ(1) - neighbor_cells(k)%cell_r(1))**2 + &
291 : (PQ(2) - neighbor_cells(k)%cell_r(2))**2 + &
292 : (PQ(3) - neighbor_cells(k)%cell_r(3))**2 > r_cutoff_sq_guard) CYCLE
293 : END IF
294 1014087992 : C11 = temp_CC + neighbor_cells(k)%cell_r(:)
295 1014087992 : tmp_D = temp_DD + neighbor_cells(k)%cell_r(:)
296 : ELSE
297 53940110 : C11 = rc
298 53940110 : tmp_D = rd
299 : END IF
300 1229848432 : Q = (Zeta_C*C11 + Zeta_D*tmp_D)*EtaInv
301 307462108 : rpq2 = (P(1) - Q(1))**2 + (P(2) - Q(2))**2 + (P(3) - Q(3))**2
302 307462108 : IF (has_r_cutoff) THEN
303 235870529 : IF (rpq2 > r_cutoff_sq) CYCLE
304 : END IF
305 307462091 : nproducts = nproducts + 1
306 :
307 : ! allocate size as needed,
308 : ! updating the global size estimate to make this a rare event in longer simulations
309 307462091 : IF (nproducts > SIZE(product_list)) THEN
310 580 : !$OMP ATOMIC READ
311 : tmp_i4 = pgf_product_list_size
312 580 : tmp_i4 = MAX(pgf_product_list_size, (3*nproducts + 1)/2)
313 580 : !$OMP ATOMIC WRITE
314 : pgf_product_list_size = tmp_i4
315 2130040 : ALLOCATE (tmp_product_list(SIZE(product_list)))
316 2101040 : tmp_product_list(:) = product_list
317 580 : DEALLOCATE (product_list)
318 3181296 : ALLOCATE (product_list(tmp_i4))
319 2101040 : product_list(1:SIZE(tmp_product_list)) = tmp_product_list
320 580 : DEALLOCATE (tmp_product_list)
321 : END IF
322 :
323 307462091 : T = Rho*rpq2
324 529680173 : SELECT CASE (potential_parameter%potential_type)
325 : CASE (do_potential_truncated)
326 222218082 : R = potential_parameter%cutoff_radius*SQRT(Rho)
327 222218082 : CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
328 222218082 : IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
329 222218082 : factor = 2.0_dp*Pi*RhoInv
330 : CASE (do_potential_TShPSC)
331 1160672 : R = potential_parameter%cutoff_radius*SQRT(Rho)
332 25534784 : product_list(nproducts)%Fm = 0.0_dp
333 1160672 : CALL trunc_CS_poly_n20(product_list(nproducts)%Fm(1), R, T, m_max)
334 1160672 : factor = 2.0_dp*Pi*RhoInv
335 : CASE (do_potential_coulomb)
336 69074202 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
337 69074202 : factor = 2.0_dp*Pi*RhoInv
338 : CASE (do_potential_short)
339 6542868 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
340 6542868 : omega2 = potential_parameter%omega**2
341 6542868 : omega_corr2 = omega2/(omega2 + Rho)
342 6542868 : omega_corr = SQRT(omega_corr2)
343 6542868 : T = T*omega_corr2
344 6542868 : CALL fgamma(m_max, T, Fm)
345 6542868 : tmp = -omega_corr
346 29789508 : DO l = 1, m_max + 1
347 23246640 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp
348 29789508 : tmp = tmp*omega_corr2
349 : END DO
350 6542868 : factor = 2.0_dp*Pi*RhoInv
351 : CASE (do_potential_long)
352 1092553 : omega2 = potential_parameter%omega**2
353 1092553 : omega_corr2 = omega2/(omega2 + Rho)
354 1092553 : omega_corr = SQRT(omega_corr2)
355 1092553 : T = T*omega_corr2
356 1092553 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
357 1092553 : tmp = omega_corr
358 3846231 : DO l = 1, m_max + 1
359 2753678 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*tmp
360 3846231 : tmp = tmp*omega_corr2
361 : END DO
362 1092553 : factor = 2.0_dp*Pi*RhoInv
363 : CASE (do_potential_mix_cl)
364 1343210 : CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
365 1343210 : omega2 = potential_parameter%omega**2
366 1343210 : omega_corr2 = omega2/(omega2 + Rho)
367 1343210 : omega_corr = SQRT(omega_corr2)
368 1343210 : T = T*omega_corr2
369 1343210 : CALL fgamma(m_max, T, Fm)
370 1343210 : tmp = omega_corr
371 5215847 : DO l = 1, m_max + 1
372 : product_list(nproducts)%Fm(l) = &
373 : product_list(nproducts)%Fm(l)*potential_parameter%scale_coulomb &
374 3872637 : + Fm(l)*tmp*potential_parameter%scale_longrange
375 5215847 : tmp = tmp*omega_corr2
376 : END DO
377 1343210 : factor = 2.0_dp*Pi*RhoInv
378 : CASE (do_potential_mix_cl_trunc)
379 :
380 : ! truncated
381 5948890 : R = potential_parameter%cutoff_radius*SQRT(rho)
382 5948890 : CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
383 5948890 : IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
384 :
385 : ! Coulomb
386 5948890 : CALL fgamma(m_max, T, Fm)
387 :
388 23047065 : DO l = 1, m_max + 1
389 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)* &
390 : (potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
391 23047065 : Fm(l)*potential_parameter%scale_longrange
392 : END DO
393 :
394 : ! longrange
395 5948890 : omega2 = potential_parameter%omega**2
396 5948890 : omega_corr2 = omega2/(omega2 + Rho)
397 5948890 : omega_corr = SQRT(omega_corr2)
398 5948890 : T = T*omega_corr2
399 5948890 : CALL fgamma(m_max, T, Fm)
400 5948890 : tmp = omega_corr
401 23047065 : DO l = 1, m_max + 1
402 17098175 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
403 23047065 : tmp = tmp*omega_corr2
404 : END DO
405 5948890 : factor = 2.0_dp*Pi*RhoInv
406 :
407 : CASE (do_potential_gaussian)
408 0 : omega2 = potential_parameter%omega**2
409 0 : T = -omega2*T/(Rho + omega2)
410 0 : tmp = 1.0_dp
411 0 : DO l = 1, m_max + 1
412 0 : product_list(nproducts)%Fm(l) = EXP(T)*tmp
413 0 : tmp = tmp*omega2/(Rho + omega2)
414 : END DO
415 0 : factor = (Pi/(Rho + omega2))**(1.5_dp)
416 : CASE (do_potential_mix_lg)
417 29282 : omega2 = potential_parameter%omega**2
418 29282 : omega_corr2 = omega2/(omega2 + Rho)
419 29282 : omega_corr = SQRT(omega_corr2)
420 29282 : T = T*omega_corr2
421 29282 : CALL fgamma(m_max, T, Fm)
422 29282 : tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange
423 122452 : DO l = 1, m_max + 1
424 93170 : Fm(l) = Fm(l)*tmp
425 122452 : tmp = tmp*omega_corr2
426 : END DO
427 : T = Rho*rpq2
428 29282 : T = -omega2*T/(Rho + omega2)
429 29282 : tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
430 122452 : DO l = 1, m_max + 1
431 93170 : product_list(nproducts)%Fm(l) = EXP(T)*tmp + Fm(l)
432 122452 : tmp = tmp*omega2/(Rho + omega2)
433 : END DO
434 52332 : factor = 1.0_dp
435 : CASE (do_potential_id)
436 52332 : num = list1%zeta*list1%zetb
437 52332 : den = list1%zeta + list1%zetb
438 209328 : ssss = -num/den*SUM((ra - rb)**2)
439 :
440 52332 : num = den*Zeta_C
441 52332 : den = den + Zeta_C
442 209328 : ssss = ssss - num/den*SUM((P - rc)**2)
443 :
444 209328 : G(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + Zeta_C*rc(:))/den
445 52332 : num = den*Zeta_D
446 52332 : den = den + Zeta_D
447 209328 : ssss = ssss - num/den*SUM((G - rd)**2)
448 :
449 1151304 : product_list(nproducts)%Fm(:) = EXP(ssss)
450 52332 : factor = 1.0_dp
451 307514423 : IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
452 : END SELECT
453 :
454 307462091 : tmp = (Pi*ZetapEtaInv)**3
455 307462091 : factor = factor*S1234*SQRT(tmp)
456 :
457 1102349504 : DO l = 1, m_max + 1
458 1102349504 : product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*factor
459 : END DO
460 :
461 1229848364 : W = (Zeta1*P + Eta*Q)*ZetapEtaInv
462 1229848364 : product_list(nproducts)%ra = ra
463 1229848364 : product_list(nproducts)%rb = rb
464 1229848364 : product_list(nproducts)%rc = C11
465 1229848364 : product_list(nproducts)%rd = tmp_D
466 307462091 : product_list(nproducts)%ZetapEtaInv = ZetapEtaInv
467 307462091 : product_list(nproducts)%Rho = Rho
468 307462091 : product_list(nproducts)%RhoInv = RhoInv
469 1229848364 : product_list(nproducts)%P = P
470 1229848364 : product_list(nproducts)%Q = Q
471 1229848364 : product_list(nproducts)%W = W
472 1229848364 : product_list(nproducts)%AB = ra - rb
473 3496296033 : product_list(nproducts)%CD = C11 - tmp_D
474 : END DO
475 : END DO
476 : END DO
477 :
478 92588479 : END SUBROUTINE build_pgf_product_list
479 :
480 : ! **************************************************************************************************
481 : !> \brief ...
482 : !> \param npgfa ...
483 : !> \param npgfb ...
484 : !> \param list ...
485 : !> \param zeta ...
486 : !> \param zetb ...
487 : !> \param screen1 ...
488 : !> \param screen2 ...
489 : !> \param pgf ...
490 : !> \param R_pgf ...
491 : !> \param log10_pmax ...
492 : !> \param log10_eps_schwarz ...
493 : !> \param ra ...
494 : !> \param rb ...
495 : !> \param nelements ...
496 : !> \param neighbor_cells ...
497 : !> \param nimages ...
498 : !> \param do_periodic ...
499 : !> \param cell ...
500 : ! **************************************************************************************************
501 20917836 : SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, &
502 : log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
503 20917836 : neighbor_cells, nimages, do_periodic, cell)
504 : INTEGER, INTENT(IN) :: npgfa, npgfb
505 : TYPE(hfx_pgf_list), DIMENSION(npgfa*npgfb) :: list
506 : REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
507 : REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
508 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2)
509 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
510 : POINTER :: pgf, R_pgf
511 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), &
512 : rb(3)
513 : INTEGER, INTENT(OUT) :: nelements
514 : TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
515 : INTEGER :: nimages(npgfa*npgfb)
516 : LOGICAL, INTENT(IN) :: do_periodic
517 : TYPE(cell_type), POINTER :: cell
518 :
519 : REAL(dp), PARAMETER :: min_img_eps = 1.0E-12_dp
520 :
521 : INTEGER :: element_counter, i, i_max, ipgf, j, jpgf
522 : REAL(dp) :: AB(3), im_B(3), min_img_delta, pgf_max, &
523 : rab2, Zeta1, Zeta_A, Zeta_B, ZetaInv
524 :
525 20917836 : i_max = SIZE(neighbor_cells)
526 20917836 : IF (do_periodic) THEN
527 8004402 : IF (cell%orthorhombic .AND. screen1(1) <= 0.0_dp .AND. i_max > 1) THEN
528 10961260 : IF (ALL(cell%perd == 1) .AND. ALL(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
529 6225784 : AB = ra - rb
530 1556446 : rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
531 1556446 : min_img_delta = min_image_gap_sq(AB, cell%hmat)
532 1556446 : IF (min_img_delta >= 0.0_dp) THEN
533 1556446 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
534 0 : nelements = 0
535 0 : RETURN
536 : END IF
537 1556446 : IF (screen1(1)*(rab2 + min_img_delta)*(1.0_dp - min_img_eps) + &
538 824058 : screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) i_max = 1
539 : END IF
540 : END IF
541 : END IF
542 : ELSE
543 12913434 : i_max = MIN(i_max, 1)
544 : END IF
545 77661623 : nimages = 0
546 : ! ** inner loop may never be reached
547 20917836 : nelements = npgfa*npgfb
548 115381922 : DO i = 1, i_max
549 94464086 : IF (do_periodic) THEN
550 326202608 : im_B = rb + neighbor_cells(i)%cell_r(:)
551 : ELSE
552 12913434 : im_B = rb
553 : END IF
554 377856344 : AB = ra - im_B
555 94464086 : rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
556 94464086 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
557 : element_counter = 0
558 69169424 : DO ipgf = 1, npgfa
559 137686166 : DO jpgf = 1, npgfb
560 68516742 : element_counter = element_counter + 1
561 68516742 : pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
562 68516742 : IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
563 : CYCLE
564 : END IF
565 57720747 : nimages(element_counter) = nimages(element_counter) + 1
566 57720747 : list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
567 230882988 : list(element_counter)%image_list(nimages(element_counter))%ra = ra
568 230882988 : list(element_counter)%image_list(nimages(element_counter))%rb = im_B
569 57720747 : list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
570 :
571 57720747 : Zeta_A = zeta(ipgf)
572 57720747 : Zeta_B = zetb(jpgf)
573 57720747 : Zeta1 = Zeta_A + Zeta_B
574 57720747 : ZetaInv = 1.0_dp/Zeta1
575 :
576 57720747 : IF (nimages(element_counter) == 1) THEN
577 50311536 : list(element_counter)%ipgf = ipgf
578 50311536 : list(element_counter)%jpgf = jpgf
579 50311536 : list(element_counter)%zetaInv = ZetaInv
580 50311536 : list(element_counter)%zetapzetb = Zeta1
581 50311536 : list(element_counter)%zeta = Zeta_A
582 50311536 : list(element_counter)%zetb = Zeta_B
583 : END IF
584 :
585 57720747 : list(element_counter)%image_list(nimages(element_counter))%S1234 = (-Zeta_A*Zeta_B*ZetaInv*rab2)
586 230882988 : list(element_counter)%image_list(nimages(element_counter))%P = (Zeta_A*ra + Zeta_B*im_B)*ZetaInv
587 : list(element_counter)%image_list(nimages(element_counter))%R = &
588 57720747 : MAX(0.0_dp, R_pgf(jpgf, ipgf)%x(1)*rab2 + R_pgf(jpgf, ipgf)%x(2))
589 230882988 : list(element_counter)%image_list(nimages(element_counter))%ra = ra
590 230882988 : list(element_counter)%image_list(nimages(element_counter))%rb = im_B
591 57720747 : list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
592 284159570 : list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
593 : END DO
594 : END DO
595 115381922 : nelements = MAX(nelements, element_counter)
596 : END DO
597 77661623 : DO j = 1, nelements
598 77661623 : list(j)%nimages = nimages(j)
599 : END DO
600 : ! ** Remove unused elements
601 :
602 : element_counter = 0
603 77661623 : DO j = 1, nelements
604 56743787 : IF (list(j)%nimages == 0) CYCLE
605 50311536 : element_counter = element_counter + 1
606 50311536 : list(element_counter)%nimages = list(j)%nimages
607 50311536 : list(element_counter)%zetapzetb = list(j)%zetapzetb
608 50311536 : list(element_counter)%ZetaInv = list(j)%ZetaInv
609 50311536 : list(element_counter)%zeta = list(j)%zeta
610 50311536 : list(element_counter)%zetb = list(j)%zetb
611 50311536 : list(element_counter)%ipgf = list(j)%ipgf
612 50311536 : list(element_counter)%jpgf = list(j)%jpgf
613 128950119 : DO i = 1, list(j)%nimages
614 114464534 : list(element_counter)%image_list(i) = list(j)%image_list(i)
615 : END DO
616 : END DO
617 :
618 20917836 : nelements = element_counter
619 :
620 : END SUBROUTINE build_pair_list_pgf
621 :
622 : ! **************************************************************************************************
623 : !> \brief ...
624 : !> \param natom ...
625 : !> \param list ...
626 : !> \param set_list ...
627 : !> \param i_start ...
628 : !> \param i_end ...
629 : !> \param j_start ...
630 : !> \param j_end ...
631 : !> \param kind_of ...
632 : !> \param basis_parameter ...
633 : !> \param particle_set ...
634 : !> \param do_periodic ...
635 : !> \param coeffs_set ...
636 : !> \param coeffs_kind ...
637 : !> \param coeffs_kind_max0 ...
638 : !> \param log10_eps_schwarz ...
639 : !> \param cell ...
640 : !> \param pmax_blocks ...
641 : !> \param atomic_pair_list ...
642 : ! **************************************************************************************************
643 52678956 : SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
644 3801644 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
645 3801644 : pmax_blocks, atomic_pair_list)
646 :
647 : INTEGER, INTENT(IN) :: natom
648 : TYPE(pair_list_type), INTENT(INOUT) :: list
649 : TYPE(pair_set_list_type), DIMENSION(:), &
650 : INTENT(OUT) :: set_list
651 : INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
652 : kind_of(*)
653 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
654 : POINTER :: basis_parameter
655 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
656 : POINTER :: particle_set
657 : LOGICAL, INTENT(IN) :: do_periodic
658 : TYPE(hfx_screen_coeff_type), &
659 : DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
660 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
661 : INTENT(IN) :: coeffs_kind
662 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
663 : TYPE(cell_type), POINTER :: cell
664 : REAL(dp), INTENT(IN) :: pmax_blocks
665 : LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
666 :
667 : INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
668 : n_element, nset_ij, nseta, nsetb
669 : REAL(KIND=dp) :: rab2
670 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
671 :
672 3801644 : n_element = 0
673 3801644 : nset_ij = 0
674 :
675 7603288 : DO iatom = i_start, i_end
676 11404932 : DO jatom = j_start, j_end
677 3801644 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
678 :
679 3667994 : ikind = kind_of(iatom)
680 3667994 : nseta = basis_parameter(ikind)%nset
681 14671976 : ra = particle_set(iatom)%r(:)
682 :
683 3667994 : IF (jatom < iatom) CYCLE
684 3667994 : jkind = kind_of(jatom)
685 3667994 : nsetb = basis_parameter(jkind)%nset
686 14671976 : rb = particle_set(jatom)%r(:)
687 :
688 3667994 : IF (do_periodic) THEN
689 4695560 : temp = rb - ra
690 1173890 : pbc_B = pbc(temp, cell)
691 4695560 : B11 = ra + pbc_B
692 1173890 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
693 : ELSE
694 2494104 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
695 2494104 : B11 = rb ! ra - rb
696 : END IF
697 3667994 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
698 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
699 :
700 3661524 : n_element = n_element + 1
701 10984572 : list%elements(n_element)%pair = [iatom, jatom]
702 10984572 : list%elements(n_element)%kind_pair = [ikind, jkind]
703 14646096 : list%elements(n_element)%r1 = ra
704 14646096 : list%elements(n_element)%r2 = B11
705 3661524 : list%elements(n_element)%dist2 = rab2
706 : ! build a list of guaranteed overlapping sets
707 3661524 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
708 12364554 : DO iset = 1, nseta
709 36779176 : DO jset = 1, nsetb
710 24414622 : IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
711 : coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
712 23210674 : nset_ij = nset_ij + 1
713 79539000 : set_list(nset_ij)%pair = [iset, jset]
714 : END DO
715 : END DO
716 7603288 : list%elements(n_element)%set_bounds(2) = nset_ij
717 : END DO
718 : END DO
719 :
720 3801644 : list%n_element = n_element
721 :
722 3801644 : END SUBROUTINE build_pair_list
723 :
724 : ! **************************************************************************************************
725 : !> \brief ...
726 : !> \param natom ...
727 : !> \param atomic_pair_list ...
728 : !> \param kind_of ...
729 : !> \param basis_parameter ...
730 : !> \param particle_set ...
731 : !> \param do_periodic ...
732 : !> \param coeffs_kind ...
733 : !> \param coeffs_kind_max0 ...
734 : !> \param log10_eps_schwarz ...
735 : !> \param cell ...
736 : !> \param blocks ...
737 : ! **************************************************************************************************
738 6238 : SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
739 6238 : do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
740 : blocks)
741 : INTEGER, INTENT(IN) :: natom
742 : LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list
743 : INTEGER, INTENT(IN) :: kind_of(*)
744 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
745 : POINTER :: basis_parameter
746 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
747 : POINTER :: particle_set
748 : LOGICAL, INTENT(IN) :: do_periodic
749 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
750 : INTENT(IN) :: coeffs_kind
751 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
752 : TYPE(cell_type), POINTER :: cell
753 : TYPE(hfx_block_range_type), DIMENSION(:), &
754 : INTENT(IN), POINTER :: blocks
755 :
756 : INTEGER :: iatom, iatom_end, iatom_start, iblock, &
757 : ikind, jatom, jatom_end, jatom_start, &
758 : jblock, jkind, nseta, nsetb
759 : REAL(KIND=dp) :: rab2
760 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
761 :
762 93482 : atomic_pair_list = .FALSE.
763 :
764 25080 : DO iblock = 1, SIZE(blocks)
765 18842 : iatom_start = blocks(iblock)%istart
766 18842 : iatom_end = blocks(iblock)%iend
767 93482 : DO jblock = 1, SIZE(blocks)
768 68402 : jatom_start = blocks(jblock)%istart
769 68402 : jatom_end = blocks(jblock)%iend
770 :
771 155646 : DO iatom = iatom_start, iatom_end
772 68402 : ikind = kind_of(iatom)
773 68402 : nseta = basis_parameter(ikind)%nset
774 273608 : ra = particle_set(iatom)%r(:)
775 205206 : DO jatom = jatom_start, jatom_end
776 68402 : IF (jatom < iatom) CYCLE
777 43622 : jkind = kind_of(jatom)
778 43622 : nsetb = basis_parameter(jkind)%nset
779 174488 : rb = particle_set(jatom)%r(:)
780 :
781 43622 : IF (do_periodic) THEN
782 50504 : temp = rb - ra
783 12626 : pbc_B = pbc(temp, cell)
784 50504 : B11 = ra + pbc_B
785 12626 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
786 : ELSE
787 30996 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
788 30996 : B11 = rb ! ra - rb
789 : END IF
790 43622 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
791 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) CYCLE
792 :
793 43166 : atomic_pair_list(jatom, iatom) = .TRUE.
794 136804 : atomic_pair_list(iatom, jatom) = .TRUE.
795 : END DO
796 : END DO
797 : END DO
798 : END DO
799 :
800 6238 : END SUBROUTINE build_atomic_pair_list
801 :
802 : ! **************************************************************************************************
803 : !> \brief ...
804 : !> \param natom ...
805 : !> \param list ...
806 : !> \param set_list ...
807 : !> \param i_start ...
808 : !> \param i_end ...
809 : !> \param j_start ...
810 : !> \param j_end ...
811 : !> \param kind_of ...
812 : !> \param basis_parameter ...
813 : !> \param particle_set ...
814 : !> \param do_periodic ...
815 : !> \param coeffs_set ...
816 : !> \param coeffs_kind ...
817 : !> \param coeffs_kind_max0 ...
818 : !> \param log10_eps_schwarz ...
819 : !> \param cell ...
820 : !> \param pmax_blocks ...
821 : !> \param atomic_pair_list ...
822 : !> \param skip_atom_symmetry ...
823 : ! **************************************************************************************************
824 595568 : SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
825 2994 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
826 2994 : pmax_blocks, atomic_pair_list, skip_atom_symmetry)
827 :
828 : INTEGER, INTENT(IN) :: natom
829 : TYPE(pair_list_type_mp2), INTENT(INOUT) :: list
830 : TYPE(pair_set_list_type), DIMENSION(:), &
831 : INTENT(OUT) :: set_list
832 : INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end, &
833 : kind_of(*)
834 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
835 : POINTER :: basis_parameter
836 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
837 : POINTER :: particle_set
838 : LOGICAL, INTENT(IN) :: do_periodic
839 : TYPE(hfx_screen_coeff_type), &
840 : DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
841 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
842 : INTENT(IN) :: coeffs_kind
843 : REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz
844 : TYPE(cell_type), POINTER :: cell
845 : REAL(dp), INTENT(IN) :: pmax_blocks
846 : LOGICAL, DIMENSION(natom, natom), INTENT(IN) :: atomic_pair_list
847 : LOGICAL, INTENT(IN), OPTIONAL :: skip_atom_symmetry
848 :
849 : INTEGER :: iatom, ikind, iset, jatom, jkind, jset, &
850 : n_element, nset_ij, nseta, nsetb
851 : LOGICAL :: my_skip_atom_symmetry
852 : REAL(KIND=dp) :: rab2
853 : REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp
854 :
855 2994 : n_element = 0
856 2994 : nset_ij = 0
857 :
858 2994 : my_skip_atom_symmetry = .FALSE.
859 2994 : IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry
860 :
861 6076 : DO iatom = i_start, i_end
862 9454 : DO jatom = j_start, j_end
863 3378 : IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
864 :
865 3378 : ikind = kind_of(iatom)
866 3378 : nseta = basis_parameter(ikind)%nset
867 13512 : ra = particle_set(iatom)%r(:)
868 :
869 3378 : IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) CYCLE
870 3230 : jkind = kind_of(jatom)
871 3230 : nsetb = basis_parameter(jkind)%nset
872 12920 : rb = particle_set(jatom)%r(:)
873 :
874 3230 : IF (do_periodic) THEN
875 0 : temp = rb - ra
876 0 : pbc_B = pbc(temp, cell)
877 0 : B11 = ra + pbc_B
878 0 : rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
879 : ELSE
880 3230 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
881 3230 : B11 = rb ! ra - rb
882 : END IF
883 3230 : IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
884 : coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
885 :
886 3230 : n_element = n_element + 1
887 9690 : list%elements(n_element)%pair = [iatom, jatom]
888 9690 : list%elements(n_element)%kind_pair = [ikind, jkind]
889 12920 : list%elements(n_element)%r1 = ra
890 12920 : list%elements(n_element)%r2 = B11
891 3230 : list%elements(n_element)%dist2 = rab2
892 : ! build a list of guaranteed overlapping sets
893 3230 : list%elements(n_element)%set_bounds(1) = nset_ij + 1
894 13988 : DO iset = 1, nseta
895 52586 : DO jset = 1, nsetb
896 38598 : IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
897 : coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
898 38526 : nset_ij = nset_ij + 1
899 126408 : set_list(nset_ij)%pair = [iset, jset]
900 : END DO
901 : END DO
902 6460 : list%elements(n_element)%set_bounds(2) = nset_ij
903 : END DO
904 : END DO
905 :
906 2994 : list%n_element = n_element
907 :
908 2994 : END SUBROUTINE build_pair_list_mp2
909 :
910 : ! **************************************************************************************************
911 : !> \brief Smallest increase of |r|**2 when r is moved to any other periodic image of an
912 : !> orthorhombic, fully periodic cell, i.e. min over n /= 0 of |r - hmat*n|**2 - |r|**2.
913 : !> Exact if r is a minimum-image vector (|r(i)| <= |hmat(i,i)|/2); negative otherwise.
914 : !> \param r displacement vector
915 : !> \param hmat cell matrix; its diagonal entries are signed cell-vector components
916 : !> \return the squared-distance gap to the nearest other periodic image
917 : ! **************************************************************************************************
918 27145334 : PURE FUNCTION min_image_gap_sq(r, hmat) RESULT(gap)
919 : REAL(dp), DIMENSION(3), INTENT(IN) :: r
920 : REAL(dp), DIMENSION(3, 3), INTENT(IN) :: hmat
921 : REAL(dp) :: gap
922 :
923 : REAL(dp), DIMENSION(3) :: cell_len
924 :
925 108581336 : cell_len = [ABS(hmat(1, 1)), ABS(hmat(2, 2)), ABS(hmat(3, 3))]
926 108581336 : gap = MINVAL(cell_len*(cell_len - 2.0_dp*ABS(r)))
927 :
928 27145334 : END FUNCTION min_image_gap_sq
929 :
930 0 : END MODULE hfx_pair_list_methods
|
