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 Four-center integral evaluation with an encapsulated Libint backend.
10 : !> \par History
11 : !> 11.2006 created [Manuel Guidon]
12 : !> 11.2019 Fixed potential_id initial values (A. Bussy)
13 : !> \author Manuel Guidon
14 : ! **************************************************************************************************
15 : MODULE four_center_integral_interface
16 :
17 : USE cell_types, ONLY: cell_type,&
18 : real_to_scaled
19 : USE four_center_integral_preparation,ONLY: prepare_eri4_pair,&
20 : prepare_eri4_primitive,&
21 : prepare_eri4_product
22 : USE four_center_integral_types, ONLY: eri4_pgf_image_type,&
23 : eri4_pgf_pair_type,&
24 : eri4_pgf_product_type
25 : USE hfx_contraction_methods, ONLY: contract
26 : USE input_constants, ONLY: 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 libint_wrapper, ONLY: &
33 : cp_libint_cleanup_eri, cp_libint_cleanup_eri1, cp_libint_get_derivs, cp_libint_get_eris, &
34 : cp_libint_init_eri, cp_libint_init_eri1, cp_libint_set_contrdepth, &
35 : cp_libint_set_params_eri, cp_libint_set_params_eri_deriv, cp_libint_t
36 : USE orbital_pointers, ONLY: nco
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 : PUBLIC :: four_center_integral_init, four_center_integral_cleanup, &
42 : evaluate_eri, &
43 : evaluate_deriv_eri, &
44 : evaluate_eri_screen
45 :
46 : !> \brief Context owning the energy and first-derivative backend workspaces.
47 : !> \note Intrinsic copies share backend storage and must not be cleaned up independently.
48 : TYPE, PUBLIC :: four_center_integral_context_type
49 : PRIVATE
50 : TYPE(cp_libint_t) :: eri, deriv
51 : END TYPE four_center_integral_context_type
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Initialize the four-center energy and first-derivative workspaces.
57 : !> \param context owning context; global backend initialization is provided by init_cp2k
58 : !> \param max_am maximum angular momentum on any center
59 : ! **************************************************************************************************
60 1622 : SUBROUTINE four_center_integral_init(context, max_am)
61 : TYPE(four_center_integral_context_type), &
62 : INTENT(INOUT) :: context
63 : INTEGER, INTENT(IN) :: max_am
64 :
65 1622 : CALL cp_libint_init_eri(context%eri, max_am)
66 1622 : CALL cp_libint_init_eri1(context%deriv, max_am)
67 1622 : CALL cp_libint_set_contrdepth(context%eri, 1)
68 1622 : CALL cp_libint_set_contrdepth(context%deriv, 1)
69 1622 : END SUBROUTINE four_center_integral_init
70 :
71 : ! **************************************************************************************************
72 : !> \brief Release the workspaces of an owning four-center context.
73 : !> \param context original initialized owner; borrowed copies must not be passed here
74 : ! **************************************************************************************************
75 1622 : SUBROUTINE four_center_integral_cleanup(context)
76 : TYPE(four_center_integral_context_type), &
77 : INTENT(INOUT) :: context
78 :
79 1622 : CALL cp_libint_cleanup_eri(context%eri)
80 1622 : CALL cp_libint_cleanup_eri1(context%deriv)
81 1622 : END SUBROUTINE four_center_integral_cleanup
82 :
83 : ! **************************************************************************************************
84 : !> \brief Return the canonical quartet order required by the Libint ERI kernels.
85 : !> \param la angular momentum of the original first center
86 : !> \param lb angular momentum of the original second center
87 : !> \param lc angular momentum of the original third center
88 : !> \param ld angular momentum of the original fourth center
89 : !> \param permutation original center index for each canonical center
90 : ! **************************************************************************************************
91 265482746 : PURE SUBROUTINE get_eri_permutation(la, lb, lc, ld, permutation)
92 : INTEGER, INTENT(IN) :: la, lb, lc, ld
93 : INTEGER, INTENT(OUT) :: permutation(4)
94 :
95 265482746 : permutation = [1, 2, 3, 4]
96 443863121 : IF (la < lb) permutation(1:2) = [2, 1]
97 471979235 : IF (lc < ld) permutation(3:4) = [4, 3]
98 795628946 : IF (la + lb > lc + ld) permutation = permutation([3, 4, 1, 2])
99 265482746 : END SUBROUTINE get_eri_permutation
100 :
101 : ! **************************************************************************************************
102 : !> \brief Restore original Cartesian axes and center labels from a canonical ERI block.
103 : !> \param source canonical block, with the fourth Cartesian index varying fastest
104 : !> \param destination original block, with the first Cartesian index varying fastest
105 : !> \param ncart Cartesian shell sizes in original center order
106 : !> \param permutation original center index for each canonical center
107 : !> \param ncomponents one integral component or twelve center-derivative components
108 : !> \param nplanes independent planes, one for integrals/forces or three for virial
109 : ! **************************************************************************************************
110 1285574733 : PURE SUBROUTINE reorder_eri_block(source, destination, ncart, permutation, ncomponents, nplanes)
111 : INTEGER, INTENT(IN) :: ncart(4), permutation(4), ncomponents, &
112 : nplanes
113 : REAL(dp), INTENT(OUT) :: destination(PRODUCT(ncart), ncomponents, nplanes)
114 : REAL(dp), INTENT(IN) :: source(PRODUCT(ncart), ncomponents, nplanes)
115 :
116 : INTEGER :: component, extent(4), i, isource, j, k, &
117 : l, output_component, p1, p2, p3, p4, &
118 : plane, stride(4)
119 :
120 714208185 : stride = [1, ncart(1), ncart(1)*ncart(2), ncart(1)*ncart(2)*ncart(3)]
121 1428416370 : stride = stride(permutation)
122 714208185 : extent = ncart(permutation)
123 : ! Read the Libint block contiguously and restore the original Fortran layout with strides.
124 287847856 : DO plane = 1, nplanes
125 843754906 : DO component = 1, ncomponents
126 555907050 : output_component = component
127 555907050 : IF (ncomponents == 12) THEN
128 448255452 : output_component = 3*(permutation((component - 1)/3 + 1) - 1) + MOD(component - 1, 3) + 1
129 : END IF
130 555907050 : isource = 0
131 555907050 : p1 = 1
132 1657115065 : DO i = 1, extent(1)
133 956201796 : p2 = p1
134 2034332538 : DO j = 1, extent(2)
135 1078130742 : p3 = p2
136 4971677542 : DO k = 1, extent(3)
137 3893546800 : p4 = p3
138 12055061046 : DO l = 1, extent(4)
139 8161514246 : isource = isource + 1
140 8161514246 : destination(p4, output_component, plane) = source(isource, component, plane)
141 12055061046 : p4 = p4 + stride(4)
142 : END DO
143 4971677542 : p3 = p3 + stride(3)
144 : END DO
145 2034332538 : p2 = p2 + stride(2)
146 : END DO
147 1512108846 : p1 = p1 + stride(1)
148 : END DO
149 : END DO
150 : END DO
151 142841637 : END SUBROUTINE reorder_eri_block
152 :
153 : ! **************************************************************************************************
154 : !> \brief Evaluates derivatives of electron repulsion integrals for a primitive quartet
155 : !> \param context four-center integral context
156 : !> \param nproducts ...
157 : !> \param pgf_product_list ...
158 : !> \param n_a ...
159 : !> \param n_b ...
160 : !> \param n_c ...
161 : !> \param n_d ...
162 : !> \param ncoa ...
163 : !> \param ncob ...
164 : !> \param ncoc ...
165 : !> \param ncod ...
166 : !> \param nsgfa ...
167 : !> \param nsgfb ...
168 : !> \param nsgfc ...
169 : !> \param nsgfd ...
170 : !> \param primitives ...
171 : !> \param max_contraction ...
172 : !> \param tmp_max_all ...
173 : !> \param eps_schwarz ...
174 : !> \param neris ...
175 : !> \param s_offset_a ...
176 : !> \param s_offset_b ...
177 : !> \param s_offset_c ...
178 : !> \param s_offset_d ...
179 : !> \param nl_a ...
180 : !> \param nl_b ...
181 : !> \param nl_c ...
182 : !> \param nl_d ...
183 : !> \param nsoa ...
184 : !> \param nsob ...
185 : !> \param nsoc ...
186 : !> \param nsod ...
187 : !> \param sphi_a ...
188 : !> \param sphi_b ...
189 : !> \param sphi_c ...
190 : !> \param sphi_d ...
191 : !> \param work ...
192 : !> \param work2 ...
193 : !> \param work_forces ...
194 : !> \param buffer1 ...
195 : !> \param buffer2 ...
196 : !> \param primitives_tmp ...
197 : !> \param use_virial ...
198 : !> \param work_virial ...
199 : !> \param work2_virial ...
200 : !> \param primitives_tmp_virial ...
201 : !> \param primitives_virial ...
202 : !> \param cell ...
203 : !> \param tmp_max_all_virial ...
204 : !> \par History
205 : !> 03.2007 created [Manuel Guidon]
206 : !> 08.2007 refactured permutation part [Manuel Guidon]
207 : !> \author Manuel Guidon
208 : ! **************************************************************************************************
209 34107748 : SUBROUTINE evaluate_deriv_eri(context, nproducts, pgf_product_list, &
210 : n_a, n_b, n_c, n_d, &
211 : ncoa, ncob, ncoc, ncod, &
212 : nsgfa, nsgfb, nsgfc, nsgfd, &
213 34107748 : primitives, max_contraction, tmp_max_all, &
214 : eps_schwarz, neris, &
215 : s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
216 : nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, &
217 34107748 : sphi_a, sphi_b, sphi_c, sphi_d, &
218 34107748 : work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
219 34107748 : use_virial, work_virial, work2_virial, primitives_tmp_virial, &
220 34107748 : primitives_virial, cell, tmp_max_all_virial)
221 :
222 : TYPE(four_center_integral_context_type) :: context
223 : INTEGER, INTENT(IN) :: nproducts
224 : TYPE(eri4_pgf_product_type), DIMENSION(nproducts) :: pgf_product_list
225 : INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, &
226 : ncod, nsgfa, nsgfb, nsgfc, nsgfd
227 : REAL(dp), &
228 : DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12) :: primitives
229 : REAL(dp) :: max_contraction, tmp_max_all, eps_schwarz
230 : INTEGER(int_8) :: neris
231 : INTEGER :: s_offset_a, s_offset_b, s_offset_c, &
232 : s_offset_d, nl_a, nl_b, nl_c, nl_d, &
233 : nsoa, nsob, nsoc, nsod
234 : REAL(dp), DIMENSION(ncoa, nsoa*nl_a) :: sphi_a
235 : REAL(dp), DIMENSION(ncob, nsob*nl_b) :: sphi_b
236 : REAL(dp), DIMENSION(ncoc, nsoc*nl_c) :: sphi_c
237 : REAL(dp), DIMENSION(ncod, nsod*nl_d) :: sphi_d
238 : REAL(dp) :: work(ncoa*ncob*ncoc*ncod, 12), work2(ncoa, ncob, ncoc, ncod, 12), &
239 : work_forces(ncoa*ncob*ncoc*ncod, 12)
240 : REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2
241 : REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*&
242 : nl_c, nsod*nl_d) :: primitives_tmp
243 : LOGICAL, INTENT(IN) :: use_virial
244 : REAL(dp) :: work_virial(ncoa*ncob*ncoc*ncod, 12, 3), &
245 : work2_virial(ncoa, ncob, ncoc, ncod, 12, 3)
246 : REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*&
247 : nl_c, nsod*nl_d) :: primitives_tmp_virial
248 : REAL(dp), &
249 : DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12, 3) :: primitives_virial
250 : TYPE(cell_type), POINTER :: cell
251 : REAL(dp) :: tmp_max_all_virial
252 :
253 : INTEGER :: a_mysize(1), angular_momenta(4), i, j, &
254 : k, m, m_max, mysize, n, permutation(4)
255 : LOGICAL :: direct_contract
256 : REAL(dp) :: scoord(12), tmp_max, tmp_max_virial
257 :
258 34107748 : m_max = n_a + n_b + n_c + n_d
259 34107748 : m_max = m_max + 1
260 34107748 : mysize = ncoa*ncob*ncoc*ncod
261 68215496 : a_mysize = mysize
262 :
263 4881735712 : work = 0.0_dp
264 :
265 34107748 : IF (use_virial) THEN
266 858458504 : work_virial = 0.0_dp
267 : END IF
268 :
269 34107748 : CALL get_eri_permutation(n_a, n_b, n_c, n_d, permutation)
270 170538740 : angular_momenta = [n_a, n_b, n_c, n_d]
271 306969732 : angular_momenta = angular_momenta(permutation)
272 :
273 140764572 : DO i = 1, nproducts
274 106656824 : CALL cp_libint_set_params_eri_deriv(context%deriv, pgf_product_list(i), permutation, m_max)
275 : CALL cp_libint_get_derivs(angular_momenta(4), angular_momenta(3), angular_momenta(2), &
276 106656824 : angular_momenta(1), context%deriv, work_forces, a_mysize)
277 1386538712 : DO k = 1, 12
278 15368578088 : DO j = 1, mysize
279 15261921264 : work(j, k) = work(j, k) + work_forces(j, k)
280 : END DO
281 : END DO
282 106656824 : neris = neris + 12*mysize
283 140764572 : IF (use_virial) THEN
284 103732160 : DO j = 1, 4
285 103732160 : CALL real_to_scaled(scoord(3*j - 2:3*j), pgf_product_list(i)%centers(:, permutation(j)), cell)
286 : END DO
287 269703616 : DO k = 1, 12
288 2501243008 : DO j = 1, mysize
289 9175114752 : DO m = 1, 3
290 8926157568 : work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
291 : END DO
292 : END DO
293 : END DO
294 : END IF
295 : END DO
296 :
297 443400724 : DO n = 1, 12
298 : tmp_max = 0.0_dp
299 4847627964 : DO i = 1, mysize
300 4847627964 : tmp_max = MAX(tmp_max, ABS(work(i, n)))
301 : END DO
302 409292976 : tmp_max = tmp_max*max_contraction
303 443400724 : tmp_max_all = MAX(tmp_max_all, tmp_max)
304 : END DO
305 170538740 : CALL reorder_eri_block(work, work2, [ncoa, ncob, ncoc, ncod], permutation, 12, 1)
306 34107748 : IF (use_virial) THEN
307 14069783 : DO n = 1, 12
308 : tmp_max_virial = 0.0_dp
309 284709780 : DO i = 1, mysize
310 : tmp_max_virial = MAX(tmp_max_virial, &
311 284709780 : ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
312 : END DO
313 12987492 : tmp_max_virial = tmp_max_virial*max_contraction
314 14069783 : tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
315 : END DO
316 5411455 : CALL reorder_eri_block(work_virial, work2_virial, [ncoa, ncob, ncoc, ncod], permutation, 12, 3)
317 : END IF
318 :
319 34107748 : IF (.NOT. use_virial) THEN
320 33025457 : IF (tmp_max_all < eps_schwarz) RETURN
321 : END IF
322 :
323 : direct_contract = (s_offset_a == 0 .AND. s_offset_b == 0 .AND. &
324 : s_offset_c == 0 .AND. s_offset_d == 0 .AND. &
325 : nsgfa == nsoa*nl_a .AND. nsgfb == nsob*nl_b .AND. &
326 29503435 : nsgfc == nsoc*nl_c .AND. nsgfd == nsod*nl_d)
327 :
328 29503435 : IF (tmp_max_all >= eps_schwarz) THEN
329 29343593 : IF (direct_contract) THEN
330 75431876 : DO i = 1, 12
331 : CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
332 : n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2(1, 1, 1, 1, i), &
333 : sphi_a, &
334 : sphi_b, &
335 : sphi_c, &
336 : sphi_d, &
337 : primitives(1, 1, 1, 1, i), &
338 75431876 : buffer1, buffer2)
339 : END DO
340 : ELSE
341 306034833 : DO i = 1, 12
342 19213285128 : primitives_tmp(:, :, :, :) = 0.0_dp
343 : CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
344 : n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2(1, 1, 1, 1, i), &
345 : sphi_a, &
346 : sphi_b, &
347 : sphi_c, &
348 : sphi_d, &
349 : primitives_tmp(1, 1, 1, 1), &
350 282493692 : buffer1, buffer2)
351 :
352 : primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
353 : s_offset_b + 1:s_offset_b + nsob*nl_b, &
354 : s_offset_c + 1:s_offset_c + nsoc*nl_c, &
355 : s_offset_d + 1:s_offset_d + nsod*nl_d, i) = &
356 : primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
357 : s_offset_b + 1:s_offset_b + nsob*nl_b, &
358 : s_offset_c + 1:s_offset_c + nsoc*nl_c, &
359 19236826269 : s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
360 : END DO
361 : END IF
362 : END IF
363 :
364 29503435 : IF (use_virial .AND. tmp_max_all_virial >= eps_schwarz) THEN
365 935149 : IF (direct_contract) THEN
366 3274375 : DO i = 1, 12
367 12341875 : DO m = 1, 3
368 : CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
369 : n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2_virial(1, 1, 1, 1, i, m), &
370 : sphi_a, &
371 : sphi_b, &
372 : sphi_c, &
373 : sphi_d, &
374 : primitives_virial(1, 1, 1, 1, i, m), &
375 12090000 : buffer1, buffer2)
376 : END DO
377 : END DO
378 : ELSE
379 8882562 : DO i = 1, 12
380 33480426 : DO m = 1, 3
381 4741145568 : primitives_tmp_virial(:, :, :, :) = 0.0_dp
382 : CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
383 : n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2_virial(1, 1, 1, 1, i, m), &
384 : sphi_a, &
385 : sphi_b, &
386 : sphi_c, &
387 : sphi_d, &
388 : primitives_tmp_virial(1, 1, 1, 1), &
389 24597864 : buffer1, buffer2)
390 :
391 : primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
392 : s_offset_b + 1:s_offset_b + nsob*nl_b, &
393 : s_offset_c + 1:s_offset_c + nsoc*nl_c, &
394 : s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) = &
395 : primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
396 : s_offset_b + 1:s_offset_b + nsob*nl_b, &
397 : s_offset_c + 1:s_offset_c + nsoc*nl_c, &
398 4749344856 : s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) + primitives_tmp_virial(:, :, :, :)
399 : END DO
400 : END DO
401 : END IF
402 : END IF
403 :
404 : END SUBROUTINE evaluate_deriv_eri
405 :
406 : ! **************************************************************************************************
407 : !> \brief Evaluates max-abs values of electron repulsion integrals for a primitive quartet
408 : !> \param context four-center integral context
409 : !> \param A ...
410 : !> \param B ...
411 : !> \param C ...
412 : !> \param D ...
413 : !> \param Zeta_A ...
414 : !> \param Zeta_B ...
415 : !> \param Zeta_C ...
416 : !> \param Zeta_D ...
417 : !> \param n_a ...
418 : !> \param n_b ...
419 : !> \param n_c ...
420 : !> \param n_d ...
421 : !> \param max_val ...
422 : !> \param potential_parameter ...
423 : !> \param R1 ...
424 : !> \param R2 ...
425 : !> \par History
426 : !> 03.2007 created [Manuel Guidon]
427 : !> 08.2007 refactured permutation part [Manuel Guidon]
428 : !> \author Manuel Guidon
429 : ! **************************************************************************************************
430 96382482 : SUBROUTINE evaluate_eri_screen(context, A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, &
431 : n_a, n_b, n_c, n_d, &
432 : max_val, potential_parameter, R1, R2)
433 :
434 : TYPE(four_center_integral_context_type) :: context
435 : REAL(dp), INTENT(IN) :: A(3), B(3), C(3), D(3), Zeta_A, Zeta_B, &
436 : Zeta_C, Zeta_D
437 : INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d
438 : REAL(dp), INTENT(INOUT) :: max_val
439 : TYPE(coulomb_operator_type) :: potential_parameter
440 : REAL(dp) :: R1, R2
441 :
442 : INTEGER :: a_mysize(1), angular_momenta(4), i, &
443 : m_max, mysize, permutation(4)
444 : REAL(dp) :: centers(3, 4), exponents(4), PQ(3), PQ2, &
445 : S1234
446 96382482 : REAL(dp), DIMENSION(:), POINTER :: p_work
447 : TYPE(eri4_pgf_image_type) :: image_a, image_b
448 : TYPE(eri4_pgf_pair_type) :: pair_a, pair_b
449 : TYPE(eri4_pgf_product_type) :: params, product
450 :
451 96382482 : m_max = n_a + n_b + n_c + n_d
452 96382482 : mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
453 192764964 : a_mysize = mysize
454 :
455 96382482 : CALL get_eri_permutation(n_a, n_b, n_c, n_d, permutation)
456 481912410 : angular_momenta = [n_a, n_b, n_c, n_d]
457 867442338 : angular_momenta = angular_momenta(permutation)
458 385529928 : centers(:, 1) = A
459 385529928 : centers(:, 2) = B
460 385529928 : centers(:, 3) = C
461 385529928 : centers(:, 4) = D
462 3180621906 : centers = centers(:, permutation)
463 481912410 : exponents = [Zeta_A, Zeta_B, Zeta_C, Zeta_D]
464 867442338 : exponents = exponents(permutation)
465 :
466 385529928 : PQ = centers(:, 1) - centers(:, 2)
467 : CALL prepare_eri4_pair(centers(:, 1), centers(:, 2), exponents(1), exponents(2), &
468 385529928 : DOT_PRODUCT(PQ, PQ), pair_a, image_a)
469 385529928 : PQ = centers(:, 3) - centers(:, 4)
470 : CALL prepare_eri4_pair(centers(:, 3), centers(:, 4), exponents(3), exponents(4), &
471 385529928 : DOT_PRODUCT(PQ, PQ), pair_b, image_b)
472 385529928 : PQ = image_a%P - image_b%P
473 385529928 : PQ2 = DOT_PRODUCT(PQ, PQ)
474 113759532 : SELECT CASE (potential_parameter%potential_type)
475 : CASE (do_potential_truncated, do_potential_short, do_potential_mix_cl_trunc)
476 96382482 : IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) RETURN
477 : END SELECT
478 :
479 96382482 : CALL prepare_eri4_product(pair_a, pair_b, m_max, product)
480 96382482 : S1234 = EXP(image_a%S1234 + image_b%S1234)
481 96382482 : CALL prepare_eri4_primitive(product, image_a, image_b, S1234, PQ2, potential_parameter, params)
482 96382482 : IF (m_max /= 0) THEN
483 60006726 : CALL cp_libint_set_params_eri(context%eri, params, [1, 2, 3, 4], m_max)
484 : CALL cp_libint_get_eris(angular_momenta(4), angular_momenta(3), angular_momenta(2), &
485 60006726 : angular_momenta(1), context%eri, p_work, a_mysize)
486 4154086974 : DO i = 1, mysize
487 4154086974 : max_val = MAX(max_val, ABS(p_work(i)))
488 : END DO
489 : ELSE
490 36375756 : max_val = ABS(params%Fm(1))
491 : END IF
492 :
493 12047810250 : END SUBROUTINE evaluate_eri_screen
494 :
495 : ! **************************************************************************************************
496 : !> \brief Evaluate electron repulsion integrals for a primitive quartet
497 : !> \param context four-center integral context
498 : !> \param nproducts ...
499 : !> \param pgf_product_list ...
500 : !> \param n_a ...
501 : !> \param n_b ...
502 : !> \param n_c ...
503 : !> \param n_d ...
504 : !> \param ncoa ...
505 : !> \param ncob ...
506 : !> \param ncoc ...
507 : !> \param ncod ...
508 : !> \param nsgfa ...
509 : !> \param nsgfb ...
510 : !> \param nsgfc ...
511 : !> \param nsgfd ...
512 : !> \param primitives ...
513 : !> \param max_contraction ...
514 : !> \param tmp_max ...
515 : !> \param eps_schwarz ...
516 : !> \param neris ...
517 : !> \param s_offset_a ...
518 : !> \param s_offset_b ...
519 : !> \param s_offset_c ...
520 : !> \param s_offset_d ...
521 : !> \param nl_a ...
522 : !> \param nl_b ...
523 : !> \param nl_c ...
524 : !> \param nl_d ...
525 : !> \param nsoa ...
526 : !> \param nsob ...
527 : !> \param nsoc ...
528 : !> \param nsod ...
529 : !> \param sphi_a ...
530 : !> \param sphi_b ...
531 : !> \param sphi_c ...
532 : !> \param sphi_d ...
533 : !> \param work ...
534 : !> \param work2 ...
535 : !> \param buffer1 ...
536 : !> \param buffer2 ...
537 : !> \param primitives_tmp ...
538 : !> \par History
539 : !> 11.2006 created [Manuel Guidon]
540 : !> 08.2007 refactured permutation part [Manuel Guidon]
541 : !> \author Manuel Guidon
542 : ! **************************************************************************************************
543 170446139 : SUBROUTINE evaluate_eri(context, nproducts, pgf_product_list, &
544 : n_a, n_b, n_c, n_d, &
545 : ncoa, ncob, ncoc, ncod, &
546 : nsgfa, nsgfb, nsgfc, nsgfd, &
547 170446139 : primitives, max_contraction, tmp_max, &
548 : eps_schwarz, neris, &
549 : s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
550 : nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, &
551 170446139 : sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
552 170446139 : primitives_tmp)
553 :
554 : TYPE(four_center_integral_context_type) :: context
555 : INTEGER, INTENT(IN) :: nproducts
556 : TYPE(eri4_pgf_product_type), DIMENSION(nproducts) :: pgf_product_list
557 : INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, &
558 : ncod, nsgfa, nsgfb, nsgfc, nsgfd
559 : REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitives
560 : REAL(dp) :: max_contraction, tmp_max, eps_schwarz
561 : INTEGER(int_8) :: neris
562 : INTEGER :: s_offset_a, s_offset_b, s_offset_c, &
563 : s_offset_d, nl_a, nl_b, nl_c, nl_d, &
564 : nsoa, nsob, nsoc, nsod
565 : REAL(dp), DIMENSION(ncoa, nsoa*nl_a), INTENT(IN) :: sphi_a
566 : REAL(dp), DIMENSION(ncob, nsob*nl_b), INTENT(IN) :: sphi_b
567 : REAL(dp), DIMENSION(ncoc, nsoc*nl_c), INTENT(IN) :: sphi_c
568 : REAL(dp), DIMENSION(ncod, nsod*nl_d), INTENT(IN) :: sphi_d
569 : REAL(dp), TARGET :: work(ncoa*ncob*ncoc*ncod)
570 : REAL(dp) :: work2(ncoa, ncob, ncoc, ncod)
571 : REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2
572 : REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*&
573 : nl_c, nsod*nl_d) :: primitives_tmp
574 :
575 : INTEGER :: a_mysize(1), angular_momenta(4), i, &
576 : m_max, mysize, permutation(4)
577 : LOGICAL :: direct_contract
578 170446139 : REAL(dp), DIMENSION(:), POINTER :: eri_src, p_work
579 :
580 30110522 : IF (nproducts == 0) RETURN
581 :
582 170446139 : m_max = n_a + n_b + n_c + n_d
583 170446139 : mysize = ncoa*ncob*ncoc*ncod
584 340892278 : a_mysize = mysize
585 :
586 170446139 : IF (m_max /= 0) THEN
587 255258163 : IF (nproducts > 1) work(1:mysize) = 0.0_dp
588 134992516 : CALL get_eri_permutation(n_a, n_b, n_c, n_d, permutation)
589 674962580 : angular_momenta = [n_a, n_b, n_c, n_d]
590 1214932644 : angular_momenta = angular_momenta(permutation)
591 :
592 522254581 : DO i = 1, nproducts
593 387262065 : CALL cp_libint_set_params_eri(context%eri, pgf_product_list(i), permutation, m_max)
594 : CALL cp_libint_get_eris(angular_momenta(4), angular_momenta(3), angular_momenta(2), &
595 387262065 : angular_momenta(1), context%eri, p_work, a_mysize)
596 7821928971 : IF (nproducts > 1) work(1:mysize) = work(1:mysize) + p_work(1:mysize)
597 522254581 : neris = neris + mysize
598 : END DO
599 134992516 : IF (nproducts == 1) THEN
600 128103568 : eri_src => p_work
601 : ELSE
602 6888948 : eri_src => work
603 : END IF
604 3575007665 : DO i = 1, mysize
605 3575007665 : tmp_max = MAX(tmp_max, ABS(eri_src(i)))
606 : END DO
607 134992516 : tmp_max = tmp_max*max_contraction
608 134992516 : IF (tmp_max < eps_schwarz) RETURN
609 538257990 : CALL reorder_eri_block(eri_src, work2, [ncoa, ncob, ncoc, ncod], permutation, 1, 1)
610 : ELSE
611 35453623 : work(1) = 0.0_dp
612 158838209 : DO i = 1, nproducts
613 123384586 : work(1) = work(1) + pgf_product_list(i)%Fm(1)
614 158838209 : neris = neris + mysize
615 : END DO
616 35453623 : work2(1, 1, 1, 1) = work(1)
617 35453623 : tmp_max = max_contraction*ABS(work(1))
618 35453623 : IF (tmp_max < eps_schwarz) RETURN
619 : END IF
620 :
621 : direct_contract = (s_offset_a == 0 .AND. s_offset_b == 0 .AND. &
622 : s_offset_c == 0 .AND. s_offset_d == 0 .AND. &
623 : nsgfa == nsoa*nl_a .AND. nsgfb == nsob*nl_b .AND. &
624 140335617 : nsgfc == nsoc*nl_c .AND. nsgfd == nsod*nl_d)
625 : IF (direct_contract) THEN
626 : CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
627 : n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2, &
628 : sphi_a, &
629 : sphi_b, &
630 : sphi_c, &
631 : sphi_d, &
632 : primitives, &
633 30522529 : buffer1, buffer2)
634 : ELSE
635 8109957757 : primitives_tmp = 0.0_dp
636 :
637 : CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
638 : n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2, &
639 : sphi_a, &
640 : sphi_b, &
641 : sphi_c, &
642 : sphi_d, &
643 : primitives_tmp, &
644 109813088 : buffer1, buffer2)
645 :
646 : primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
647 : s_offset_b + 1:s_offset_b + nsob*nl_b, &
648 : s_offset_c + 1:s_offset_c + nsoc*nl_c, &
649 : s_offset_d + 1:s_offset_d + nsod*nl_d) = &
650 : primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
651 : s_offset_b + 1:s_offset_b + nsob*nl_b, &
652 : s_offset_c + 1:s_offset_c + nsoc*nl_c, &
653 8109957757 : s_offset_d + 1:s_offset_d + nsod*nl_d) + primitives_tmp(:, :, :, :)
654 : END IF
655 :
656 170446139 : END SUBROUTINE evaluate_eri
657 :
658 0 : END MODULE four_center_integral_interface
|
