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 gamma, ONLY: fgamma => fgamma_0
16 : USE input_constants, ONLY: do_potential_coulomb,&
17 : do_potential_id,&
18 : do_potential_long,&
19 : do_potential_mix_cl_trunc,&
20 : do_potential_short,&
21 : do_potential_truncated
22 : USE integral_library_types, ONLY: coulomb_operator_type,&
23 : libint_potential_type => coulomb_operator_type
24 : USE kinds, ONLY: dp
25 : USE libint_wrapper, ONLY: cp_libint_get_2eri_derivs,&
26 : cp_libint_get_2eris,&
27 : cp_libint_get_3eri_derivs,&
28 : cp_libint_get_3eris,&
29 : cp_libint_set_params_eri,&
30 : cp_libint_set_params_eri_deriv,&
31 : cp_libint_t,&
32 : prim_data_f_size
33 : USE mathconstants, ONLY: pi
34 : USE orbital_pointers, ONLY: nco,&
35 : ncoset
36 : USE t_c_g0, ONLY: get_lmax_init,&
37 : t_c_g0_n
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libint_2c_3c'
44 :
45 : PUBLIC :: eri_2center, eri_3center, cutoff_screen_factor, libint_potential_type, &
46 : eri_3center_derivs, eri_2center_derivs, compare_potential_types
47 :
48 : ! For screening of integrals with a truncated potential, it is important to use a slightly larger
49 : ! cutoff radius due to the discontinuity of the truncated Coulomb potential at the cutoff radius.
50 : REAL(KIND=dp), PARAMETER :: cutoff_screen_factor = 1.0001_dp
51 :
52 : TYPE :: params_2c
53 : INTEGER :: m_max = 0
54 : REAL(dp) :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
55 : REAL(dp), DIMENSION(3) :: W = 0.0_dp
56 : REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
57 : END TYPE params_2c
58 :
59 : TYPE :: params_3c
60 : INTEGER :: m_max = 0
61 : REAL(dp) :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
62 : REAL(dp), DIMENSION(3) :: Q = 0.0_dp, W = 0.0_dp
63 : REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
64 : END TYPE params_3c
65 :
66 : ! Compatibility alias retained for existing callers. New interfaces should use the
67 : ! backend-neutral name from integral_library_types.
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian
73 : !> gaussian orbitals
74 : !> \param int_abc the integrals as array of cartesian orbitals (allocated before hand)
75 : !> \param la_min ...
76 : !> \param la_max ...
77 : !> \param npgfa ...
78 : !> \param zeta ...
79 : !> \param rpgfa ...
80 : !> \param ra ...
81 : !> \param lb_min ...
82 : !> \param lb_max ...
83 : !> \param npgfb ...
84 : !> \param zetb ...
85 : !> \param rpgfb ...
86 : !> \param rb ...
87 : !> \param lc_min ...
88 : !> \param lc_max ...
89 : !> \param npgfc ...
90 : !> \param zetc ...
91 : !> \param rpgfc ...
92 : !> \param rc ...
93 : !> \param dab ...
94 : !> \param dac ...
95 : !> \param dbc ...
96 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
97 : !> \param potential_parameter the info about the potential
98 : !> \param int_abc_ext the extremal value of int_abc, i.e., MAXVAL(ABS(int_abc))
99 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
100 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
101 : !> the latter must be initialized too
102 : ! **************************************************************************************************
103 5343808 : SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
104 5343808 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
105 5343808 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
106 : dab, dac, dbc, lib, potential_parameter, &
107 : int_abc_ext)
108 :
109 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: int_abc
110 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
111 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
112 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
113 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
114 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
115 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
116 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
117 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
118 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
119 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
120 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
121 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
122 : REAL(dp), INTENT(INOUT), OPTIONAL :: int_abc_ext
123 :
124 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, ipgf, j, &
125 : jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
126 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
127 5343808 : REAL(dp), DIMENSION(:), POINTER :: p_work
128 : TYPE(params_3c), POINTER :: params
129 :
130 5343808 : NULLIFY (params, p_work)
131 160314240 : ALLOCATE (params)
132 :
133 5343808 : dr_ab = 0.0_dp
134 5343808 : dr_bc = 0.0_dp
135 5343808 : dr_ac = 0.0_dp
136 :
137 5343808 : op = potential_parameter%potential_type
138 :
139 : IF (op == do_potential_truncated .OR. op == do_potential_short &
140 5343808 : .OR. op == do_potential_mix_cl_trunc) THEN
141 4417257 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
142 4417257 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
143 926551 : ELSE IF (op == do_potential_coulomb) THEN
144 108645 : dr_bc = 1000000.0_dp
145 108645 : dr_ac = 1000000.0_dp
146 : END IF
147 :
148 5343808 : IF (PRESENT(int_abc_ext)) THEN
149 5217266 : int_abc_ext = 0.0_dp
150 : END IF
151 :
152 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
153 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
154 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
155 :
156 : !Looping over the pgfs
157 13512118 : DO ipgf = 1, npgfa
158 8168310 : zeti = zeta(ipgf)
159 8168310 : a_start = (ipgf - 1)*ncoset(la_max)
160 :
161 30599354 : DO jpgf = 1, npgfb
162 :
163 : ! screening
164 17087236 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
165 :
166 10007878 : zetj = zetb(jpgf)
167 10007878 : b_start = (jpgf - 1)*ncoset(lb_max)
168 :
169 41174802 : DO kpgf = 1, npgfc
170 :
171 : ! screening
172 22998614 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
173 15675571 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
174 :
175 13156074 : zetk = zetc(kpgf)
176 13156074 : c_start = (kpgf - 1)*ncoset(lc_max)
177 :
178 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
179 : CALL set_params_3c(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
180 13156074 : potential_parameter=potential_parameter, params_out=params)
181 :
182 29106913 : DO li = la_min, la_max
183 15950839 : a_offset = a_start + ncoset(li - 1)
184 15950839 : ncoa = nco(li)
185 44698311 : DO lj = MAX(li, lb_min), lb_max
186 15591398 : b_offset = b_start + ncoset(lj - 1)
187 15591398 : ncob = nco(lj)
188 56762758 : DO lk = lc_min, lc_max
189 25220521 : c_offset = c_start + ncoset(lk - 1)
190 25220521 : ncoc = nco(lk)
191 :
192 25220521 : a_mysize(1) = ncoa*ncob*ncoc
193 25220521 : CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
194 :
195 40811919 : IF (PRESENT(int_abc_ext)) THEN
196 97058043 : DO k = 1, ncoc
197 72942318 : p1 = (k - 1)*ncob
198 254096419 : DO j = 1, ncob
199 157038376 : p2 = (p1 + j - 1)*ncoa
200 482409506 : DO i = 1, ncoa
201 252428812 : p3 = p2 + i
202 252428812 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
203 409467188 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
204 : END DO
205 : END DO
206 : END DO
207 : ELSE
208 4968585 : DO k = 1, ncoc
209 3863789 : p1 = (k - 1)*ncob
210 13124670 : DO j = 1, ncob
211 8156085 : p2 = (p1 + j - 1)*ncoa
212 24731615 : DO i = 1, ncoa
213 12711741 : p3 = p2 + i
214 20867826 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
215 : END DO
216 : END DO
217 : END DO
218 : END IF
219 :
220 : END DO !lk
221 : END DO !lj
222 : END DO !li
223 :
224 : !swap centers 3 and 4 to compute (c|ab) with lb < la
225 13156074 : CALL set_params_3c(lib, rb, ra, rc, params_in=params)
226 :
227 46215307 : DO lj = lb_min, lb_max
228 15971997 : b_offset = b_start + ncoset(lj - 1)
229 15971997 : ncob = nco(lj)
230 43353061 : DO li = MAX(lj + 1, la_min), la_max
231 4382450 : a_offset = a_start + ncoset(li - 1)
232 4382450 : ncoa = nco(li)
233 28046893 : DO lk = lc_min, lc_max
234 7692446 : c_offset = c_start + ncoset(lk - 1)
235 7692446 : ncoc = nco(lk)
236 :
237 7692446 : a_mysize(1) = ncoa*ncob*ncoc
238 7692446 : CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
239 :
240 12074896 : IF (PRESENT(int_abc_ext)) THEN
241 30425178 : DO k = 1, ncoc
242 23058981 : p1 = (k - 1)*ncoa
243 116231263 : DO i = 1, ncoa
244 85806085 : p2 = (p1 + i - 1)*ncob
245 217307111 : DO j = 1, ncob
246 108442045 : p3 = p2 + j
247 108442045 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
248 194248130 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
249 : END DO
250 : END DO
251 : END DO
252 : ELSE
253 1482968 : DO k = 1, ncoc
254 1156719 : p1 = (k - 1)*ncoa
255 5851277 : DO i = 1, ncoa
256 4368309 : p2 = (p1 + i - 1)*ncob
257 11122929 : DO j = 1, ncob
258 5597901 : p3 = p2 + j
259 9966210 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
260 : END DO
261 : END DO
262 : END DO
263 : END IF
264 :
265 : END DO !lk
266 : END DO !li
267 : END DO !lj
268 :
269 : END DO !kpgf
270 : END DO !jpgf
271 : END DO !ipgf
272 :
273 5343808 : DEALLOCATE (params)
274 :
275 5343808 : END SUBROUTINE eri_3center
276 :
277 : ! **************************************************************************************************
278 : !> \brief Sets the internals of the cp_libint_t object for integrals of type (k|ji)
279 : !> \param lib ..
280 : !> \param ri ...
281 : !> \param rj ...
282 : !> \param rk ...
283 : !> \param zeti ...
284 : !> \param zetj ...
285 : !> \param zetk ...
286 : !> \param li_max ...
287 : !> \param lj_max ...
288 : !> \param lk_max ...
289 : !> \param potential_parameter ...
290 : !> \param params_in external parameters to use for libint
291 : !> \param params_out returns the libint parameters computed based on the other arguments
292 : !> \note The use of params_in and params_out comes from the fact that one might have to swap
293 : !> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
294 : !> remain the same upon such a change => might avoid recomputing things over and over again
295 : ! **************************************************************************************************
296 26312148 : SUBROUTINE set_params_3c(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
297 : potential_parameter, params_in, params_out)
298 :
299 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
300 : REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
301 : REAL(dp), INTENT(IN), OPTIONAL :: zeti, zetj, zetk
302 : INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
303 : TYPE(coulomb_operator_type), INTENT(IN), OPTIONAL :: potential_parameter
304 : TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
305 :
306 : INTEGER :: l
307 : LOGICAL :: use_gamma
308 : REAL(dp) :: gammaq, omega2, omega_corr, omega_corr2, &
309 : prefac, R, S1234, T, tmp
310 26312148 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
311 : TYPE(params_3c), POINTER :: params
312 :
313 : !Assume that one of params_in or params_out is present, and that in the latter case, all
314 : !other optional arguments are here
315 :
316 : !The internal structure of libint2 is based on 4-center integrals
317 : !For 3-center, one of those is a dummy center
318 : !The integral is assumed to be (k|ji) where the centers are ordered as:
319 : !k -> 1, j -> 3 and i -> 4 (the center #2 is the dummy center)
320 :
321 : !If external parameters are given, just use them
322 26312148 : IF (PRESENT(params_in)) THEN
323 13156074 : params => params_in
324 :
325 : !If no external parameters to use, compute them
326 : ELSE
327 13156074 : params => params_out
328 :
329 : !Note: some variable of 4-center integrals simplify with a dummy center:
330 : ! P -> rk, gammap -> zetk
331 13156074 : params%m_max = li_max + lj_max + lk_max
332 13156074 : gammaq = zeti + zetj
333 13156074 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
334 13156074 : params%ZetapEtaInv = 1._dp/(zetk + gammaq)
335 :
336 52624296 : params%Q = (zeti*ri + zetj*rj)*params%EtaInv
337 52624296 : params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
338 13156074 : params%Rho = zetk*gammaq/(zetk + gammaq)
339 :
340 289433628 : params%Fm = 0.0_dp
341 13156074 : SELECT CASE (potential_parameter%potential_type)
342 : CASE (do_potential_coulomb)
343 2670048 : T = params%Rho*SUM((params%Q - rk)**2)
344 2670048 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
345 667512 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
346 :
347 667512 : CALL fgamma(params%m_max, T, params%Fm)
348 14685264 : params%Fm = prefac*params%Fm
349 : CASE (do_potential_truncated)
350 7381671 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
351 29526684 : T = params%Rho*SUM((params%Q - rk)**2)
352 29526684 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
353 7381671 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
354 :
355 7381671 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
356 7381671 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
357 7381671 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
358 162396762 : params%Fm = prefac*params%Fm
359 : CASE (do_potential_short)
360 8439944 : T = params%Rho*SUM((params%Q - rk)**2)
361 8439944 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
362 2109986 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
363 :
364 2109986 : CALL fgamma(params%m_max, T, params%Fm)
365 :
366 2109986 : omega2 = potential_parameter%omega**2
367 2109986 : omega_corr2 = omega2/(omega2 + params%Rho)
368 2109986 : omega_corr = SQRT(omega_corr2)
369 2109986 : T = T*omega_corr2
370 2109986 : ALLOCATE (Fm(prim_data_f_size))
371 :
372 2109986 : CALL fgamma(params%m_max, T, Fm)
373 2109986 : tmp = -omega_corr
374 11707924 : DO l = 1, params%m_max + 1
375 9597938 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
376 11707924 : tmp = tmp*omega_corr2
377 : END DO
378 46419692 : params%Fm = prefac*params%Fm
379 : CASE (do_potential_mix_cl_trunc)
380 1399165 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
381 5596660 : T = params%Rho*SUM((params%Q - rk)**2)
382 5596660 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
383 1399165 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
384 :
385 1399165 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
386 1399165 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
387 1399165 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
388 :
389 1399165 : ALLOCATE (Fm(prim_data_f_size))
390 1399165 : CALL fgamma(params%m_max, T, Fm)
391 5377685 : DO l = 1, params%m_max + 1
392 : params%Fm(l) = params%Fm(l) &
393 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
394 5377685 : - Fm(l)*potential_parameter%scale_longrange
395 : END DO
396 1399165 : DEALLOCATE (Fm)
397 :
398 1399165 : omega2 = potential_parameter%omega**2
399 1399165 : omega_corr2 = omega2/(omega2 + params%Rho)
400 1399165 : omega_corr = SQRT(omega_corr2)
401 1399165 : T = T*omega_corr2
402 :
403 1399165 : ALLOCATE (Fm(prim_data_f_size))
404 1399165 : CALL fgamma(params%m_max, T, Fm)
405 1399165 : tmp = omega_corr
406 5377685 : DO l = 1, params%m_max + 1
407 3978520 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
408 5377685 : tmp = tmp*omega_corr2
409 : END DO
410 30781630 : params%Fm = prefac*params%Fm
411 : CASE (do_potential_id)
412 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
413 11184180 : - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
414 1597740 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
415 :
416 35150280 : params%Fm(:) = prefac
417 : CASE DEFAULT
418 13156074 : CPABORT("Requested operator NYI")
419 : END SELECT
420 :
421 : END IF
422 :
423 : CALL cp_libint_set_params_eri(lib, rk, rk, rj, ri, params%ZetaInv, params%EtaInv, &
424 : params%ZetapEtaInv, params%Rho, rk, params%Q, params%W, &
425 26312148 : params%m_max, params%Fm)
426 :
427 26312148 : END SUBROUTINE set_params_3c
428 :
429 : ! **************************************************************************************************
430 : !> \brief Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given
431 : !> set of cartesian gaussian orbitals. Returns x,y,z derivatives for 1st and 2nd center
432 : !> \param der_abc_1 the derivatives for the 1st center (allocated before hand)
433 : !> \param der_abc_2 the derivatives for the 2nd center (allocated before hand)
434 : !> \param la_min ...
435 : !> \param la_max ...
436 : !> \param npgfa ...
437 : !> \param zeta ...
438 : !> \param rpgfa ...
439 : !> \param ra ...
440 : !> \param lb_min ...
441 : !> \param lb_max ...
442 : !> \param npgfb ...
443 : !> \param zetb ...
444 : !> \param rpgfb ...
445 : !> \param rb ...
446 : !> \param lc_min ...
447 : !> \param lc_max ...
448 : !> \param npgfc ...
449 : !> \param zetc ...
450 : !> \param rpgfc ...
451 : !> \param rc ...
452 : !> \param dab ...
453 : !> \param dac ...
454 : !> \param dbc ...
455 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
456 : !> \param potential_parameter the info about the potential
457 : !> \param der_abc_1_ext the extremal value of der_abc_1, i.e., MAXVAL(ABS(der_abc_1))
458 : !> \param der_abc_2_ext ...
459 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
460 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
461 : !> the latter must be initialized too. Note that the derivative wrt to the third center
462 : !> can be obtained via translational invariance
463 : ! **************************************************************************************************
464 345768 : SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
465 345768 : la_min, la_max, npgfa, zeta, rpgfa, ra, &
466 345768 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
467 345768 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
468 : dab, dac, dbc, lib, potential_parameter, &
469 : der_abc_1_ext, der_abc_2_ext)
470 :
471 : REAL(dp), DIMENSION(:, :, :, :), INTENT(INOUT) :: der_abc_1, der_abc_2
472 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
473 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
474 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
475 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
476 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
477 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
478 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
479 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
480 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
481 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
482 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
483 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
484 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: der_abc_1_ext, der_abc_2_ext
485 :
486 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, i_deriv, &
487 : ipgf, j, jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
488 : INTEGER, DIMENSION(3) :: permute_1, permute_2
489 : LOGICAL :: do_ext
490 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
491 : REAL(dp), DIMENSION(3) :: der_abc_1_ext_prv, der_abc_2_ext_prv
492 345768 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
493 : TYPE(params_3c), POINTER :: params
494 :
495 345768 : NULLIFY (params, p_deriv)
496 10373040 : ALLOCATE (params)
497 :
498 345768 : permute_1 = [4, 5, 6]
499 345768 : permute_2 = [7, 8, 9]
500 :
501 345768 : dr_ab = 0.0_dp
502 345768 : dr_bc = 0.0_dp
503 345768 : dr_ac = 0.0_dp
504 :
505 345768 : op = potential_parameter%potential_type
506 :
507 : IF (op == do_potential_truncated .OR. op == do_potential_short &
508 345768 : .OR. op == do_potential_mix_cl_trunc) THEN
509 119807 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
510 119807 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
511 225961 : ELSE IF (op == do_potential_coulomb) THEN
512 9988 : dr_bc = 1000000.0_dp
513 9988 : dr_ac = 1000000.0_dp
514 : END IF
515 :
516 345768 : do_ext = .FALSE.
517 345768 : IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
518 345768 : der_abc_1_ext_prv = 0.0_dp
519 345768 : der_abc_2_ext_prv = 0.0_dp
520 :
521 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
522 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
523 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
524 :
525 : !Looping over the pgfs
526 1188953 : DO ipgf = 1, npgfa
527 843185 : zeti = zeta(ipgf)
528 843185 : a_start = (ipgf - 1)*ncoset(la_max)
529 :
530 4761095 : DO jpgf = 1, npgfb
531 :
532 : ! screening
533 3572142 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
534 :
535 1134668 : zetj = zetb(jpgf)
536 1134668 : b_start = (jpgf - 1)*ncoset(lb_max)
537 :
538 8321733 : DO kpgf = 1, npgfc
539 :
540 : ! screening
541 6343880 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
542 2667010 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
543 :
544 1871344 : zetk = zetc(kpgf)
545 1871344 : c_start = (kpgf - 1)*ncoset(lc_max)
546 :
547 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
548 : CALL set_params_3c_deriv(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
549 1871344 : potential_parameter=potential_parameter, params_out=params)
550 :
551 4308805 : DO li = la_min, la_max
552 2437461 : a_offset = a_start + ncoset(li - 1)
553 2437461 : ncoa = nco(li)
554 6910844 : DO lj = MAX(li, lb_min), lb_max
555 2602039 : b_offset = b_start + ncoset(lj - 1)
556 2602039 : ncob = nco(lj)
557 8804159 : DO lk = lc_min, lc_max
558 3764659 : c_offset = c_start + ncoset(lk - 1)
559 3764659 : ncoc = nco(lk)
560 :
561 3764659 : a_mysize(1) = ncoa*ncob*ncoc
562 :
563 3764659 : CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
564 :
565 3764659 : IF (do_ext) THEN
566 15058636 : DO i_deriv = 1, 3
567 42682984 : DO k = 1, ncoc
568 27624348 : p1 = (k - 1)*ncob
569 90934137 : DO j = 1, ncob
570 52015812 : p2 = (p1 + j - 1)*ncoa
571 158121102 : DO i = 1, ncoa
572 78480942 : p3 = p2 + i
573 :
574 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
575 78480942 : p_deriv(p3, permute_2(i_deriv))
576 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
577 78480942 : ABS(p_deriv(p3, permute_2(i_deriv))))
578 :
579 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
580 78480942 : p_deriv(p3, permute_1(i_deriv))
581 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
582 130496754 : ABS(p_deriv(p3, permute_1(i_deriv))))
583 :
584 : END DO
585 : END DO
586 : END DO
587 : END DO
588 : ELSE
589 0 : DO i_deriv = 1, 3
590 0 : DO k = 1, ncoc
591 0 : p1 = (k - 1)*ncob
592 0 : DO j = 1, ncob
593 0 : p2 = (p1 + j - 1)*ncoa
594 0 : DO i = 1, ncoa
595 0 : p3 = p2 + i
596 :
597 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
598 0 : p_deriv(p3, permute_2(i_deriv))
599 :
600 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
601 0 : p_deriv(p3, permute_1(i_deriv))
602 : END DO
603 : END DO
604 : END DO
605 : END DO
606 : END IF
607 :
608 6366698 : DEALLOCATE (p_deriv)
609 : END DO !lk
610 : END DO !lj
611 : END DO !li
612 :
613 : !swap centers 3 and 4 to compute (c|ab) with lb < la
614 1871344 : CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
615 :
616 7890144 : DO lj = lb_min, lb_max
617 2446658 : b_offset = b_start + ncoset(lj - 1)
618 2446658 : ncob = nco(lj)
619 9454558 : DO li = MAX(lj + 1, la_min), la_max
620 664020 : a_offset = a_start + ncoset(li - 1)
621 664020 : ncoa = nco(li)
622 4083837 : DO lk = lc_min, lc_max
623 973159 : c_offset = c_start + ncoset(lk - 1)
624 973159 : ncoc = nco(lk)
625 :
626 973159 : a_mysize(1) = ncoa*ncob*ncoc
627 973159 : CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
628 :
629 973159 : IF (do_ext) THEN
630 3892636 : DO i_deriv = 1, 3
631 11260180 : DO k = 1, ncoc
632 7367544 : p1 = (k - 1)*ncoa
633 34472487 : DO i = 1, ncoa
634 24185466 : p2 = (p1 + i - 1)*ncob
635 59226300 : DO j = 1, ncob
636 27673290 : p3 = p2 + j
637 :
638 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
639 27673290 : p_deriv(p3, permute_1(i_deriv))
640 :
641 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
642 27673290 : ABS(p_deriv(p3, permute_1(i_deriv))))
643 :
644 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
645 27673290 : p_deriv(p3, permute_2(i_deriv))
646 :
647 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
648 51858756 : ABS(p_deriv(p3, permute_2(i_deriv))))
649 : END DO
650 : END DO
651 : END DO
652 : END DO
653 : ELSE
654 0 : DO i_deriv = 1, 3
655 0 : DO k = 1, ncoc
656 0 : p1 = (k - 1)*ncoa
657 0 : DO i = 1, ncoa
658 0 : p2 = (p1 + i - 1)*ncob
659 0 : DO j = 1, ncob
660 0 : p3 = p2 + j
661 :
662 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
663 0 : p_deriv(p3, permute_1(i_deriv))
664 :
665 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
666 0 : p_deriv(p3, permute_2(i_deriv))
667 : END DO
668 : END DO
669 : END DO
670 : END DO
671 : END IF
672 :
673 1637179 : DEALLOCATE (p_deriv)
674 : END DO !lk
675 : END DO !li
676 : END DO !lj
677 :
678 : END DO !kpgf
679 : END DO !jpgf
680 : END DO !ipgf
681 :
682 345768 : IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
683 345768 : IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
684 :
685 345768 : DEALLOCATE (params)
686 :
687 345768 : END SUBROUTINE eri_3center_derivs
688 :
689 : ! **************************************************************************************************
690 : !> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|ji)
691 : !> \param lib ..
692 : !> \param ri ...
693 : !> \param rj ...
694 : !> \param rk ...
695 : !> \param zeti ...
696 : !> \param zetj ...
697 : !> \param zetk ...
698 : !> \param li_max ...
699 : !> \param lj_max ...
700 : !> \param lk_max ...
701 : !> \param potential_parameter ...
702 : !> \param params_in ...
703 : !> \param params_out ...
704 : !> \note The use of params_in and params_out comes from the fact that one might have to swap
705 : !> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
706 : !> remain the same upon such a change => might avoid recomputing things over and over again
707 : ! **************************************************************************************************
708 3742688 : SUBROUTINE set_params_3c_deriv(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
709 : potential_parameter, params_in, params_out)
710 :
711 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
712 : REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
713 : REAL(dp), INTENT(IN) :: zeti, zetj, zetk
714 : INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
715 : TYPE(coulomb_operator_type), INTENT(IN), OPTIONAL :: potential_parameter
716 : TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
717 :
718 : INTEGER :: l
719 : LOGICAL :: use_gamma
720 : REAL(dp) :: gammaq, omega2, omega_corr, omega_corr2, &
721 : prefac, R, S1234, T, tmp
722 3742688 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
723 : TYPE(params_3c), POINTER :: params
724 :
725 3742688 : IF (PRESENT(params_in)) THEN
726 1871344 : params => params_in
727 :
728 : ELSE
729 1871344 : params => params_out
730 :
731 1871344 : params%m_max = li_max + lj_max + lk_max + 1
732 1871344 : gammaq = zeti + zetj
733 1871344 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
734 1871344 : params%ZetapEtaInv = 1._dp/(zetk + gammaq)
735 :
736 7485376 : params%Q = (zeti*ri + zetj*rj)*params%EtaInv
737 7485376 : params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
738 1871344 : params%Rho = zetk*gammaq/(zetk + gammaq)
739 :
740 41169568 : params%Fm = 0.0_dp
741 1871344 : SELECT CASE (potential_parameter%potential_type)
742 : CASE (do_potential_coulomb)
743 416424 : T = params%Rho*SUM((params%Q - rk)**2)
744 416424 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
745 104106 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
746 :
747 104106 : CALL fgamma(params%m_max, T, params%Fm)
748 2290332 : params%Fm = prefac*params%Fm
749 : CASE (do_potential_truncated)
750 1244984 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
751 4979936 : T = params%Rho*SUM((params%Q - rk)**2)
752 4979936 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
753 1244984 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
754 :
755 1244984 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
756 1244984 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
757 1244984 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
758 27389648 : params%Fm = prefac*params%Fm
759 : CASE (do_potential_short)
760 0 : T = params%Rho*SUM((params%Q - rk)**2)
761 0 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
762 0 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
763 :
764 0 : CALL fgamma(params%m_max, T, params%Fm)
765 :
766 0 : omega2 = potential_parameter%omega**2
767 0 : omega_corr2 = omega2/(omega2 + params%Rho)
768 0 : omega_corr = SQRT(omega_corr2)
769 0 : T = T*omega_corr2
770 0 : ALLOCATE (Fm(prim_data_f_size))
771 :
772 0 : CALL fgamma(params%m_max, T, Fm)
773 0 : tmp = -omega_corr
774 0 : DO l = 1, params%m_max + 1
775 0 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
776 0 : tmp = tmp*omega_corr2
777 : END DO
778 0 : params%Fm = prefac*params%Fm
779 : CASE (do_potential_mix_cl_trunc)
780 0 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
781 0 : T = params%Rho*SUM((params%Q - rk)**2)
782 0 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
783 0 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
784 :
785 0 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
786 0 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
787 0 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
788 :
789 0 : ALLOCATE (Fm(prim_data_f_size))
790 0 : CALL fgamma(params%m_max, T, Fm)
791 0 : DO l = 1, params%m_max + 1
792 : params%Fm(l) = params%Fm(l) &
793 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
794 0 : - Fm(l)*potential_parameter%scale_longrange
795 : END DO
796 0 : DEALLOCATE (Fm)
797 :
798 0 : omega2 = potential_parameter%omega**2
799 0 : omega_corr2 = omega2/(omega2 + params%Rho)
800 0 : omega_corr = SQRT(omega_corr2)
801 0 : T = T*omega_corr2
802 :
803 0 : ALLOCATE (Fm(prim_data_f_size))
804 0 : CALL fgamma(params%m_max, T, Fm)
805 0 : tmp = omega_corr
806 0 : DO l = 1, params%m_max + 1
807 0 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
808 0 : tmp = tmp*omega_corr2
809 : END DO
810 0 : params%Fm = prefac*params%Fm
811 : CASE (do_potential_id)
812 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
813 3655778 : - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
814 522254 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
815 :
816 11489588 : params%Fm(:) = prefac
817 : CASE DEFAULT
818 1871344 : CPABORT("Requested operator NYI")
819 : END SELECT
820 :
821 : END IF
822 :
823 : CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, ri, rk, &
824 : params%Q, params%W, zetk, 0.0_dp, zetj, zeti, params%ZetaInv, &
825 3742688 : params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
826 :
827 3742688 : END SUBROUTINE set_params_3c_deriv
828 :
829 : ! **************************************************************************************************
830 : !> \brief Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian
831 : !> gaussian orbitals
832 : !> \param int_ab the integrals as array of cartesian orbitals (allocated before hand)
833 : !> \param la_min ...
834 : !> \param la_max ...
835 : !> \param npgfa ...
836 : !> \param zeta ...
837 : !> \param rpgfa ...
838 : !> \param ra ...
839 : !> \param lb_min ...
840 : !> \param lb_max ...
841 : !> \param npgfb ...
842 : !> \param zetb ...
843 : !> \param rpgfb ...
844 : !> \param rb ...
845 : !> \param dab ...
846 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
847 : !> \param potential_parameter the info about the potential
848 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
849 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
850 : !> the latter must be initialized too
851 : ! **************************************************************************************************
852 552655 : SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
853 552655 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
854 : dab, lib, potential_parameter)
855 :
856 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int_ab
857 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
858 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
859 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
860 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
861 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
862 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
863 : REAL(dp), INTENT(IN) :: dab
864 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
865 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
866 :
867 : INTEGER :: a_mysize(1), a_offset, a_start, &
868 : b_offset, b_start, i, ipgf, j, jpgf, &
869 : li, lj, ncoa, ncob, p1, p2
870 : REAL(dp) :: dr_ab, zeti, zetj
871 552655 : REAL(dp), DIMENSION(:), POINTER :: p_work
872 :
873 552655 : NULLIFY (p_work)
874 :
875 552655 : dr_ab = 0.0_dp
876 :
877 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
878 282343 : potential_parameter%potential_type == do_potential_short .OR. &
879 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
880 282948 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
881 269707 : ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
882 93986 : dr_ab = 1000000.0_dp
883 : END IF
884 :
885 : !Looping over the pgfs
886 1604908 : DO ipgf = 1, npgfa
887 1052253 : zeti = zeta(ipgf)
888 1052253 : a_start = (ipgf - 1)*ncoset(la_max)
889 :
890 6212564 : DO jpgf = 1, npgfb
891 4607656 : zetj = zetb(jpgf)
892 4607656 : b_start = (jpgf - 1)*ncoset(lb_max)
893 :
894 : !screening
895 4607656 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
896 :
897 3011872 : CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
898 :
899 9525013 : DO li = la_min, la_max
900 5460888 : a_offset = a_start + ncoset(li - 1)
901 5460888 : ncoa = nco(li)
902 20734363 : DO lj = lb_min, lb_max
903 10665819 : b_offset = b_start + ncoset(lj - 1)
904 10665819 : ncob = nco(lj)
905 :
906 10665819 : a_mysize(1) = ncoa*ncob
907 10665819 : CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
908 :
909 48702366 : DO j = 1, ncob
910 32575659 : p1 = (j - 1)*ncoa
911 144569598 : DO i = 1, ncoa
912 101328120 : p2 = p1 + i
913 133903779 : int_ab(a_offset + i, b_offset + j) = p_work(p2)
914 : END DO
915 : END DO
916 :
917 : END DO
918 : END DO
919 :
920 : END DO
921 : END DO
922 :
923 552655 : END SUBROUTINE eri_2center
924 :
925 : ! **************************************************************************************************
926 : !> \brief Sets the internals of the cp_libint_t object for integrals of type (k|j)
927 : !> \param lib ..
928 : !> \param rj ...
929 : !> \param rk ...
930 : !> \param zetj ...
931 : !> \param zetk ...
932 : !> \param lj_max ...
933 : !> \param lk_max ...
934 : !> \param potential_parameter ...
935 : ! **************************************************************************************************
936 3011872 : SUBROUTINE set_params_2c(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
937 :
938 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
939 : REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
940 : REAL(dp), INTENT(IN) :: zetj, zetk
941 : INTEGER, INTENT(IN) :: lj_max, lk_max
942 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
943 :
944 : INTEGER :: l, op
945 : LOGICAL :: use_gamma
946 : REAL(dp) :: omega2, omega_corr, omega_corr2, prefac, &
947 : R, T, tmp
948 3011872 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
949 : TYPE(params_2c) :: params
950 :
951 : !The internal structure of libint2 is based on 4-center integrals
952 : !For 2-center, two of those are dummy centers
953 : !The integral is assumed to be (k|j) where the centers are ordered as:
954 : !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)
955 :
956 : !Note: some variable of 4-center integrals simplify due to dummy centers:
957 : ! P -> rk, gammap -> zetk
958 : ! Q -> rj, gammaq -> zetj
959 :
960 3011872 : op = potential_parameter%potential_type
961 3011872 : params%m_max = lj_max + lk_max
962 3011872 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
963 3011872 : params%ZetapEtaInv = 1._dp/(zetk + zetj)
964 :
965 12047488 : params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
966 3011872 : params%Rho = zetk*zetj/(zetk + zetj)
967 :
968 66261184 : params%Fm = 0.0_dp
969 : SELECT CASE (op)
970 : CASE (do_potential_coulomb)
971 482472 : T = params%Rho*SUM((rj - rk)**2)
972 120618 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
973 120618 : CALL fgamma(params%m_max, T, params%Fm)
974 2653596 : params%Fm = prefac*params%Fm
975 : CASE (do_potential_truncated)
976 220572 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
977 882288 : T = params%Rho*SUM((rj - rk)**2)
978 220572 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
979 :
980 220572 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
981 220572 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
982 220572 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
983 4852584 : params%Fm = prefac*params%Fm
984 : CASE (do_potential_short)
985 10017192 : T = params%Rho*SUM((rj - rk)**2)
986 2504298 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
987 :
988 2504298 : CALL fgamma(params%m_max, T, params%Fm)
989 :
990 2504298 : omega2 = potential_parameter%omega**2
991 2504298 : omega_corr2 = omega2/(omega2 + params%Rho)
992 2504298 : omega_corr = SQRT(omega_corr2)
993 2504298 : T = T*omega_corr2
994 2504298 : ALLOCATE (Fm(prim_data_f_size))
995 :
996 2504298 : CALL fgamma(params%m_max, T, Fm)
997 2504298 : tmp = -omega_corr
998 12305810 : DO l = 1, params%m_max + 1
999 9801512 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
1000 12305810 : tmp = tmp*omega_corr2
1001 : END DO
1002 55094556 : params%Fm = prefac*params%Fm
1003 : CASE (do_potential_mix_cl_trunc)
1004 61464 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
1005 245856 : T = params%Rho*SUM((rj - rk)**2)
1006 61464 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1007 :
1008 61464 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
1009 61464 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
1010 61464 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
1011 :
1012 61464 : ALLOCATE (Fm(prim_data_f_size))
1013 61464 : CALL fgamma(params%m_max, T, Fm)
1014 223681 : DO l = 1, params%m_max + 1
1015 : params%Fm(l) = params%Fm(l) &
1016 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
1017 223681 : - Fm(l)*potential_parameter%scale_longrange
1018 : END DO
1019 61464 : DEALLOCATE (Fm)
1020 :
1021 61464 : omega2 = potential_parameter%omega**2
1022 61464 : omega_corr2 = omega2/(omega2 + params%Rho)
1023 61464 : omega_corr = SQRT(omega_corr2)
1024 61464 : T = T*omega_corr2
1025 :
1026 61464 : ALLOCATE (Fm(prim_data_f_size))
1027 61464 : CALL fgamma(params%m_max, T, Fm)
1028 61464 : tmp = omega_corr
1029 223681 : DO l = 1, params%m_max + 1
1030 162217 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
1031 223681 : tmp = tmp*omega_corr2
1032 : END DO
1033 1352208 : params%Fm = prefac*params%Fm
1034 : CASE (do_potential_id)
1035 :
1036 419680 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
1037 2308240 : params%Fm(:) = prefac
1038 : CASE DEFAULT
1039 3011872 : CPABORT("Requested operator NYI")
1040 : END SELECT
1041 :
1042 : CALL cp_libint_set_params_eri(lib, rk, rk, rj, rj, params%ZetaInv, params%EtaInv, &
1043 : params%ZetapEtaInv, params%Rho, rk, rj, params%W, &
1044 3011872 : params%m_max, params%Fm)
1045 :
1046 78308672 : END SUBROUTINE set_params_2c
1047 :
1048 : ! **************************************************************************************************
1049 : !> \brief Helper function to compare Coulomb operator types
1050 : !> \param potential1 first potential
1051 : !> \param potential2 second potential
1052 : !> \return Boolean whether both potentials are equal
1053 : ! **************************************************************************************************
1054 14400 : PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
1055 : TYPE(coulomb_operator_type), INTENT(IN) :: potential1, potential2
1056 : LOGICAL :: equals
1057 :
1058 14400 : IF (potential1%potential_type /= potential2%potential_type) THEN
1059 : equals = .FALSE.
1060 : ELSE
1061 13572 : equals = .TRUE.
1062 736 : SELECT CASE (potential1%potential_type)
1063 : CASE (do_potential_short, do_potential_long)
1064 736 : IF (potential1%omega /= potential2%omega) equals = .FALSE.
1065 : CASE (do_potential_truncated)
1066 14 : IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
1067 : CASE (do_potential_mix_cl_trunc)
1068 2 : IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
1069 2 : IF (potential1%omega /= potential2%omega) equals = .FALSE.
1070 2 : IF (potential1%scale_coulomb /= potential2%scale_coulomb) equals = .FALSE.
1071 13574 : IF (potential1%scale_longrange /= potential2%scale_longrange) equals = .FALSE.
1072 : END SELECT
1073 : END IF
1074 :
1075 14400 : END FUNCTION compare_potential_types
1076 :
1077 : !> \brief Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given
1078 : !> set of cartesian gaussian orbitals. Returns the derivatives wrt to the first center
1079 : !> \param der_ab the derivatives as array of cartesian orbitals (allocated before hand)
1080 : !> \param la_min ...
1081 : !> \param la_max ...
1082 : !> \param npgfa ...
1083 : !> \param zeta ...
1084 : !> \param rpgfa ...
1085 : !> \param ra ...
1086 : !> \param lb_min ...
1087 : !> \param lb_max ...
1088 : !> \param npgfb ...
1089 : !> \param zetb ...
1090 : !> \param rpgfb ...
1091 : !> \param rb ...
1092 : !> \param dab ...
1093 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
1094 : !> \param potential_parameter the info about the potential
1095 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
1096 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
1097 : !> the latter must be initialized too
1098 : ! **************************************************************************************************
1099 144023 : SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
1100 144023 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
1101 : dab, lib, potential_parameter)
1102 :
1103 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: der_ab
1104 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
1105 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1106 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
1107 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
1108 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1109 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
1110 : REAL(dp), INTENT(IN) :: dab
1111 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
1112 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
1113 :
1114 : INTEGER :: a_mysize(1), a_offset, a_start, &
1115 : b_offset, b_start, i, i_deriv, ipgf, &
1116 : j, jpgf, li, lj, ncoa, ncob, p1, p2
1117 : INTEGER, DIMENSION(3) :: permute
1118 : REAL(dp) :: dr_ab, zeti, zetj
1119 144023 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
1120 :
1121 144023 : NULLIFY (p_deriv)
1122 :
1123 144023 : permute = [4, 5, 6]
1124 :
1125 144023 : dr_ab = 0.0_dp
1126 :
1127 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
1128 80255 : potential_parameter%potential_type == do_potential_short .OR. &
1129 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
1130 64460 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
1131 79563 : ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
1132 9831 : dr_ab = 1000000.0_dp
1133 : END IF
1134 :
1135 : !Looping over the pgfs
1136 640202 : DO ipgf = 1, npgfa
1137 496179 : zeti = zeta(ipgf)
1138 496179 : a_start = (ipgf - 1)*ncoset(la_max)
1139 :
1140 3871477 : DO jpgf = 1, npgfb
1141 3231275 : zetj = zetb(jpgf)
1142 3231275 : b_start = (jpgf - 1)*ncoset(lb_max)
1143 :
1144 : !screening
1145 3231275 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
1146 :
1147 2175084 : CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
1148 :
1149 6319188 : DO li = la_min, la_max
1150 3647925 : a_offset = a_start + ncoset(li - 1)
1151 3647925 : ncoa = nco(li)
1152 13047347 : DO lj = lb_min, lb_max
1153 6168147 : b_offset = b_start + ncoset(lj - 1)
1154 6168147 : ncob = nco(lj)
1155 :
1156 6168147 : a_mysize(1) = ncoa*ncob
1157 6168147 : CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
1158 :
1159 24672588 : DO i_deriv = 1, 3
1160 75351867 : DO j = 1, ncob
1161 50679279 : p1 = (j - 1)*ncoa
1162 208445424 : DO i = 1, ncoa
1163 139261704 : p2 = p1 + i
1164 189940983 : der_ab(a_offset + i, b_offset + j, i_deriv) = p_deriv(p2, permute(i_deriv))
1165 : END DO
1166 : END DO
1167 : END DO
1168 :
1169 9816072 : DEALLOCATE (p_deriv)
1170 : END DO
1171 : END DO
1172 :
1173 : END DO
1174 : END DO
1175 :
1176 144023 : END SUBROUTINE eri_2center_derivs
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|j)
1180 : !> \param lib ..
1181 : !> \param rj ...
1182 : !> \param rk ...
1183 : !> \param zetj ...
1184 : !> \param zetk ...
1185 : !> \param lj_max ...
1186 : !> \param lk_max ...
1187 : !> \param potential_parameter ...
1188 : ! **************************************************************************************************
1189 2175084 : SUBROUTINE set_params_2c_deriv(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
1190 :
1191 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
1192 : REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
1193 : REAL(dp), INTENT(IN) :: zetj, zetk
1194 : INTEGER, INTENT(IN) :: lj_max, lk_max
1195 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
1196 :
1197 : INTEGER :: l, op
1198 : LOGICAL :: use_gamma
1199 : REAL(dp) :: omega2, omega_corr, omega_corr2, prefac, &
1200 : R, T, tmp
1201 2175084 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
1202 : TYPE(params_2c) :: params
1203 :
1204 : !The internal structure of libint2 is based on 4-center integrals
1205 : !For 2-center, two of those are dummy centers
1206 : !The integral is assumed to be (k|j) where the centers are ordered as:
1207 : !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)
1208 :
1209 : !Note: some variable of 4-center integrals simplify due to dummy centers:
1210 : ! P -> rk, gammap -> zetk
1211 : ! Q -> rj, gammaq -> zetj
1212 :
1213 2175084 : op = potential_parameter%potential_type
1214 2175084 : params%m_max = lj_max + lk_max + 1
1215 2175084 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
1216 2175084 : params%ZetapEtaInv = 1._dp/(zetk + zetj)
1217 :
1218 8700336 : params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
1219 2175084 : params%Rho = zetk*zetj/(zetk + zetj)
1220 :
1221 47851848 : params%Fm = 0.0_dp
1222 : SELECT CASE (op)
1223 : CASE (do_potential_coulomb)
1224 102916 : T = params%Rho*SUM((rj - rk)**2)
1225 25729 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1226 25729 : CALL fgamma(params%m_max, T, params%Fm)
1227 566038 : params%Fm = prefac*params%Fm
1228 : CASE (do_potential_truncated)
1229 61941 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
1230 247764 : T = params%Rho*SUM((rj - rk)**2)
1231 61941 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1232 :
1233 61941 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
1234 61941 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
1235 61941 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
1236 1362702 : params%Fm = prefac*params%Fm
1237 : CASE (do_potential_short)
1238 8096600 : T = params%Rho*SUM((rj - rk)**2)
1239 2024150 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1240 :
1241 2024150 : CALL fgamma(params%m_max, T, params%Fm)
1242 :
1243 2024150 : omega2 = potential_parameter%omega**2
1244 2024150 : omega_corr2 = omega2/(omega2 + params%Rho)
1245 2024150 : omega_corr = SQRT(omega_corr2)
1246 2024150 : T = T*omega_corr2
1247 2024150 : ALLOCATE (Fm(prim_data_f_size))
1248 :
1249 2024150 : CALL fgamma(params%m_max, T, Fm)
1250 2024150 : tmp = -omega_corr
1251 11369176 : DO l = 1, params%m_max + 1
1252 9345026 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
1253 11369176 : tmp = tmp*omega_corr2
1254 : END DO
1255 44531300 : params%Fm = prefac*params%Fm
1256 : CASE (do_potential_mix_cl_trunc)
1257 2578 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
1258 10312 : T = params%Rho*SUM((rj - rk)**2)
1259 2578 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1260 :
1261 2578 : CPASSERT(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
1262 2578 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
1263 2578 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
1264 :
1265 2578 : ALLOCATE (Fm(prim_data_f_size))
1266 2578 : CALL fgamma(params%m_max, T, Fm)
1267 8234 : DO l = 1, params%m_max + 1
1268 : params%Fm(l) = params%Fm(l) &
1269 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
1270 8234 : - Fm(l)*potential_parameter%scale_longrange
1271 : END DO
1272 2578 : DEALLOCATE (Fm)
1273 :
1274 2578 : omega2 = potential_parameter%omega**2
1275 2578 : omega_corr2 = omega2/(omega2 + params%Rho)
1276 2578 : omega_corr = SQRT(omega_corr2)
1277 2578 : T = T*omega_corr2
1278 :
1279 2578 : ALLOCATE (Fm(prim_data_f_size))
1280 2578 : CALL fgamma(params%m_max, T, Fm)
1281 2578 : tmp = omega_corr
1282 8234 : DO l = 1, params%m_max + 1
1283 5656 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
1284 8234 : tmp = tmp*omega_corr2
1285 : END DO
1286 56716 : params%Fm = prefac*params%Fm
1287 : CASE (do_potential_id)
1288 :
1289 242744 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
1290 1335092 : params%Fm(:) = prefac
1291 : CASE DEFAULT
1292 2175084 : CPABORT("Requested operator NYI")
1293 : END SELECT
1294 :
1295 : CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, rj, rk, rj, params%W, zetk, 0.0_dp, &
1296 : zetj, 0.0_dp, params%ZetaInv, params%EtaInv, &
1297 : params%ZetapEtaInv, params%Rho, &
1298 2175084 : params%m_max, params%Fm)
1299 :
1300 56552184 : END SUBROUTINE set_params_2c_deriv
1301 :
1302 0 : END MODULE libint_2c_3c
|
