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 Prepare and evaluate screened four-center Gaussian integrals and derivatives.
10 : ! **************************************************************************************************
11 : MODULE four_center_integral_methods
12 : USE cell_types, ONLY: cell_type,&
13 : pbc
14 : USE four_center_integral_interface, ONLY: evaluate_deriv_eri,&
15 : evaluate_eri,&
16 : four_center_integral_context_type
17 : USE four_center_integral_preparation,ONLY: prepare_eri4_pair,&
18 : prepare_eri4_primitive,&
19 : prepare_eri4_product
20 : USE four_center_integral_types, ONLY: eri4_cell_type,&
21 : eri4_pgf_image_type,&
22 : eri4_pgf_pair_type,&
23 : eri4_pgf_product_type,&
24 : eri4_screen_coeff_type
25 : USE input_constants, ONLY: do_potential_TShPSC,&
26 : do_potential_mix_cl_trunc,&
27 : do_potential_short,&
28 : do_potential_truncated
29 : USE integral_library_types, ONLY: coulomb_operator_type
30 : USE kinds, ONLY: dp,&
31 : int_8
32 : USE orbital_pointers, ONLY: nco,&
33 : nso
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 : PRIVATE
38 :
39 : PUBLIC :: coulomb4, forces4, min_image_gap_sq, pgf_product_list_size
40 :
41 : ! Initial capacity estimate, shared by the per-thread product buffers.
42 : INTEGER, SAVE :: pgf_product_list_size = 128
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief ...
48 : !> \param list1 ...
49 : !> \param list2 ...
50 : !> \param product_list ...
51 : !> \param nproducts ...
52 : !> \param log10_pmax ...
53 : !> \param log10_eps_schwarz ...
54 : !> \param neighbor_cells ...
55 : !> \param cell ...
56 : !> \param potential_parameter ...
57 : !> \param m_max ...
58 : !> \param do_periodic ...
59 : ! **************************************************************************************************
60 90548112 : SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
61 : log10_pmax, log10_eps_schwarz, neighbor_cells, &
62 : cell, potential_parameter, m_max, do_periodic)
63 :
64 : TYPE(eri4_pgf_pair_type) :: list1, list2
65 : TYPE(eri4_pgf_product_type), ALLOCATABLE, &
66 : DIMENSION(:), INTENT(INOUT) :: product_list
67 : INTEGER, INTENT(OUT) :: nproducts
68 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
69 : TYPE(eri4_cell_type), DIMENSION(:), POINTER :: neighbor_cells
70 : TYPE(cell_type), POINTER :: cell
71 : TYPE(coulomb_operator_type) :: potential_parameter
72 : INTEGER, INTENT(IN) :: m_max
73 : LOGICAL, INTENT(IN) :: do_periodic
74 :
75 : REAL(dp), PARAMETER :: r_cutoff_eps = 1.0E-12_dp
76 :
77 : INTEGER :: i, j, k, k_max, nimages1, nimages2, &
78 : tmp_i4
79 : LOGICAL :: has_r_cutoff
80 : REAL(dp) :: C11(3), EtaInv, min_img_delta, P(3), pgf_max_1, pgf_max_2, PQ(3), pq2, Q(3), R1, &
81 : R2, r_cutoff_sq, r_cutoff_sq_guard, rc(3), rcut_eff, rd(3), rpq2, S1234, S1234a, S1234b, &
82 : shift(3), temp(3), temp_CC(3), temp_DD(3), tmp_D(3), Zeta_C, Zeta_D
83 : TYPE(eri4_pgf_image_type) :: shifted_image
84 : TYPE(eri4_pgf_product_type) :: product
85 : TYPE(eri4_pgf_product_type), ALLOCATABLE, &
86 90548112 : DIMENSION(:) :: tmp_product_list
87 :
88 90548112 : nimages1 = list1%nimages
89 90548112 : nimages2 = list2%nimages
90 90548112 : nproducts = 0
91 90548112 : IF (nimages1 == 0 .OR. nimages2 == 0) RETURN
92 90548112 : CALL prepare_eri4_product(list1, list2, m_max, product)
93 90548112 : EtaInv = list2%ZetaInv
94 90548112 : Zeta_C = list2%zeta
95 90548112 : Zeta_D = list2%zetb
96 90548112 : temp_CC = 0.0_dp
97 90548112 : temp_DD = 0.0_dp
98 90548112 : r_cutoff_sq = 0.0_dp
99 90548112 : r_cutoff_sq_guard = 0.0_dp
100 102457055 : SELECT CASE (potential_parameter%potential_type)
101 : CASE (do_potential_truncated, do_potential_short, do_potential_mix_cl_trunc)
102 11908943 : has_r_cutoff = .TRUE.
103 11908943 : rcut_eff = potential_parameter%cutoff_radius
104 : CASE (do_potential_TShPSC)
105 1253728 : has_r_cutoff = .TRUE.
106 1253728 : rcut_eff = potential_parameter%cutoff_radius*2.0_dp
107 : CASE DEFAULT
108 : has_r_cutoff = .FALSE.
109 90548112 : rcut_eff = 0.0_dp
110 : END SELECT
111 192892128 : DO i = 1, nimages1
112 409376064 : P = list1%image_list(i)%P
113 102344016 : R1 = list1%image_list(i)%R
114 102344016 : S1234a = list1%image_list(i)%S1234
115 102344016 : pgf_max_1 = list1%image_list(i)%pgf_max
116 345665876 : DO j = 1, nimages2
117 152773748 : pgf_max_2 = list2%image_list(j)%pgf_max
118 152773748 : IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE
119 410613544 : Q = list2%image_list(j)%P
120 102653386 : R2 = list2%image_list(j)%R
121 102653386 : S1234b = list2%image_list(j)%S1234
122 410613544 : rc = list2%image_list(j)%ra
123 410613544 : rd = list2%image_list(j)%rb
124 :
125 102653386 : IF (has_r_cutoff) THEN
126 32694262 : r_cutoff_sq = (R1 + R2 + rcut_eff)**2
127 32694262 : r_cutoff_sq_guard = r_cutoff_sq*(1.0_dp + r_cutoff_eps)
128 : END IF
129 102653386 : k_max = SIZE(neighbor_cells)
130 102653386 : IF (do_periodic) THEN
131 194610400 : temp = P - Q
132 48652600 : PQ = pbc(temp, cell)
133 194456116 : IF (has_r_cutoff .AND. cell%orthorhombic .AND. ALL(cell%perd == 1) .AND. k_max >= 1) THEN
134 103088020 : IF (ALL(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
135 25772005 : pq2 = PQ(1)**2 + PQ(2)**2 + PQ(3)**2
136 25772005 : IF (pq2 > r_cutoff_sq_guard) CYCLE
137 25601762 : min_img_delta = min_image_gap_sq(PQ, cell%hmat)
138 25601762 : IF (min_img_delta >= 0.0_dp .AND. pq2 + min_img_delta > r_cutoff_sq_guard) k_max = 1
139 : END IF
140 : END IF
141 193929428 : shift = -PQ + temp
142 193929428 : temp_CC = rc + shift
143 193929428 : temp_DD = rd + shift
144 : END IF
145 :
146 102483143 : S1234 = EXP(S1234a + S1234b)
147 :
148 2624978216 : DO k = 1, k_max
149 2420151057 : IF (do_periodic) THEN
150 2366150271 : IF (has_r_cutoff .AND. k_max > 1) THEN
151 2337342937 : IF ((PQ(1) - neighbor_cells(k)%cell_r(1))**2 + &
152 : (PQ(2) - neighbor_cells(k)%cell_r(2))**2 + &
153 : (PQ(3) - neighbor_cells(k)%cell_r(3))**2 > r_cutoff_sq_guard) CYCLE
154 : END IF
155 1007900524 : C11 = temp_CC + neighbor_cells(k)%cell_r(:)
156 1007900524 : tmp_D = temp_DD + neighbor_cells(k)%cell_r(:)
157 : ELSE
158 54000786 : C11 = rc
159 54000786 : tmp_D = rd
160 : END IF
161 1223903668 : Q = (Zeta_C*C11 + Zeta_D*tmp_D)*EtaInv
162 305975917 : rpq2 = (P(1) - Q(1))**2 + (P(2) - Q(2))**2 + (P(3) - Q(3))**2
163 305975917 : IF (has_r_cutoff) THEN
164 236016793 : IF (rpq2 > r_cutoff_sq) CYCLE
165 : END IF
166 305975900 : nproducts = nproducts + 1
167 :
168 : ! allocate size as needed,
169 : ! updating the global size estimate to make this a rare event in longer simulations
170 305975900 : IF (nproducts > SIZE(product_list)) THEN
171 582 : !$OMP ATOMIC READ
172 : tmp_i4 = pgf_product_list_size
173 582 : tmp_i4 = MAX(pgf_product_list_size, (3*nproducts + 1)/2)
174 582 : !$OMP ATOMIC WRITE
175 : pgf_product_list_size = tmp_i4
176 2135798 : ALLOCATE (tmp_product_list(SIZE(product_list)))
177 2105534 : tmp_product_list(:) = product_list
178 582 : DEALLOCATE (product_list)
179 3189304 : ALLOCATE (product_list(tmp_i4))
180 2105534 : product_list(1:SIZE(tmp_product_list)) = tmp_product_list
181 582 : DEALLOCATE (tmp_product_list)
182 : END IF
183 :
184 1223903600 : shifted_image%ra = C11
185 1223903600 : shifted_image%rb = tmp_D
186 1223903600 : shifted_image%P = Q
187 : CALL prepare_eri4_primitive(product, list1%image_list(i), shifted_image, S1234, rpq2, &
188 2572924805 : potential_parameter, product_list(nproducts))
189 : END DO
190 : END DO
191 : END DO
192 :
193 5795079168 : END SUBROUTINE build_pgf_product_list
194 :
195 : ! **************************************************************************************************
196 : !> \brief ...
197 : !> \param npgfa ...
198 : !> \param npgfb ...
199 : !> \param list ...
200 : !> \param zeta ...
201 : !> \param zetb ...
202 : !> \param screen1 ...
203 : !> \param screen2 ...
204 : !> \param pgf ...
205 : !> \param R_pgf ...
206 : !> \param log10_pmax ...
207 : !> \param log10_eps_schwarz ...
208 : !> \param ra ...
209 : !> \param rb ...
210 : !> \param nelements ...
211 : !> \param neighbor_cells ...
212 : !> \param nimages ...
213 : !> \param do_periodic ...
214 : !> \param cell ...
215 : ! **************************************************************************************************
216 20488502 : SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, &
217 : log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
218 20488502 : neighbor_cells, nimages, do_periodic, cell)
219 : INTEGER, INTENT(IN) :: npgfa, npgfb
220 : TYPE(eri4_pgf_pair_type), DIMENSION(npgfa*npgfb) :: list
221 : REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
222 : REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
223 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2)
224 : TYPE(eri4_screen_coeff_type), DIMENSION(:, :), &
225 : POINTER :: pgf, R_pgf
226 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), &
227 : rb(3)
228 : INTEGER, INTENT(OUT) :: nelements
229 : TYPE(eri4_cell_type), DIMENSION(:), POINTER :: neighbor_cells
230 : INTEGER :: nimages(npgfa*npgfb)
231 : LOGICAL, INTENT(IN) :: do_periodic
232 : TYPE(cell_type), POINTER :: cell
233 :
234 : REAL(dp), PARAMETER :: min_img_eps = 1.0E-12_dp
235 :
236 : INTEGER :: element_counter, i, i_max, ipgf, j, jpgf
237 : REAL(dp) :: AB(3), im_B(3), min_img_delta, pgf_max, &
238 : rab2
239 :
240 20488502 : i_max = SIZE(neighbor_cells)
241 20488502 : IF (do_periodic) THEN
242 7569374 : IF (cell%orthorhombic .AND. screen1(1) <= 0.0_dp .AND. i_max > 1) THEN
243 10973790 : IF (ALL(cell%perd == 1) .AND. ALL(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
244 6232944 : AB = ra - rb
245 1558236 : rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
246 1558236 : min_img_delta = min_image_gap_sq(AB, cell%hmat)
247 1558236 : IF (min_img_delta >= 0.0_dp) THEN
248 1558236 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
249 0 : nelements = 0
250 0 : RETURN
251 : END IF
252 1558236 : IF (screen1(1)*(rab2 + min_img_delta)*(1.0_dp - min_img_eps) + &
253 825064 : screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) i_max = 1
254 : END IF
255 : END IF
256 : END IF
257 : ELSE
258 12919128 : i_max = MIN(i_max, 1)
259 : END IF
260 75559037 : nimages = 0
261 : ! ** inner loop may never be reached
262 20488502 : nelements = npgfa*npgfb
263 114584778 : DO i = 1, i_max
264 94096276 : IF (do_periodic) THEN
265 324708592 : im_B = rb + neighbor_cells(i)%cell_r(:)
266 : ELSE
267 12919128 : im_B = rb
268 : END IF
269 376385104 : AB = ra - im_B
270 94096276 : rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
271 94096276 : IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
272 : element_counter = 0
273 67864133 : DO ipgf = 1, npgfa
274 134718564 : DO jpgf = 1, npgfb
275 66854431 : element_counter = element_counter + 1
276 66854431 : pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
277 66854431 : IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
278 : CYCLE
279 : END IF
280 56404230 : nimages(element_counter) = nimages(element_counter) + 1
281 56404230 : list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
282 : CALL prepare_eri4_pair(ra, im_B, zeta(ipgf), zetb(jpgf), rab2, &
283 56404230 : list(element_counter), list(element_counter)%image_list(nimages(element_counter)))
284 56404230 : IF (nimages(element_counter) == 1) THEN
285 48989710 : list(element_counter)%ipgf = ipgf
286 48989710 : list(element_counter)%jpgf = jpgf
287 : END IF
288 :
289 : list(element_counter)%image_list(nimages(element_counter))%R = &
290 56404230 : MAX(0.0_dp, R_pgf(jpgf, ipgf)%x(1)*rab2 + R_pgf(jpgf, ipgf)%x(2))
291 277667492 : list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
292 : END DO
293 : END DO
294 114584778 : nelements = MAX(nelements, element_counter)
295 : END DO
296 75559037 : DO j = 1, nelements
297 75559037 : list(j)%nimages = nimages(j)
298 : END DO
299 : ! ** Remove unused elements
300 :
301 : element_counter = 0
302 75559037 : DO j = 1, nelements
303 55070535 : IF (list(j)%nimages == 0) CYCLE
304 48989710 : element_counter = element_counter + 1
305 48989710 : list(element_counter)%nimages = list(j)%nimages
306 48989710 : list(element_counter)%zetapzetb = list(j)%zetapzetb
307 48989710 : list(element_counter)%ZetaInv = list(j)%ZetaInv
308 48989710 : list(element_counter)%zeta = list(j)%zeta
309 48989710 : list(element_counter)%zetb = list(j)%zetb
310 48989710 : list(element_counter)%ipgf = list(j)%ipgf
311 48989710 : list(element_counter)%jpgf = list(j)%jpgf
312 125882442 : DO i = 1, list(j)%nimages
313 111474765 : list(element_counter)%image_list(i) = list(j)%image_list(i)
314 : END DO
315 : END DO
316 :
317 20488502 : nelements = element_counter
318 :
319 : END SUBROUTINE build_pair_list_pgf
320 :
321 : ! **************************************************************************************************
322 : !> \brief Smallest increase of |r|**2 when r is moved to any other periodic image of an
323 : !> orthorhombic, fully periodic cell, i.e. min over n /= 0 of |r - hmat*n|**2 - |r|**2.
324 : !> Exact if r is a minimum-image vector (|r(i)| <= |hmat(i,i)|/2); negative otherwise.
325 : !> \param r displacement vector
326 : !> \param hmat cell matrix; its diagonal entries are signed cell-vector components
327 : !> \return the squared-distance gap to the nearest other periodic image
328 : ! **************************************************************************************************
329 27160062 : PURE FUNCTION min_image_gap_sq(r, hmat) RESULT(gap)
330 : REAL(dp), DIMENSION(3), INTENT(IN) :: r
331 : REAL(dp), DIMENSION(3, 3), INTENT(IN) :: hmat
332 : REAL(dp) :: gap
333 :
334 : REAL(dp), DIMENSION(3) :: cell_len
335 :
336 108640248 : cell_len = [ABS(hmat(1, 1)), ABS(hmat(2, 2)), ABS(hmat(3, 3))]
337 108640248 : gap = MINVAL(cell_len*(cell_len - 2.0_dp*ABS(r)))
338 :
339 27160062 : END FUNCTION min_image_gap_sq
340 :
341 : ! **************************************************************************************************
342 : !> \brief calculates two-electron integrals of a quartet/shell using the library
343 : !> lib_int in the periodic case
344 : !> \param context four-center integral context
345 : !> \param ra ...
346 : !> \param rb ...
347 : !> \param rc ...
348 : !> \param rd ...
349 : !> \param npgfa ...
350 : !> \param npgfb ...
351 : !> \param npgfc ...
352 : !> \param npgfd ...
353 : !> \param la_min ...
354 : !> \param la_max ...
355 : !> \param lb_min ...
356 : !> \param lb_max ...
357 : !> \param lc_min ...
358 : !> \param lc_max ...
359 : !> \param ld_min ...
360 : !> \param ld_max ...
361 : !> \param nsgfa ...
362 : !> \param nsgfb ...
363 : !> \param nsgfc ...
364 : !> \param nsgfd ...
365 : !> \param sphi_a_u1 ...
366 : !> \param sphi_a_u2 ...
367 : !> \param sphi_a_u3 ...
368 : !> \param sphi_b_u1 ...
369 : !> \param sphi_b_u2 ...
370 : !> \param sphi_b_u3 ...
371 : !> \param sphi_c_u1 ...
372 : !> \param sphi_c_u2 ...
373 : !> \param sphi_c_u3 ...
374 : !> \param sphi_d_u1 ...
375 : !> \param sphi_d_u2 ...
376 : !> \param sphi_d_u3 ...
377 : !> \param zeta ...
378 : !> \param zetb ...
379 : !> \param zetc ...
380 : !> \param zetd ...
381 : !> \param primitive_integrals array of primitive_integrals
382 : !> \param potential_parameter contains info for libint
383 : !> \param neighbor_cells Periodic images
384 : !> \param screen1 set based coefficients for near field screening
385 : !> \param screen2 set based coefficients for near field screening
386 : !> \param eps_schwarz threshold
387 : !> \param max_contraction_val maximum multiplication factor for cart -> sph
388 : !> \param cart_estimate maximum calculated integral value
389 : !> \param cell cell
390 : !> \param neris_tmp counter for calculated cart integrals
391 : !> \param log10_pmax logarithm of initial p matrix max element
392 : !> \param log10_eps_schwarz log of threshold
393 : !> \param R1_pgf coefficients for radii of product distribution function
394 : !> \param R2_pgf coefficients for radii of product distribution function
395 : !> \param pgf1 schwarz coefficients pgf basid
396 : !> \param pgf2 schwarz coefficients pgf basid
397 : !> \param pgf_list_ij ...
398 : !> \param pgf_list_kl ...
399 : !> \param pgf_product_list ...
400 : !> \param nsgfl_a ...
401 : !> \param nsgfl_b ...
402 : !> \param nsgfl_c ...
403 : !> \param nsgfl_d ...
404 : !> \param sphi_a_ext ...
405 : !> \param sphi_b_ext ...
406 : !> \param sphi_c_ext ...
407 : !> \param sphi_d_ext ...
408 : !> \param ee_work ...
409 : !> \param ee_work2 ...
410 : !> \param ee_buffer1 ...
411 : !> \param ee_buffer2 ...
412 : !> \param ee_primitives_tmp ...
413 : !> \param nimages ...
414 : !> \param do_periodic ...
415 : !> \par History
416 : !> 11.2006 created [Manuel Guidon]
417 : !> 02.2009 completely rewritten screening part [Manuel Guidon]
418 : !> \author Manuel Guidon
419 : ! **************************************************************************************************
420 8403174 : SUBROUTINE coulomb4(context, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, &
421 : la_min, la_max, lb_min, lb_max, &
422 : lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, &
423 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
424 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
425 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
426 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
427 8403174 : zeta, zetb, zetc, zetd, &
428 8403174 : primitive_integrals, &
429 : potential_parameter, neighbor_cells, &
430 : screen1, screen2, eps_schwarz, max_contraction_val, &
431 : cart_estimate, cell, neris_tmp, log10_pmax, &
432 : log10_eps_schwarz, R1_pgf, R2_pgf, pgf1, pgf2, &
433 : pgf_list_ij, pgf_list_kl, &
434 : pgf_product_list, &
435 8403174 : nsgfl_a, nsgfl_b, nsgfl_c, &
436 8403174 : nsgfl_d, &
437 8403174 : sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext, &
438 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
439 : nimages, do_periodic)
440 :
441 : TYPE(four_center_integral_context_type) :: context
442 : REAL(dp), INTENT(IN) :: ra(3), rb(3), rc(3), rd(3)
443 : INTEGER, INTENT(IN) :: npgfa, npgfb, npgfc, npgfd, la_min, la_max, lb_min, lb_max, lc_min, &
444 : lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, &
445 : sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, &
446 : sphi_d_u3
447 : REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
448 : REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
449 : REAL(dp), DIMENSION(1:npgfc), INTENT(IN) :: zetc
450 : REAL(dp), DIMENSION(1:npgfd), INTENT(IN) :: zetd
451 : REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitive_integrals
452 : TYPE(coulomb_operator_type) :: potential_parameter
453 : TYPE(eri4_cell_type), DIMENSION(:), POINTER :: neighbor_cells
454 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2), eps_schwarz, &
455 : max_contraction_val
456 : REAL(dp) :: cart_estimate
457 : TYPE(cell_type), POINTER :: cell
458 : INTEGER(int_8) :: neris_tmp
459 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
460 : TYPE(eri4_screen_coeff_type), DIMENSION(:, :), &
461 : POINTER :: R1_pgf, R2_pgf, pgf1, pgf2
462 : TYPE(eri4_pgf_pair_type), DIMENSION(*) :: pgf_list_ij, pgf_list_kl
463 : TYPE(eri4_pgf_product_type), ALLOCATABLE, &
464 : DIMENSION(:), INTENT(INOUT) :: pgf_product_list
465 : INTEGER, DIMENSION(0:), INTENT(IN) :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d
466 : REAL(dp), INTENT(IN) :: sphi_a_ext(sphi_a_u1, sphi_a_u2, sphi_a_u3), &
467 : sphi_b_ext(sphi_b_u1, sphi_b_u2, sphi_b_u3), sphi_c_ext(sphi_c_u1, sphi_c_u2, sphi_c_u3), &
468 : sphi_d_ext(sphi_d_u1, sphi_d_u2, sphi_d_u3)
469 : REAL(dp), DIMENSION(*) :: ee_work, ee_work2, ee_buffer1, &
470 : ee_buffer2, ee_primitives_tmp
471 : INTEGER, DIMENSION(*) :: nimages
472 : LOGICAL, INTENT(IN) :: do_periodic
473 :
474 : INTEGER :: ipgf, jpgf, kpgf, la, lb, lc, ld, list_ij, list_kl, lpgf, max_l, ncoa, ncob, &
475 : ncoc, ncod, nelements_ij, nelements_kl, nproducts, nsgfla, nsgflb, nsgflc, nsgfld, nsoa, &
476 : nsob, nsoc, nsod, s_offset_a, s_offset_b, s_offset_c, s_offset_d
477 : REAL(dp) :: tmp_max
478 :
479 8403174 : cart_estimate = 0.0_dp
480 8403174 : neris_tmp = 0
481 562108655 : primitive_integrals = 0.0_dp
482 :
483 : CALL build_pair_list_pgf(npgfa, npgfb, pgf_list_ij, zeta, zetb, screen1, screen2, &
484 : pgf1, R1_pgf, log10_pmax, log10_eps_schwarz, ra, rb, &
485 : nelements_ij, &
486 8403174 : neighbor_cells, nimages, do_periodic, cell)
487 8403174 : IF (nelements_ij == 0) RETURN
488 : CALL build_pair_list_pgf(npgfc, npgfd, pgf_list_kl, zetc, zetd, screen2, screen1, &
489 : pgf2, R2_pgf, log10_pmax, log10_eps_schwarz, rc, rd, &
490 : nelements_kl, &
491 8403174 : neighbor_cells, nimages, do_periodic, cell)
492 8403174 : IF (nelements_kl == 0) RETURN
493 :
494 8403174 : max_l = la_max + lb_max + lc_max + ld_max
495 :
496 28912185 : DO list_ij = 1, nelements_ij
497 20509011 : ipgf = pgf_list_ij(list_ij)%ipgf
498 20509011 : jpgf = pgf_list_ij(list_ij)%jpgf
499 :
500 104611152 : DO list_kl = 1, nelements_kl
501 75698967 : kpgf = pgf_list_kl(list_kl)%ipgf
502 75698967 : lpgf = pgf_list_kl(list_kl)%jpgf
503 :
504 : CALL build_pgf_product_list(pgf_list_ij(list_ij), pgf_list_kl(list_kl), pgf_product_list, &
505 : nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, &
506 75698967 : potential_parameter, max_l, do_periodic)
507 75698967 : IF (nproducts == 0) CYCLE
508 :
509 68082275 : s_offset_a = 0
510 245659153 : DO la = la_min, la_max
511 81368900 : s_offset_b = 0
512 81368900 : ncoa = nco(la)
513 81368900 : nsgfla = nsgfl_a(la)
514 81368900 : nsoa = nso(la)
515 :
516 174156519 : DO lb = lb_min, lb_max
517 92787619 : s_offset_c = 0
518 92787619 : ncob = nco(lb)
519 92787619 : nsgflb = nsgfl_b(lb)
520 92787619 : nsob = nso(lb)
521 :
522 219723419 : DO lc = lc_min, lc_max
523 126935800 : s_offset_d = 0
524 126935800 : ncoc = nco(lc)
525 126935800 : nsgflc = nsgfl_c(lc)
526 126935800 : nsoc = nso(lc)
527 :
528 297381939 : DO ld = ld_min, ld_max
529 170446139 : ncod = nco(ld)
530 170446139 : nsgfld = nsgfl_d(ld)
531 170446139 : nsod = nso(ld)
532 :
533 170446139 : tmp_max = 0.0_dp
534 : CALL evaluate_eri(context, nproducts, pgf_product_list, &
535 : la, lb, lc, ld, &
536 : ncoa, ncob, ncoc, ncod, &
537 : nsgfa, nsgfb, nsgfc, nsgfd, &
538 : primitive_integrals, &
539 : max_contraction_val, tmp_max, eps_schwarz, &
540 : neris_tmp, &
541 : s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
542 : nsgfla, nsgflb, nsgflc, nsgfld, nsoa, nsob, nsoc, nsod, &
543 : sphi_a_ext(1, la + 1, ipgf), &
544 : sphi_b_ext(1, lb + 1, jpgf), &
545 : sphi_c_ext(1, lc + 1, kpgf), &
546 : sphi_d_ext(1, ld + 1, lpgf), &
547 170446139 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp)
548 170446139 : cart_estimate = MAX(tmp_max, cart_estimate)
549 297381939 : s_offset_d = s_offset_d + nsod*nsgfld
550 : END DO !ld
551 219723419 : s_offset_c = s_offset_c + nsoc*nsgflc
552 : END DO !lc
553 174156519 : s_offset_b = s_offset_b + nsob*nsgflb
554 : END DO !lb
555 157067867 : s_offset_a = s_offset_a + nsoa*nsgfla
556 : END DO !la
557 : END DO
558 : END DO
559 :
560 : END SUBROUTINE coulomb4
561 :
562 : ! **************************************************************************************************
563 : !> \brief ...
564 : !> \param context four-center integral context
565 : !> \param ra ...
566 : !> \param rb ...
567 : !> \param rc ...
568 : !> \param rd ...
569 : !> \param npgfa ...
570 : !> \param npgfb ...
571 : !> \param npgfc ...
572 : !> \param npgfd ...
573 : !> \param la_min ...
574 : !> \param la_max ...
575 : !> \param lb_min ...
576 : !> \param lb_max ...
577 : !> \param lc_min ...
578 : !> \param lc_max ...
579 : !> \param ld_min ...
580 : !> \param ld_max ...
581 : !> \param nsgfa ...
582 : !> \param nsgfb ...
583 : !> \param nsgfc ...
584 : !> \param nsgfd ...
585 : !> \param sphi_a_u1 ...
586 : !> \param sphi_a_u2 ...
587 : !> \param sphi_a_u3 ...
588 : !> \param sphi_b_u1 ...
589 : !> \param sphi_b_u2 ...
590 : !> \param sphi_b_u3 ...
591 : !> \param sphi_c_u1 ...
592 : !> \param sphi_c_u2 ...
593 : !> \param sphi_c_u3 ...
594 : !> \param sphi_d_u1 ...
595 : !> \param sphi_d_u2 ...
596 : !> \param sphi_d_u3 ...
597 : !> \param zeta ...
598 : !> \param zetb ...
599 : !> \param zetc ...
600 : !> \param zetd ...
601 : !> \param primitive_integrals ...
602 : !> \param potential_parameter ...
603 : !> \param neighbor_cells ...
604 : !> \param screen1 ...
605 : !> \param screen2 ...
606 : !> \param eps_schwarz ...
607 : !> \param max_contraction_val ...
608 : !> \param cart_estimate ...
609 : !> \param cell ...
610 : !> \param neris_tmp ...
611 : !> \param log10_pmax ...
612 : !> \param log10_eps_schwarz ...
613 : !> \param R1_pgf ...
614 : !> \param R2_pgf ...
615 : !> \param pgf1 ...
616 : !> \param pgf2 ...
617 : !> \param pgf_list_ij ...
618 : !> \param pgf_list_kl ...
619 : !> \param pgf_product_list ...
620 : !> \param nsgfl_a ...
621 : !> \param nsgfl_b ...
622 : !> \param nsgfl_c ...
623 : !> \param nsgfl_d ...
624 : !> \param sphi_a_ext ...
625 : !> \param sphi_b_ext ...
626 : !> \param sphi_c_ext ...
627 : !> \param sphi_d_ext ...
628 : !> \param ede_work ...
629 : !> \param ede_work2 ...
630 : !> \param ede_work_forces ...
631 : !> \param ede_buffer1 ...
632 : !> \param ede_buffer2 ...
633 : !> \param ede_primitives_tmp ...
634 : !> \param nimages ...
635 : !> \param do_periodic ...
636 : !> \param use_virial ...
637 : !> \param ede_work_virial ...
638 : !> \param ede_work2_virial ...
639 : !> \param ede_primitives_tmp_virial ...
640 : !> \param primitive_integrals_virial ...
641 : !> \param cart_estimate_virial ...
642 : ! **************************************************************************************************
643 1841077 : SUBROUTINE forces4(context, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, &
644 : la_min, la_max, lb_min, lb_max, &
645 : lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, &
646 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
647 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
648 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
649 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
650 1841077 : zeta, zetb, zetc, zetd, &
651 1841077 : primitive_integrals, &
652 : potential_parameter, neighbor_cells, &
653 : screen1, screen2, eps_schwarz, max_contraction_val, &
654 : cart_estimate, cell, neris_tmp, log10_pmax, &
655 : log10_eps_schwarz, R1_pgf, R2_pgf, pgf1, pgf2, &
656 : pgf_list_ij, pgf_list_kl, &
657 : pgf_product_list, &
658 1841077 : nsgfl_a, nsgfl_b, nsgfl_c, &
659 1841077 : nsgfl_d, &
660 1841077 : sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext, &
661 : ede_work, ede_work2, ede_work_forces, &
662 : ede_buffer1, ede_buffer2, ede_primitives_tmp, &
663 : nimages, do_periodic, use_virial, ede_work_virial, &
664 1841077 : ede_work2_virial, ede_primitives_tmp_virial, primitive_integrals_virial, &
665 : cart_estimate_virial)
666 :
667 : TYPE(four_center_integral_context_type) :: context
668 : REAL(dp), INTENT(IN) :: ra(3), rb(3), rc(3), rd(3)
669 : INTEGER, INTENT(IN) :: npgfa, npgfb, npgfc, npgfd, la_min, la_max, lb_min, lb_max, lc_min, &
670 : lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, &
671 : sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, &
672 : sphi_d_u3
673 : REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
674 : REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
675 : REAL(dp), DIMENSION(1:npgfc), INTENT(IN) :: zetc
676 : REAL(dp), DIMENSION(1:npgfd), INTENT(IN) :: zetd
677 : REAL(dp), &
678 : DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12) :: primitive_integrals
679 : TYPE(coulomb_operator_type) :: potential_parameter
680 : TYPE(eri4_cell_type), DIMENSION(:), POINTER :: neighbor_cells
681 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2), eps_schwarz, &
682 : max_contraction_val
683 : REAL(dp) :: cart_estimate
684 : TYPE(cell_type), POINTER :: cell
685 : INTEGER(int_8) :: neris_tmp
686 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
687 : TYPE(eri4_screen_coeff_type), DIMENSION(:, :), &
688 : POINTER :: R1_pgf, R2_pgf, pgf1, pgf2
689 : TYPE(eri4_pgf_pair_type), DIMENSION(*) :: pgf_list_ij, pgf_list_kl
690 : TYPE(eri4_pgf_product_type), ALLOCATABLE, &
691 : DIMENSION(:), INTENT(INOUT) :: pgf_product_list
692 : INTEGER, DIMENSION(0:), INTENT(IN) :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d
693 : REAL(dp), INTENT(IN) :: sphi_a_ext(sphi_a_u1, sphi_a_u2, sphi_a_u3), &
694 : sphi_b_ext(sphi_b_u1, sphi_b_u2, sphi_b_u3), sphi_c_ext(sphi_c_u1, sphi_c_u2, sphi_c_u3), &
695 : sphi_d_ext(sphi_d_u1, sphi_d_u2, sphi_d_u3)
696 : REAL(dp), DIMENSION(*) :: ede_work, ede_work2, ede_work_forces, &
697 : ede_buffer1, ede_buffer2, &
698 : ede_primitives_tmp
699 : INTEGER, DIMENSION(*) :: nimages
700 : LOGICAL, INTENT(IN) :: do_periodic, use_virial
701 : REAL(dp), DIMENSION(*) :: ede_work_virial, ede_work2_virial, &
702 : ede_primitives_tmp_virial
703 : REAL(dp), &
704 : DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12, 3) :: primitive_integrals_virial
705 : REAL(dp) :: cart_estimate_virial
706 :
707 : INTEGER :: ipgf, jpgf, kpgf, la, lb, lc, ld, list_ij, list_kl, lpgf, max_l, ncoa, ncob, &
708 : ncoc, ncod, nelements_ij, nelements_kl, nproducts, nsgfla, nsgflb, nsgflc, nsgfld, nsoa, &
709 : nsob, nsoc, nsod, s_offset_a, s_offset_b, s_offset_c, s_offset_d
710 : REAL(dp) :: tmp_max, tmp_max_virial
711 :
712 1841077 : cart_estimate = 0.0_dp
713 978203425 : primitive_integrals = 0.0_dp
714 1841077 : IF (use_virial) THEN
715 288753616 : primitive_integrals_virial = 0.0_dp
716 203542 : cart_estimate_virial = 0.0_dp
717 : END IF
718 1841077 : neris_tmp = 0_int_8
719 :
720 : CALL build_pair_list_pgf(npgfa, npgfb, pgf_list_ij, zeta, zetb, screen1, screen2, &
721 : pgf1, R1_pgf, log10_pmax, log10_eps_schwarz, ra, rb, &
722 : nelements_ij, &
723 1841077 : neighbor_cells, nimages, do_periodic, cell)
724 1841077 : IF (nelements_ij == 0) RETURN
725 : CALL build_pair_list_pgf(npgfc, npgfd, pgf_list_kl, zetc, zetd, screen2, screen1, &
726 : pgf2, R2_pgf, log10_pmax, log10_eps_schwarz, rc, rd, &
727 : nelements_kl, &
728 1841077 : neighbor_cells, nimages, do_periodic, cell)
729 1841077 : IF (nelements_kl == 0) RETURN
730 :
731 1841077 : max_l = la_max + lb_max + lc_max + ld_max + 1
732 5048161 : DO list_ij = 1, nelements_ij
733 3207084 : ipgf = pgf_list_ij(list_ij)%ipgf
734 3207084 : jpgf = pgf_list_ij(list_ij)%jpgf
735 :
736 19897306 : DO list_kl = 1, nelements_kl
737 14849145 : kpgf = pgf_list_kl(list_kl)%ipgf
738 14849145 : lpgf = pgf_list_kl(list_kl)%jpgf
739 :
740 : CALL build_pgf_product_list(pgf_list_ij(list_ij), pgf_list_kl(list_kl), pgf_product_list, &
741 : nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, &
742 14849145 : potential_parameter, max_l, do_periodic)
743 14849145 : IF (nproducts == 0) CYCLE
744 14092867 : s_offset_a = 0
745 49161209 : DO la = la_min, la_max
746 17012113 : s_offset_b = 0
747 17012113 : ncoa = nco(la)
748 17012113 : nsgfla = nsgfl_a(la)
749 17012113 : nsoa = nso(la)
750 :
751 35137194 : DO lb = lb_min, lb_max
752 18125081 : s_offset_c = 0
753 18125081 : ncob = nco(lb)
754 18125081 : nsgflb = nsgfl_b(lb)
755 18125081 : nsob = nso(lb)
756 :
757 44642483 : DO lc = lc_min, lc_max
758 26517402 : s_offset_d = 0
759 26517402 : ncoc = nco(lc)
760 26517402 : nsgflc = nsgfl_c(lc)
761 26517402 : nsoc = nso(lc)
762 :
763 60625150 : DO ld = ld_min, ld_max
764 34107748 : ncod = nco(ld)
765 34107748 : nsgfld = nsgfl_d(ld)
766 34107748 : nsod = nso(ld)
767 34107748 : tmp_max = 0.0_dp
768 34107748 : tmp_max_virial = 0.0_dp
769 : CALL evaluate_deriv_eri(context, nproducts, pgf_product_list, &
770 : la, lb, lc, ld, &
771 : ncoa, ncob, ncoc, ncod, &
772 : nsgfa, nsgfb, nsgfc, nsgfd, &
773 : primitive_integrals, &
774 : max_contraction_val, tmp_max, eps_schwarz, neris_tmp, &
775 : s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
776 : nsgfla, nsgflb, nsgflc, nsgfld, nsoa, nsob, nsoc, nsod, &
777 : sphi_a_ext(1, la + 1, ipgf), &
778 : sphi_b_ext(1, lb + 1, jpgf), &
779 : sphi_c_ext(1, lc + 1, kpgf), &
780 : sphi_d_ext(1, ld + 1, lpgf), &
781 : ede_work, ede_work2, ede_work_forces, &
782 : ede_buffer1, ede_buffer2, ede_primitives_tmp, use_virial, &
783 : ede_work_virial, ede_work2_virial, ede_primitives_tmp_virial, &
784 34107748 : primitive_integrals_virial, cell, tmp_max_virial)
785 34107748 : cart_estimate = MAX(tmp_max, cart_estimate)
786 34107748 : IF (use_virial) cart_estimate_virial = MAX(tmp_max_virial, cart_estimate_virial)
787 60625150 : s_offset_d = s_offset_d + nsod*nsgfld
788 : END DO !ld
789 44642483 : s_offset_c = s_offset_c + nsoc*nsgflc
790 : END DO !lc
791 35137194 : s_offset_b = s_offset_b + nsob*nsgflb
792 : END DO !lb
793 31861258 : s_offset_a = s_offset_a + nsoa*nsgfla
794 : END DO !la
795 : END DO
796 : END DO
797 :
798 : END SUBROUTINE forces4
799 :
800 : END MODULE four_center_integral_methods
|
