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 2- and 3-center electron repulsion integral routines based on libint2
10 : !> Currently available operators: Coulomb, Truncated Coulomb, Short Range (erfc), Overlap
11 : !> \author A. Bussy (05.2019)
12 : ! **************************************************************************************************
13 :
14 : MODULE libint_2c_3c
15 : USE input_constants, ONLY: do_potential_coulomb,&
16 : do_potential_id,&
17 : do_potential_mix_cl_trunc,&
18 : do_potential_short,&
19 : do_potential_truncated
20 : USE integral_library_types, ONLY: compare_potential_types,&
21 : coulomb_operator_type,&
22 : libint_potential_type => coulomb_operator_type
23 : USE kinds, ONLY: dp
24 : USE libint_wrapper, ONLY: &
25 : cp_libint_get_2eri_derivs, cp_libint_get_2eris, cp_libint_get_3eri_derivs, &
26 : cp_libint_get_3eris, cp_libint_set_params_eri, cp_libint_set_params_eri_deriv, &
27 : cp_libint_t, eri_primitive_data_type, prim_data_f_size
28 : USE mathconstants, ONLY: pi
29 : USE orbital_pointers, ONLY: nco,&
30 : ncoset
31 : USE r12_fundamentals, ONLY: compute_r12_fundamentals
32 : USE t_c_g0, ONLY: get_lmax_init
33 : #include "./base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libint_2c_3c'
39 :
40 : PUBLIC :: eri_2center, eri_3center, cutoff_screen_factor, libint_potential_type, &
41 : eri_3center_derivs, eri_2center_derivs, compare_potential_types
42 :
43 : ! For screening of integrals with a truncated potential, it is important to use a slightly larger
44 : ! cutoff radius due to the discontinuity of the truncated Coulomb potential at the cutoff radius.
45 : REAL(KIND=dp), PARAMETER :: cutoff_screen_factor = 1.0001_dp
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian
51 : !> gaussian orbitals
52 : !> \param int_abc the integrals as array of cartesian orbitals (allocated before hand)
53 : !> \param la_min ...
54 : !> \param la_max ...
55 : !> \param npgfa ...
56 : !> \param zeta ...
57 : !> \param rpgfa ...
58 : !> \param ra ...
59 : !> \param lb_min ...
60 : !> \param lb_max ...
61 : !> \param npgfb ...
62 : !> \param zetb ...
63 : !> \param rpgfb ...
64 : !> \param rb ...
65 : !> \param lc_min ...
66 : !> \param lc_max ...
67 : !> \param npgfc ...
68 : !> \param zetc ...
69 : !> \param rpgfc ...
70 : !> \param rc ...
71 : !> \param dab ...
72 : !> \param dac ...
73 : !> \param dbc ...
74 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
75 : !> \param potential_parameter the info about the potential
76 : !> \param int_abc_ext the extremal value of int_abc, i.e., MAXVAL(ABS(int_abc))
77 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
78 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
79 : !> the latter must be initialized too
80 : ! **************************************************************************************************
81 9668738 : SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
82 9668738 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
83 9668738 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
84 : dab, dac, dbc, lib, potential_parameter, &
85 : int_abc_ext)
86 :
87 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: int_abc
88 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
89 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
90 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
91 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
92 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
93 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
94 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
95 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
96 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
97 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
98 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
99 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
100 : REAL(dp), INTENT(INOUT), OPTIONAL :: int_abc_ext
101 :
102 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, ipgf, j, &
103 : jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
104 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
105 9668738 : REAL(dp), DIMENSION(:), POINTER :: p_work
106 : TYPE(eri_primitive_data_type) :: params
107 :
108 9668738 : NULLIFY (p_work)
109 :
110 9668738 : dr_ab = 0.0_dp
111 9668738 : dr_bc = 0.0_dp
112 9668738 : dr_ac = 0.0_dp
113 :
114 9668738 : op = potential_parameter%potential_type
115 :
116 : IF (op == do_potential_truncated .OR. op == do_potential_short &
117 1269736 : .OR. op == do_potential_mix_cl_trunc) THEN
118 8690143 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
119 8690143 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
120 978595 : ELSE IF (op == do_potential_coulomb) THEN
121 158765 : dr_bc = 1000000.0_dp
122 158765 : dr_ac = 1000000.0_dp
123 : END IF
124 :
125 9668738 : IF (PRESENT(int_abc_ext)) THEN
126 9497196 : int_abc_ext = 0.0_dp
127 : END IF
128 :
129 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
130 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
131 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
132 :
133 : !Looping over the pgfs
134 30114868 : DO ipgf = 1, npgfa
135 20446130 : zeti = zeta(ipgf)
136 20446130 : a_start = (ipgf - 1)*ncoset(la_max)
137 :
138 86207078 : DO jpgf = 1, npgfb
139 :
140 : ! screening
141 56092210 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
142 :
143 26190984 : zetj = zetb(jpgf)
144 26190984 : b_start = (jpgf - 1)*ncoset(lb_max)
145 :
146 85844343 : DO kpgf = 1, npgfc
147 :
148 : ! screening
149 39207229 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
150 28673481 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
151 :
152 23890360 : zetk = zetc(kpgf)
153 23890360 : c_start = (kpgf - 1)*ncoset(lc_max)
154 :
155 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
156 : CALL prepare_3c_primitive(ra, rb, rc, zeti, zetj, zetk, &
157 23890360 : la_max + lb_max + lc_max, potential_parameter, params)
158 23890360 : CALL cp_libint_set_params_eri(lib, params, [1, 2, 3, 4], params%prepared_max_order)
159 :
160 56153381 : DO li = la_min, la_max
161 32263021 : a_offset = a_start + ncoset(li - 1)
162 32263021 : ncoa = nco(li)
163 90703755 : DO lj = MAX(li, lb_min), lb_max
164 34550374 : b_offset = b_start + ncoset(lj - 1)
165 34550374 : ncob = nco(lj)
166 110994170 : DO lk = lc_min, lc_max
167 44180775 : c_offset = c_start + ncoset(lk - 1)
168 44180775 : ncoc = nco(lk)
169 :
170 44180775 : a_mysize(1) = ncoa*ncob*ncoc
171 44180775 : CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
172 :
173 78731149 : IF (PRESENT(int_abc_ext)) THEN
174 159421579 : DO k = 1, ncoc
175 116465516 : p1 = (k - 1)*ncob
176 424762066 : DO j = 1, ncob
177 265340487 : p2 = (p1 + j - 1)*ncoa
178 835441538 : DO i = 1, ncoa
179 453635535 : p3 = p2 + i
180 453635535 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
181 718976022 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
182 : END DO
183 : END DO
184 : END DO
185 : ELSE
186 5588697 : DO k = 1, ncoc
187 4363985 : p1 = (k - 1)*ncob
188 15128846 : DO j = 1, ncob
189 9540149 : p2 = (p1 + j - 1)*ncoa
190 29175759 : DO i = 1, ncoa
191 15271625 : p3 = p2 + i
192 24811774 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
193 : END DO
194 : END DO
195 : END DO
196 : END IF
197 :
198 : END DO !lk
199 : END DO !lj
200 : END DO !li
201 :
202 : !swap centers 3 and 4 to compute (c|ab) with lb < la
203 23890360 : CALL cp_libint_set_params_eri(lib, params, [1, 2, 4, 3], params%prepared_max_order)
204 :
205 112256621 : DO lj = lb_min, lb_max
206 32274051 : b_offset = b_start + ncoset(lj - 1)
207 32274051 : ncob = nco(lj)
208 82851190 : DO li = MAX(lj + 1, la_min), la_max
209 11369910 : a_offset = a_start + ncoset(li - 1)
210 11369910 : ncoa = nco(li)
211 58324069 : DO lk = lc_min, lc_max
212 14680108 : c_offset = c_start + ncoset(lk - 1)
213 14680108 : ncoc = nco(lk)
214 :
215 14680108 : a_mysize(1) = ncoa*ncob*ncoc
216 14680108 : CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
217 :
218 26050018 : IF (PRESENT(int_abc_ext)) THEN
219 52991527 : DO k = 1, ncoc
220 38684442 : p1 = (k - 1)*ncoa
221 198382643 : DO i = 1, ncoa
222 145391116 : p2 = (p1 + i - 1)*ncob
223 375956234 : DO j = 1, ncob
224 191880676 : p3 = p2 + j
225 191880676 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
226 337271792 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
227 : END DO
228 : END DO
229 : END DO
230 : ELSE
231 1726525 : DO k = 1, ncoc
232 1353502 : p1 = (k - 1)*ncoa
233 6998668 : DO i = 1, ncoa
234 5272143 : p2 = (p1 + i - 1)*ncob
235 13612648 : DO j = 1, ncob
236 6987003 : p3 = p2 + j
237 12259146 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
238 : END DO
239 : END DO
240 : END DO
241 : END IF
242 :
243 : END DO !lk
244 : END DO !li
245 : END DO !lj
246 :
247 : END DO !kpgf
248 : END DO !jpgf
249 : END DO !ipgf
250 :
251 541449328 : END SUBROUTINE eri_3center
252 :
253 : ! **************************************************************************************************
254 : !> \brief Prepare geometry and radial seeds for 3-center integrals at the requested order.
255 : !> \param ri ...
256 : !> \param rj ...
257 : !> \param rk ...
258 : !> \param zeti ...
259 : !> \param zetj ...
260 : !> \param zetk ...
261 : !> \param m_max maximum radial order, including any derivative increment
262 : !> \param potential_parameter radial operator parameters
263 : !> \param params reusable primitive parameter cache
264 : ! **************************************************************************************************
265 25762487 : SUBROUTINE prepare_3c_primitive(ri, rj, rk, zeti, zetj, zetk, m_max, potential_parameter, params)
266 : REAL(dp), INTENT(IN) :: ri(3), rj(3), rk(3), zeti, zetj, zetk
267 : INTEGER, INTENT(IN) :: m_max
268 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
269 : TYPE(eri_primitive_data_type), INTENT(INOUT) :: params
270 :
271 : REAL(dp) :: fundamentals_work(prim_data_f_size), gammaq, normalization, prefac, S1234, T
272 :
273 25762487 : CPASSERT(m_max < prim_data_f_size)
274 :
275 : ! Libint layout: (k, dummy k | j, i).
276 103049948 : params%centers(:, 1) = rk
277 103049948 : params%centers(:, 2) = rk
278 103049948 : params%centers(:, 3) = rj
279 103049948 : params%centers(:, 4) = ri
280 128812435 : params%exponents = [zetk, 0.0_dp, zetj, zeti]
281 103049948 : params%P = rk
282 :
283 25762487 : params%prepared_max_order = m_max
284 25762487 : gammaq = zeti + zetj
285 25762487 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
286 25762487 : params%ZetapEtaInv = 1._dp/(zetk + gammaq)
287 :
288 103049948 : params%Q = (zeti*ri + zetj*rj)*params%EtaInv
289 103049948 : params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
290 25762487 : params%Rho = zetk*gammaq/(zetk + gammaq)
291 :
292 566774714 : params%Fm = 0.0_dp
293 49396307 : SELECT CASE (potential_parameter%potential_type)
294 : CASE (do_potential_coulomb, do_potential_truncated, do_potential_short, do_potential_mix_cl_trunc)
295 23633820 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
296 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
297 20590854 : CPASSERT(get_lmax_init() >= params%prepared_max_order)
298 : END IF
299 94535280 : T = params%Rho*SUM((params%Q - rk)**2)
300 94535280 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
301 : CALL compute_r12_fundamentals(params%Fm, normalization, params%prepared_max_order, potential_parameter%potential_type, &
302 : params%Rho, T, fundamentals_work, &
303 : omega=potential_parameter%omega, &
304 : cutoff_radius=potential_parameter%cutoff_radius, &
305 : scale_coulomb=potential_parameter%scale_coulomb, &
306 : scale_longrange=potential_parameter%scale_longrange, &
307 23633820 : scale_gaussian=0.0_dp)
308 23633820 : prefac = normalization*SQRT((pi*params%ZetapEtaInv)**3)*S1234
309 519944040 : params%Fm = prefac*params%Fm
310 : CASE (do_potential_id)
311 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
312 14900669 : - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
313 2128667 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
314 :
315 46830674 : params%Fm(:) = prefac
316 : CASE DEFAULT
317 25762487 : CPABORT("Requested operator NYI")
318 : END SELECT
319 25762487 : END SUBROUTINE prepare_3c_primitive
320 :
321 : ! **************************************************************************************************
322 : !> \brief Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given
323 : !> set of cartesian gaussian orbitals. Returns x,y,z derivatives for 1st and 2nd center
324 : !> \param der_abc_1 the derivatives for the 1st center (allocated before hand)
325 : !> \param der_abc_2 the derivatives for the 2nd center (allocated before hand)
326 : !> \param la_min ...
327 : !> \param la_max ...
328 : !> \param npgfa ...
329 : !> \param zeta ...
330 : !> \param rpgfa ...
331 : !> \param ra ...
332 : !> \param lb_min ...
333 : !> \param lb_max ...
334 : !> \param npgfb ...
335 : !> \param zetb ...
336 : !> \param rpgfb ...
337 : !> \param rb ...
338 : !> \param lc_min ...
339 : !> \param lc_max ...
340 : !> \param npgfc ...
341 : !> \param zetc ...
342 : !> \param rpgfc ...
343 : !> \param rc ...
344 : !> \param dab ...
345 : !> \param dac ...
346 : !> \param dbc ...
347 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
348 : !> \param potential_parameter the info about the potential
349 : !> \param der_abc_1_ext the extremal value of der_abc_1, i.e., MAXVAL(ABS(der_abc_1))
350 : !> \param der_abc_2_ext ...
351 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
352 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
353 : !> the latter must be initialized too. Note that the derivative wrt to the third center
354 : !> can be obtained via translational invariance
355 : ! **************************************************************************************************
356 345920 : SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
357 345920 : la_min, la_max, npgfa, zeta, rpgfa, ra, &
358 345920 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
359 345920 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
360 : dab, dac, dbc, lib, potential_parameter, &
361 : der_abc_1_ext, der_abc_2_ext)
362 :
363 : REAL(dp), DIMENSION(:, :, :, :), INTENT(INOUT) :: der_abc_1, der_abc_2
364 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
365 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
366 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
367 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
368 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
369 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
370 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
371 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
372 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
373 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
374 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
375 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
376 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: der_abc_1_ext, der_abc_2_ext
377 :
378 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, i_deriv, &
379 : ipgf, j, jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
380 : INTEGER, DIMENSION(3) :: permute_1, permute_2
381 : LOGICAL :: do_ext
382 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
383 : REAL(dp), DIMENSION(3) :: der_abc_1_ext_prv, der_abc_2_ext_prv
384 345920 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
385 : TYPE(eri_primitive_data_type) :: params
386 :
387 345920 : NULLIFY (p_deriv)
388 :
389 345920 : permute_1 = [4, 5, 6]
390 345920 : permute_2 = [7, 8, 9]
391 :
392 345920 : dr_ab = 0.0_dp
393 345920 : dr_bc = 0.0_dp
394 345920 : dr_ac = 0.0_dp
395 :
396 345920 : op = potential_parameter%potential_type
397 :
398 : IF (op == do_potential_truncated .OR. op == do_potential_short &
399 226113 : .OR. op == do_potential_mix_cl_trunc) THEN
400 119807 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
401 119807 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
402 226113 : ELSE IF (op == do_potential_coulomb) THEN
403 9988 : dr_bc = 1000000.0_dp
404 9988 : dr_ac = 1000000.0_dp
405 : END IF
406 :
407 345920 : do_ext = .FALSE.
408 345920 : IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
409 345920 : der_abc_1_ext_prv = 0.0_dp
410 345920 : der_abc_2_ext_prv = 0.0_dp
411 :
412 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
413 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
414 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
415 :
416 : !Looping over the pgfs
417 1189313 : DO ipgf = 1, npgfa
418 843393 : zeti = zeta(ipgf)
419 843393 : a_start = (ipgf - 1)*ncoset(la_max)
420 :
421 4761715 : DO jpgf = 1, npgfb
422 :
423 : ! screening
424 3572402 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
425 :
426 1134892 : zetj = zetb(jpgf)
427 1134892 : b_start = (jpgf - 1)*ncoset(lb_max)
428 :
429 8323219 : DO kpgf = 1, npgfc
430 :
431 : ! screening
432 6344934 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
433 2667789 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
434 :
435 1872127 : zetk = zetc(kpgf)
436 1872127 : c_start = (kpgf - 1)*ncoset(lc_max)
437 :
438 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
439 : CALL prepare_3c_primitive(ra, rb, rc, zeti, zetj, zetk, &
440 1872127 : la_max + lb_max + lc_max + 1, potential_parameter, params)
441 1872127 : CALL cp_libint_set_params_eri_deriv(lib, params, [1, 2, 3, 4], params%prepared_max_order)
442 :
443 4310369 : DO li = la_min, la_max
444 2438242 : a_offset = a_start + ncoset(li - 1)
445 2438242 : ncoa = nco(li)
446 6913203 : DO lj = MAX(li, lb_min), lb_max
447 2602834 : b_offset = b_start + ncoset(lj - 1)
448 2602834 : ncob = nco(lj)
449 8806530 : DO lk = lc_min, lc_max
450 3765454 : c_offset = c_start + ncoset(lk - 1)
451 3765454 : ncoc = nco(lk)
452 :
453 3765454 : a_mysize(1) = ncoa*ncob*ncoc
454 :
455 3765454 : CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
456 :
457 3765454 : IF (do_ext) THEN
458 15061816 : DO i_deriv = 1, 3
459 42689449 : DO k = 1, ncoc
460 27627633 : p1 = (k - 1)*ncob
461 90942462 : DO j = 1, ncob
462 52018467 : p2 = (p1 + j - 1)*ncoa
463 158128833 : DO i = 1, ncoa
464 78482733 : p3 = p2 + i
465 :
466 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
467 78482733 : p_deriv(p3, permute_2(i_deriv))
468 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
469 78482733 : ABS(p_deriv(p3, permute_2(i_deriv))))
470 :
471 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
472 78482733 : p_deriv(p3, permute_1(i_deriv))
473 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
474 130501200 : ABS(p_deriv(p3, permute_1(i_deriv))))
475 :
476 : END DO
477 : END DO
478 : END DO
479 : END DO
480 : ELSE
481 0 : DO i_deriv = 1, 3
482 0 : DO k = 1, ncoc
483 0 : p1 = (k - 1)*ncob
484 0 : DO j = 1, ncob
485 0 : p2 = (p1 + j - 1)*ncoa
486 0 : DO i = 1, ncoa
487 0 : p3 = p2 + i
488 :
489 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
490 0 : p_deriv(p3, permute_2(i_deriv))
491 :
492 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
493 0 : p_deriv(p3, permute_1(i_deriv))
494 : END DO
495 : END DO
496 : END DO
497 : END DO
498 : END IF
499 :
500 6368288 : DEALLOCATE (p_deriv)
501 : END DO !lk
502 : END DO !lj
503 : END DO !li
504 :
505 : !swap centers 3 and 4 to compute (c|ab) with lb < la
506 1872127 : CALL cp_libint_set_params_eri_deriv(lib, params, [1, 2, 4, 3], params%prepared_max_order)
507 :
508 7891974 : DO lj = lb_min, lb_max
509 2447445 : b_offset = b_start + ncoset(lj - 1)
510 2447445 : ncob = nco(lj)
511 9456389 : DO li = MAX(lj + 1, la_min), la_max
512 664010 : a_offset = a_start + ncoset(li - 1)
513 664010 : ncoa = nco(li)
514 4084604 : DO lk = lc_min, lc_max
515 973149 : c_offset = c_start + ncoset(lk - 1)
516 973149 : ncoc = nco(lk)
517 :
518 973149 : a_mysize(1) = ncoa*ncob*ncoc
519 973149 : CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
520 :
521 973149 : IF (do_ext) THEN
522 3892596 : DO i_deriv = 1, 3
523 11259984 : DO k = 1, ncoc
524 7367388 : p1 = (k - 1)*ncoa
525 34471833 : DO i = 1, ncoa
526 24184998 : p2 = (p1 + i - 1)*ncob
527 59225208 : DO j = 1, ncob
528 27672822 : p3 = p2 + j
529 :
530 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
531 27672822 : p_deriv(p3, permute_1(i_deriv))
532 :
533 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
534 27672822 : ABS(p_deriv(p3, permute_1(i_deriv))))
535 :
536 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
537 27672822 : p_deriv(p3, permute_2(i_deriv))
538 :
539 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
540 51857820 : ABS(p_deriv(p3, permute_2(i_deriv))))
541 : END DO
542 : END DO
543 : END DO
544 : END DO
545 : ELSE
546 0 : DO i_deriv = 1, 3
547 0 : DO k = 1, ncoc
548 0 : p1 = (k - 1)*ncoa
549 0 : DO i = 1, ncoa
550 0 : p2 = (p1 + i - 1)*ncob
551 0 : DO j = 1, ncob
552 0 : p3 = p2 + j
553 :
554 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
555 0 : p_deriv(p3, permute_1(i_deriv))
556 :
557 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
558 0 : p_deriv(p3, permute_2(i_deriv))
559 : END DO
560 : END DO
561 : END DO
562 : END DO
563 : END IF
564 :
565 1637159 : DEALLOCATE (p_deriv)
566 : END DO !lk
567 : END DO !li
568 : END DO !lj
569 :
570 : END DO !kpgf
571 : END DO !jpgf
572 : END DO !ipgf
573 :
574 345920 : IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
575 345920 : IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
576 :
577 19371520 : END SUBROUTINE eri_3center_derivs
578 :
579 : ! **************************************************************************************************
580 : !> \brief Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian
581 : !> gaussian orbitals
582 : !> \param int_ab the integrals as array of cartesian orbitals (allocated before hand)
583 : !> \param la_min ...
584 : !> \param la_max ...
585 : !> \param npgfa ...
586 : !> \param zeta ...
587 : !> \param rpgfa ...
588 : !> \param ra ...
589 : !> \param lb_min ...
590 : !> \param lb_max ...
591 : !> \param npgfb ...
592 : !> \param zetb ...
593 : !> \param rpgfb ...
594 : !> \param rb ...
595 : !> \param dab ...
596 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
597 : !> \param potential_parameter the info about the potential
598 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
599 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
600 : !> the latter must be initialized too
601 : ! **************************************************************************************************
602 596701 : SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
603 596701 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
604 : dab, lib, potential_parameter)
605 :
606 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int_ab
607 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
608 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
609 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
610 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
611 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
612 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
613 : REAL(dp), INTENT(IN) :: dab
614 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
615 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
616 :
617 : INTEGER :: a_mysize(1), a_offset, a_start, &
618 : b_offset, b_start, i, ipgf, j, jpgf, &
619 : li, lj, ncoa, ncob, p1, p2
620 : REAL(dp) :: dr_ab, zeti, zetj
621 596701 : REAL(dp), DIMENSION(:), POINTER :: p_work
622 : TYPE(eri_primitive_data_type) :: params
623 :
624 596701 : NULLIFY (p_work)
625 :
626 596701 : dr_ab = 0.0_dp
627 :
628 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
629 313423 : potential_parameter%potential_type == do_potential_short .OR. &
630 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
631 295914 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
632 300787 : ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
633 124370 : dr_ab = 1000000.0_dp
634 : END IF
635 :
636 : !Looping over the pgfs
637 1695576 : DO ipgf = 1, npgfa
638 1098875 : zeti = zeta(ipgf)
639 1098875 : a_start = (ipgf - 1)*ncoset(la_max)
640 :
641 6361038 : DO jpgf = 1, npgfb
642 4665462 : zetj = zetb(jpgf)
643 4665462 : b_start = (jpgf - 1)*ncoset(lb_max)
644 :
645 : !screening
646 4665462 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
647 :
648 : CALL prepare_2c_primitive(ra, rb, zeti, zetj, la_max + lb_max, &
649 3049683 : potential_parameter, params)
650 3049683 : CALL cp_libint_set_params_eri(lib, params, [1, 2, 3, 4], params%prepared_max_order)
651 :
652 9647647 : DO li = la_min, la_max
653 5499089 : a_offset = a_start + ncoset(li - 1)
654 5499089 : ncoa = nco(li)
655 20869017 : DO lj = lb_min, lb_max
656 10704466 : b_offset = b_start + ncoset(lj - 1)
657 10704466 : ncob = nco(lj)
658 :
659 10704466 : a_mysize(1) = ncoa*ncob
660 10704466 : CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
661 :
662 48899708 : DO j = 1, ncob
663 32696153 : p1 = (j - 1)*ncoa
664 145169550 : DO i = 1, ncoa
665 101768931 : p2 = p1 + i
666 134465084 : int_ab(a_offset + i, b_offset + j) = p_work(p2)
667 : END DO
668 : END DO
669 :
670 : END DO
671 : END DO
672 :
673 : END DO
674 : END DO
675 :
676 33415256 : END SUBROUTINE eri_2center
677 :
678 : ! **************************************************************************************************
679 : !> \brief Prepare geometry and radial seeds for 2-center integrals at the requested order.
680 : !> \param rj ...
681 : !> \param rk ...
682 : !> \param zetj ...
683 : !> \param zetk ...
684 : !> \param m_max maximum radial order, including any derivative increment
685 : !> \param potential_parameter radial operator parameters
686 : !> \param params reusable primitive parameter cache
687 : ! **************************************************************************************************
688 5225304 : SUBROUTINE prepare_2c_primitive(rj, rk, zetj, zetk, m_max, potential_parameter, params)
689 : REAL(dp), INTENT(IN) :: rj(3), rk(3), zetj, zetk
690 : INTEGER, INTENT(IN) :: m_max
691 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
692 : TYPE(eri_primitive_data_type), INTENT(INOUT) :: params
693 :
694 : REAL(dp) :: fundamentals_work(prim_data_f_size), normalization, prefac, T
695 :
696 5225304 : CPASSERT(m_max < prim_data_f_size)
697 :
698 : ! Libint layout: (k, dummy k | j, dummy j).
699 20901216 : params%centers(:, 1) = rk
700 20901216 : params%centers(:, 2) = rk
701 20901216 : params%centers(:, 3) = rj
702 20901216 : params%centers(:, 4) = rj
703 26126520 : params%exponents = [zetk, 0.0_dp, zetj, 0.0_dp]
704 20901216 : params%P = rk
705 20901216 : params%Q = rj
706 :
707 5225304 : params%prepared_max_order = m_max
708 5225304 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
709 5225304 : params%ZetapEtaInv = 1._dp/(zetk + zetj)
710 :
711 20901216 : params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
712 5225304 : params%Rho = zetk*zetj/(zetk + zetj)
713 :
714 114956688 : params%Fm = 0.0_dp
715 10282162 : SELECT CASE (potential_parameter%potential_type)
716 : CASE (do_potential_coulomb, do_potential_truncated, do_potential_short, do_potential_mix_cl_trunc)
717 5056858 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
718 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
719 351319 : CPASSERT(get_lmax_init() >= params%prepared_max_order)
720 : END IF
721 20227432 : T = params%Rho*SUM((rj - rk)**2)
722 : CALL compute_r12_fundamentals(params%Fm, normalization, params%prepared_max_order, potential_parameter%potential_type, &
723 : params%Rho, T, fundamentals_work, &
724 : omega=potential_parameter%omega, &
725 : cutoff_radius=potential_parameter%cutoff_radius, &
726 : scale_coulomb=potential_parameter%scale_coulomb, &
727 : scale_longrange=potential_parameter%scale_longrange, &
728 5056858 : scale_gaussian=0.0_dp)
729 5056858 : prefac = normalization*SQRT((pi*params%ZetapEtaInv)**3)
730 111250876 : params%Fm = prefac*params%Fm
731 : CASE (do_potential_id)
732 :
733 673784 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
734 3705812 : params%Fm(:) = prefac
735 : CASE DEFAULT
736 5225304 : CPABORT("Requested operator NYI")
737 : END SELECT
738 5225304 : END SUBROUTINE prepare_2c_primitive
739 :
740 : !> \brief Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given
741 : !> set of cartesian gaussian orbitals. Returns the derivatives wrt to the first center
742 : !> \param der_ab the derivatives as array of cartesian orbitals (allocated before hand)
743 : !> \param la_min ...
744 : !> \param la_max ...
745 : !> \param npgfa ...
746 : !> \param zeta ...
747 : !> \param rpgfa ...
748 : !> \param ra ...
749 : !> \param lb_min ...
750 : !> \param lb_max ...
751 : !> \param npgfb ...
752 : !> \param zetb ...
753 : !> \param rpgfb ...
754 : !> \param rb ...
755 : !> \param dab ...
756 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
757 : !> \param potential_parameter the info about the potential
758 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
759 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
760 : !> the latter must be initialized too
761 : ! **************************************************************************************************
762 144239 : SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
763 144239 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
764 : dab, lib, potential_parameter)
765 :
766 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: der_ab
767 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
768 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
769 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
770 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
771 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
772 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
773 : REAL(dp), INTENT(IN) :: dab
774 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
775 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
776 :
777 : INTEGER :: a_mysize(1), a_offset, a_start, &
778 : b_offset, b_start, i, i_deriv, ipgf, &
779 : j, jpgf, li, lj, ncoa, ncob, p1, p2
780 : INTEGER, DIMENSION(3) :: permute
781 : REAL(dp) :: dr_ab, zeti, zetj
782 144239 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
783 : TYPE(eri_primitive_data_type) :: params
784 :
785 144239 : NULLIFY (p_deriv)
786 :
787 144239 : permute = [4, 5, 6]
788 :
789 144239 : dr_ab = 0.0_dp
790 :
791 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
792 80471 : potential_parameter%potential_type == do_potential_short .OR. &
793 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
794 64460 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
795 79779 : ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
796 9831 : dr_ab = 1000000.0_dp
797 : END IF
798 :
799 : !Looping over the pgfs
800 641174 : DO ipgf = 1, npgfa
801 496935 : zeti = zeta(ipgf)
802 496935 : a_start = (ipgf - 1)*ncoset(la_max)
803 :
804 3875095 : DO jpgf = 1, npgfb
805 3233921 : zetj = zetb(jpgf)
806 3233921 : b_start = (jpgf - 1)*ncoset(lb_max)
807 :
808 : !screening
809 3233921 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
810 :
811 : CALL prepare_2c_primitive(ra, rb, zeti, zetj, la_max + lb_max + 1, &
812 2175621 : potential_parameter, params)
813 2175621 : CALL cp_libint_set_params_eri_deriv(lib, params, [1, 2, 3, 4], params%prepared_max_order)
814 :
815 6321018 : DO li = la_min, la_max
816 3648462 : a_offset = a_start + ncoset(li - 1)
817 3648462 : ncoa = nco(li)
818 13051067 : DO lj = lb_min, lb_max
819 6168684 : b_offset = b_start + ncoset(lj - 1)
820 6168684 : ncob = nco(lj)
821 :
822 6168684 : a_mysize(1) = ncoa*ncob
823 6168684 : CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
824 :
825 24674736 : DO i_deriv = 1, 3
826 75355947 : DO j = 1, ncob
827 50681211 : p1 = (j - 1)*ncoa
828 208451331 : DO i = 1, ncoa
829 139264068 : p2 = p1 + i
830 189945279 : der_ab(a_offset + i, b_offset + j, i_deriv) = p_deriv(p2, permute(i_deriv))
831 : END DO
832 : END DO
833 : END DO
834 :
835 9817146 : DEALLOCATE (p_deriv)
836 : END DO
837 : END DO
838 :
839 : END DO
840 : END DO
841 :
842 8077384 : END SUBROUTINE eri_2center_derivs
843 :
844 : END MODULE libint_2c_3c
|
