(git:f7b9015)
Loading...
Searching...
No Matches
four_center_integral_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
11MODULE 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
44CONTAINS
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 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 DIMENSION(:) :: tmp_product_list
87
88 nimages1 = list1%nimages
89 nimages2 = list2%nimages
90 nproducts = 0
91 IF (nimages1 == 0 .OR. nimages2 == 0) RETURN
92 CALL prepare_eri4_product(list1, list2, m_max, product)
93 etainv = list2%ZetaInv
94 zeta_c = list2%zeta
95 zeta_d = list2%zetb
96 temp_cc = 0.0_dp
97 temp_dd = 0.0_dp
98 r_cutoff_sq = 0.0_dp
99 r_cutoff_sq_guard = 0.0_dp
100 SELECT CASE (potential_parameter%potential_type)
101 CASE (do_potential_truncated, do_potential_short, do_potential_mix_cl_trunc)
102 has_r_cutoff = .true.
103 rcut_eff = potential_parameter%cutoff_radius
104 CASE (do_potential_tshpsc)
105 has_r_cutoff = .true.
106 rcut_eff = potential_parameter%cutoff_radius*2.0_dp
107 CASE DEFAULT
108 has_r_cutoff = .false.
109 rcut_eff = 0.0_dp
110 END SELECT
111 DO i = 1, nimages1
112 p = list1%image_list(i)%P
113 r1 = list1%image_list(i)%R
114 s1234a = list1%image_list(i)%S1234
115 pgf_max_1 = list1%image_list(i)%pgf_max
116 DO j = 1, nimages2
117 pgf_max_2 = list2%image_list(j)%pgf_max
118 IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) cycle
119 q = list2%image_list(j)%P
120 r2 = list2%image_list(j)%R
121 s1234b = list2%image_list(j)%S1234
122 rc = list2%image_list(j)%ra
123 rd = list2%image_list(j)%rb
124
125 IF (has_r_cutoff) THEN
126 r_cutoff_sq = (r1 + r2 + rcut_eff)**2
127 r_cutoff_sq_guard = r_cutoff_sq*(1.0_dp + r_cutoff_eps)
128 END IF
129 k_max = SIZE(neighbor_cells)
130 IF (do_periodic) THEN
131 temp = p - q
132 pq = pbc(temp, cell)
133 IF (has_r_cutoff .AND. cell%orthorhombic .AND. all(cell%perd == 1) .AND. k_max >= 1) THEN
134 IF (all(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
135 pq2 = pq(1)**2 + pq(2)**2 + pq(3)**2
136 IF (pq2 > r_cutoff_sq_guard) cycle
137 min_img_delta = min_image_gap_sq(pq, cell%hmat)
138 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 shift = -pq + temp
142 temp_cc = rc + shift
143 temp_dd = rd + shift
144 END IF
145
146 s1234 = exp(s1234a + s1234b)
147
148 DO k = 1, k_max
149 IF (do_periodic) THEN
150 IF (has_r_cutoff .AND. k_max > 1) THEN
151 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 c11 = temp_cc + neighbor_cells(k)%cell_r(:)
156 tmp_d = temp_dd + neighbor_cells(k)%cell_r(:)
157 ELSE
158 c11 = rc
159 tmp_d = rd
160 END IF
161 q = (zeta_c*c11 + zeta_d*tmp_d)*etainv
162 rpq2 = (p(1) - q(1))**2 + (p(2) - q(2))**2 + (p(3) - q(3))**2
163 IF (has_r_cutoff) THEN
164 IF (rpq2 > r_cutoff_sq) cycle
165 END IF
166 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 IF (nproducts > SIZE(product_list)) THEN
171!$OMP ATOMIC READ
172 tmp_i4 = pgf_product_list_size
173 tmp_i4 = max(pgf_product_list_size, (3*nproducts + 1)/2)
174!$OMP ATOMIC WRITE
175 pgf_product_list_size = tmp_i4
176 ALLOCATE (tmp_product_list(SIZE(product_list)))
177 tmp_product_list(:) = product_list
178 DEALLOCATE (product_list)
179 ALLOCATE (product_list(tmp_i4))
180 product_list(1:SIZE(tmp_product_list)) = tmp_product_list
181 DEALLOCATE (tmp_product_list)
182 END IF
183
184 shifted_image%ra = c11
185 shifted_image%rb = tmp_d
186 shifted_image%P = q
187 CALL prepare_eri4_primitive(product, list1%image_list(i), shifted_image, s1234, rpq2, &
188 potential_parameter, product_list(nproducts))
189 END DO
190 END DO
191 END DO
192
193 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 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 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 i_max = SIZE(neighbor_cells)
241 IF (do_periodic) THEN
242 IF (cell%orthorhombic .AND. screen1(1) <= 0.0_dp .AND. i_max > 1) THEN
243 IF (all(cell%perd == 1) .AND. all(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
244 ab = ra - rb
245 rab2 = ab(1)**2 + ab(2)**2 + ab(3)**2
246 min_img_delta = min_image_gap_sq(ab, cell%hmat)
247 IF (min_img_delta >= 0.0_dp) THEN
248 IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
249 nelements = 0
250 RETURN
251 END IF
252 IF (screen1(1)*(rab2 + min_img_delta)*(1.0_dp - min_img_eps) + &
253 screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) i_max = 1
254 END IF
255 END IF
256 END IF
257 ELSE
258 i_max = min(i_max, 1)
259 END IF
260 nimages = 0
261 ! ** inner loop may never be reached
262 nelements = npgfa*npgfb
263 DO i = 1, i_max
264 IF (do_periodic) THEN
265 im_b = rb + neighbor_cells(i)%cell_r(:)
266 ELSE
267 im_b = rb
268 END IF
269 ab = ra - im_b
270 rab2 = ab(1)**2 + ab(2)**2 + ab(3)**2
271 IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) cycle
272 element_counter = 0
273 DO ipgf = 1, npgfa
274 DO jpgf = 1, npgfb
275 element_counter = element_counter + 1
276 pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
277 IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
278 cycle
279 END IF
280 nimages(element_counter) = nimages(element_counter) + 1
281 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 list(element_counter), list(element_counter)%image_list(nimages(element_counter)))
284 IF (nimages(element_counter) == 1) THEN
285 list(element_counter)%ipgf = ipgf
286 list(element_counter)%jpgf = jpgf
287 END IF
288
289 list(element_counter)%image_list(nimages(element_counter))%R = &
290 max(0.0_dp, r_pgf(jpgf, ipgf)%x(1)*rab2 + r_pgf(jpgf, ipgf)%x(2))
291 list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
292 END DO
293 END DO
294 nelements = max(nelements, element_counter)
295 END DO
296 DO j = 1, nelements
297 list(j)%nimages = nimages(j)
298 END DO
299 ! ** Remove unused elements
300
301 element_counter = 0
302 DO j = 1, nelements
303 IF (list(j)%nimages == 0) cycle
304 element_counter = element_counter + 1
305 list(element_counter)%nimages = list(j)%nimages
306 list(element_counter)%zetapzetb = list(j)%zetapzetb
307 list(element_counter)%ZetaInv = list(j)%ZetaInv
308 list(element_counter)%zeta = list(j)%zeta
309 list(element_counter)%zetb = list(j)%zetb
310 list(element_counter)%ipgf = list(j)%ipgf
311 list(element_counter)%jpgf = list(j)%jpgf
312 DO i = 1, list(j)%nimages
313 list(element_counter)%image_list(i) = list(j)%image_list(i)
314 END DO
315 END DO
316
317 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 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 cell_len = [abs(hmat(1, 1)), abs(hmat(2, 2)), abs(hmat(3, 3))]
337 gap = minval(cell_len*(cell_len - 2.0_dp*abs(r)))
338
339 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 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 zeta, zetb, zetc, zetd, &
428 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 nsgfl_a, nsgfl_b, nsgfl_c, &
436 nsgfl_d, &
437 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 cart_estimate = 0.0_dp
480 neris_tmp = 0
481 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 neighbor_cells, nimages, do_periodic, cell)
487 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 neighbor_cells, nimages, do_periodic, cell)
492 IF (nelements_kl == 0) RETURN
493
494 max_l = la_max + lb_max + lc_max + ld_max
495
496 DO list_ij = 1, nelements_ij
497 ipgf = pgf_list_ij(list_ij)%ipgf
498 jpgf = pgf_list_ij(list_ij)%jpgf
499
500 DO list_kl = 1, nelements_kl
501 kpgf = pgf_list_kl(list_kl)%ipgf
502 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 potential_parameter, max_l, do_periodic)
507 IF (nproducts == 0) cycle
508
509 s_offset_a = 0
510 DO la = la_min, la_max
511 s_offset_b = 0
512 ncoa = nco(la)
513 nsgfla = nsgfl_a(la)
514 nsoa = nso(la)
515
516 DO lb = lb_min, lb_max
517 s_offset_c = 0
518 ncob = nco(lb)
519 nsgflb = nsgfl_b(lb)
520 nsob = nso(lb)
521
522 DO lc = lc_min, lc_max
523 s_offset_d = 0
524 ncoc = nco(lc)
525 nsgflc = nsgfl_c(lc)
526 nsoc = nso(lc)
527
528 DO ld = ld_min, ld_max
529 ncod = nco(ld)
530 nsgfld = nsgfl_d(ld)
531 nsod = nso(ld)
532
533 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 ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp)
548 cart_estimate = max(tmp_max, cart_estimate)
549 s_offset_d = s_offset_d + nsod*nsgfld
550 END DO !ld
551 s_offset_c = s_offset_c + nsoc*nsgflc
552 END DO !lc
553 s_offset_b = s_offset_b + nsob*nsgflb
554 END DO !lb
555 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 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 zeta, zetb, zetc, zetd, &
651 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 nsgfl_a, nsgfl_b, nsgfl_c, &
659 nsgfl_d, &
660 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 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 cart_estimate = 0.0_dp
713 primitive_integrals = 0.0_dp
714 IF (use_virial) THEN
715 primitive_integrals_virial = 0.0_dp
716 cart_estimate_virial = 0.0_dp
717 END IF
718 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 neighbor_cells, nimages, do_periodic, cell)
724 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 neighbor_cells, nimages, do_periodic, cell)
729 IF (nelements_kl == 0) RETURN
730
731 max_l = la_max + lb_max + lc_max + ld_max + 1
732 DO list_ij = 1, nelements_ij
733 ipgf = pgf_list_ij(list_ij)%ipgf
734 jpgf = pgf_list_ij(list_ij)%jpgf
735
736 DO list_kl = 1, nelements_kl
737 kpgf = pgf_list_kl(list_kl)%ipgf
738 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 potential_parameter, max_l, do_periodic)
743 IF (nproducts == 0) cycle
744 s_offset_a = 0
745 DO la = la_min, la_max
746 s_offset_b = 0
747 ncoa = nco(la)
748 nsgfla = nsgfl_a(la)
749 nsoa = nso(la)
750
751 DO lb = lb_min, lb_max
752 s_offset_c = 0
753 ncob = nco(lb)
754 nsgflb = nsgfl_b(lb)
755 nsob = nso(lb)
756
757 DO lc = lc_min, lc_max
758 s_offset_d = 0
759 ncoc = nco(lc)
760 nsgflc = nsgfl_c(lc)
761 nsoc = nso(lc)
762
763 DO ld = ld_min, ld_max
764 ncod = nco(ld)
765 nsgfld = nsgfl_d(ld)
766 nsod = nso(ld)
767 tmp_max = 0.0_dp
768 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 primitive_integrals_virial, cell, tmp_max_virial)
785 cart_estimate = max(tmp_max, cart_estimate)
786 IF (use_virial) cart_estimate_virial = max(tmp_max_virial, cart_estimate_virial)
787 s_offset_d = s_offset_d + nsod*nsgfld
788 END DO !ld
789 s_offset_c = s_offset_c + nsoc*nsgflc
790 END DO !lc
791 s_offset_b = s_offset_b + nsob*nsgflb
792 END DO !lb
793 s_offset_a = s_offset_a + nsoa*nsgfla
794 END DO !la
795 END DO
796 END DO
797
798 END SUBROUTINE forces4
799
800END MODULE four_center_integral_methods
cell_types::pbc
Definition cell_types.F:105
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
four_center_integral_interface
Four-center integral evaluation with an encapsulated Libint backend.
Definition four_center_integral_interface.F:15
four_center_integral_interface::evaluate_deriv_eri
subroutine, public evaluate_deriv_eri(context, nproducts, pgf_product_list, n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd, primitives, max_contraction, tmp_max_all, eps_schwarz, neris, s_offset_a, s_offset_b, s_offset_c, s_offset_d, nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, sphi_a, sphi_b, sphi_c, sphi_d, work, work2, work_forces, buffer1, buffer2, primitives_tmp, use_virial, work_virial, work2_virial, primitives_tmp_virial, primitives_virial, cell, tmp_max_all_virial)
Evaluates derivatives of electron repulsion integrals for a primitive quartet.
Definition four_center_integral_interface.F:221
four_center_integral_interface::evaluate_eri
subroutine, public evaluate_eri(context, nproducts, pgf_product_list, n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd, primitives, max_contraction, tmp_max, eps_schwarz, neris, s_offset_a, s_offset_b, s_offset_c, s_offset_d, nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, primitives_tmp)
Evaluate electron repulsion integrals for a primitive quartet.
Definition four_center_integral_interface.F:553
four_center_integral_methods
Prepare and evaluate screened four-center Gaussian integrals and derivatives.
Definition four_center_integral_methods.F:11
four_center_integral_methods::forces4
subroutine, public forces4(context, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, la_min, la_max, lb_min, lb_max, lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, zeta, zetb, zetc, zetd, primitive_integrals, potential_parameter, neighbor_cells, screen1, screen2, eps_schwarz, max_contraction_val, cart_estimate, cell, neris_tmp, log10_pmax, log10_eps_schwarz, r1_pgf, r2_pgf, pgf1, pgf2, pgf_list_ij, pgf_list_kl, pgf_product_list, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d, sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext, ede_work, ede_work2, ede_work_forces, ede_buffer1, ede_buffer2, ede_primitives_tmp, nimages, do_periodic, use_virial, ede_work_virial, ede_work2_virial, ede_primitives_tmp_virial, primitive_integrals_virial, cart_estimate_virial)
...
Definition four_center_integral_methods.F:666
four_center_integral_methods::min_image_gap_sq
pure real(dp) function, public min_image_gap_sq(r, hmat)
Smallest increase of |r|**2 when r is moved to any other periodic image of an orthorhombic,...
Definition four_center_integral_methods.F:330
four_center_integral_methods::coulomb4
subroutine, public coulomb4(context, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, la_min, la_max, lb_min, lb_max, lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, zeta, zetb, zetc, zetd, primitive_integrals, potential_parameter, neighbor_cells, screen1, screen2, eps_schwarz, max_contraction_val, cart_estimate, cell, neris_tmp, log10_pmax, log10_eps_schwarz, r1_pgf, r2_pgf, pgf1, pgf2, pgf_list_ij, pgf_list_kl, pgf_product_list, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d, sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext, ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, nimages, do_periodic)
calculates two-electron integrals of a quartet/shell using the library lib_int in the periodic case
Definition four_center_integral_methods.F:440
four_center_integral_methods::pgf_product_list_size
integer, save, public pgf_product_list_size
Definition four_center_integral_methods.F:42
four_center_integral_preparation
Gaussian-product geometry and radial seeds shared by four-center evaluation and screening.
Definition four_center_integral_preparation.F:11
four_center_integral_preparation::prepare_eri4_primitive
subroutine, public prepare_eri4_primitive(product, image_a, image_b, attenuation, pq2, potential_parameter, params)
Complete a prepared quartet for one selected image, including normalized radial seeds.
Definition four_center_integral_preparation.F:89
four_center_integral_preparation::prepare_eri4_pair
subroutine, public prepare_eri4_pair(ra, rb, zeta, zetb, rab2, pair, image)
Prepare one Gaussian pair, retaining the caller's image topology and screening metadata.
Definition four_center_integral_preparation.F:40
four_center_integral_preparation::prepare_eri4_product
subroutine, public prepare_eri4_product(pair_a, pair_b, m_max, product)
Prepare the image-independent quartet data once per combination of Gaussian pairs.
Definition four_center_integral_preparation.F:63
four_center_integral_types
Primitive Gaussian pairs and products for four-center integrals. The fixed radial workspace capacity ...
Definition four_center_integral_types.F:12
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:856
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:856
input_constants::do_potential_mix_cl_trunc
integer, parameter, public do_potential_mix_cl_trunc
Definition input_constants.F:856
input_constants::do_potential_tshpsc
integer, parameter, public do_potential_tshpsc
Definition input_constants.F:856
integral_library_types
Library choices for electronic integral APIs.
Definition integral_library_types.F:11
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
four_center_integral_interface::four_center_integral_context_type
Context owning the energy and first-derivative backend workspaces.
Definition four_center_integral_interface.F:48
four_center_integral_types::eri4_cell_type
Definition four_center_integral_types.F:24
four_center_integral_types::eri4_pgf_image_type
Definition four_center_integral_types.F:35
four_center_integral_types::eri4_pgf_pair_type
Definition four_center_integral_types.F:46
four_center_integral_types::eri4_screen_coeff_type
Definition four_center_integral_types.F:30
integral_library_types::coulomb_operator_type
Definition integral_library_types.F:39