(git:24d69ee)
Loading...
Searching...
No Matches
libint_2c_3c.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 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
14MODULE 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
69CONTAINS
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 SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
104 lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
105 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 REAL(dp), DIMENSION(:), POINTER :: p_work
128 TYPE(params_3c), POINTER :: params
129
130 NULLIFY (params, p_work)
131 ALLOCATE (params)
132
133 dr_ab = 0.0_dp
134 dr_bc = 0.0_dp
135 dr_ac = 0.0_dp
136
137 op = potential_parameter%potential_type
138
139 IF (op == do_potential_truncated .OR. op == do_potential_short &
140 .OR. op == do_potential_mix_cl_trunc) THEN
141 dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
142 dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
143 ELSE IF (op == do_potential_coulomb) THEN
144 dr_bc = 1000000.0_dp
145 dr_ac = 1000000.0_dp
146 END IF
147
148 IF (PRESENT(int_abc_ext)) THEN
149 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 DO ipgf = 1, npgfa
158 zeti = zeta(ipgf)
159 a_start = (ipgf - 1)*ncoset(la_max)
160
161 DO jpgf = 1, npgfb
162
163 ! screening
164 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
165
166 zetj = zetb(jpgf)
167 b_start = (jpgf - 1)*ncoset(lb_max)
168
169 DO kpgf = 1, npgfc
170
171 ! screening
172 IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) cycle
173 IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) cycle
174
175 zetk = zetc(kpgf)
176 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 potential_parameter=potential_parameter, params_out=params)
181
182 DO li = la_min, la_max
183 a_offset = a_start + ncoset(li - 1)
184 ncoa = nco(li)
185 DO lj = max(li, lb_min), lb_max
186 b_offset = b_start + ncoset(lj - 1)
187 ncob = nco(lj)
188 DO lk = lc_min, lc_max
189 c_offset = c_start + ncoset(lk - 1)
190 ncoc = nco(lk)
191
192 a_mysize(1) = ncoa*ncob*ncoc
193 CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
194
195 IF (PRESENT(int_abc_ext)) THEN
196 DO k = 1, ncoc
197 p1 = (k - 1)*ncob
198 DO j = 1, ncob
199 p2 = (p1 + j - 1)*ncoa
200 DO i = 1, ncoa
201 p3 = p2 + i
202 int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
203 int_abc_ext = max(int_abc_ext, abs(p_work(p3)))
204 END DO
205 END DO
206 END DO
207 ELSE
208 DO k = 1, ncoc
209 p1 = (k - 1)*ncob
210 DO j = 1, ncob
211 p2 = (p1 + j - 1)*ncoa
212 DO i = 1, ncoa
213 p3 = p2 + i
214 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 CALL set_params_3c(lib, rb, ra, rc, params_in=params)
226
227 DO lj = lb_min, lb_max
228 b_offset = b_start + ncoset(lj - 1)
229 ncob = nco(lj)
230 DO li = max(lj + 1, la_min), la_max
231 a_offset = a_start + ncoset(li - 1)
232 ncoa = nco(li)
233 DO lk = lc_min, lc_max
234 c_offset = c_start + ncoset(lk - 1)
235 ncoc = nco(lk)
236
237 a_mysize(1) = ncoa*ncob*ncoc
238 CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
239
240 IF (PRESENT(int_abc_ext)) THEN
241 DO k = 1, ncoc
242 p1 = (k - 1)*ncoa
243 DO i = 1, ncoa
244 p2 = (p1 + i - 1)*ncob
245 DO j = 1, ncob
246 p3 = p2 + j
247 int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
248 int_abc_ext = max(int_abc_ext, abs(p_work(p3)))
249 END DO
250 END DO
251 END DO
252 ELSE
253 DO k = 1, ncoc
254 p1 = (k - 1)*ncoa
255 DO i = 1, ncoa
256 p2 = (p1 + i - 1)*ncob
257 DO j = 1, ncob
258 p3 = p2 + j
259 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 DEALLOCATE (params)
274
275 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 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 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 IF (PRESENT(params_in)) THEN
323 params => params_in
324
325 !If no external parameters to use, compute them
326 ELSE
327 params => params_out
328
329 !Note: some variable of 4-center integrals simplify with a dummy center:
330 ! P -> rk, gammap -> zetk
331 params%m_max = li_max + lj_max + lk_max
332 gammaq = zeti + zetj
333 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
334 params%ZetapEtaInv = 1._dp/(zetk + gammaq)
335
336 params%Q = (zeti*ri + zetj*rj)*params%EtaInv
337 params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
338 params%Rho = zetk*gammaq/(zetk + gammaq)
339
340 params%Fm = 0.0_dp
341 SELECT CASE (potential_parameter%potential_type)
342 CASE (do_potential_coulomb)
343 t = params%Rho*sum((params%Q - rk)**2)
344 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
345 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
346
347 CALL fgamma(params%m_max, t, params%Fm)
348 params%Fm = prefac*params%Fm
349 CASE (do_potential_truncated)
350 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
351 t = params%Rho*sum((params%Q - rk)**2)
352 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
353 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
354
355 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
356 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
357 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
358 params%Fm = prefac*params%Fm
359 CASE (do_potential_short)
360 t = params%Rho*sum((params%Q - rk)**2)
361 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
362 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
363
364 CALL fgamma(params%m_max, t, params%Fm)
365
366 omega2 = potential_parameter%omega**2
367 omega_corr2 = omega2/(omega2 + params%Rho)
368 omega_corr = sqrt(omega_corr2)
369 t = t*omega_corr2
370 ALLOCATE (fm(prim_data_f_size))
371
372 CALL fgamma(params%m_max, t, fm)
373 tmp = -omega_corr
374 DO l = 1, params%m_max + 1
375 params%Fm(l) = params%Fm(l) + fm(l)*tmp
376 tmp = tmp*omega_corr2
377 END DO
378 params%Fm = prefac*params%Fm
379 CASE (do_potential_mix_cl_trunc)
380 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
381 t = params%Rho*sum((params%Q - rk)**2)
382 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
383 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
384
385 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
386 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
387 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
388
389 ALLOCATE (fm(prim_data_f_size))
390 CALL fgamma(params%m_max, t, fm)
391 DO l = 1, params%m_max + 1
392 params%Fm(l) = params%Fm(l) &
393 *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
394 - fm(l)*potential_parameter%scale_longrange
395 END DO
396 DEALLOCATE (fm)
397
398 omega2 = potential_parameter%omega**2
399 omega_corr2 = omega2/(omega2 + params%Rho)
400 omega_corr = sqrt(omega_corr2)
401 t = t*omega_corr2
402
403 ALLOCATE (fm(prim_data_f_size))
404 CALL fgamma(params%m_max, t, fm)
405 tmp = omega_corr
406 DO l = 1, params%m_max + 1
407 params%Fm(l) = params%Fm(l) + fm(l)*tmp*potential_parameter%scale_longrange
408 tmp = tmp*omega_corr2
409 END DO
410 params%Fm = prefac*params%Fm
411 CASE (do_potential_id)
412 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2) &
413 - gammaq*zetk*params%ZetapEtaInv*sum((params%Q - rk)**2))
414 prefac = sqrt((pi*params%ZetapEtaInv)**3)*s1234
415
416 params%Fm(:) = prefac
417 CASE DEFAULT
418 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 params%m_max, params%Fm)
426
427 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 SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
465 la_min, la_max, npgfa, zeta, rpgfa, ra, &
466 lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
467 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 REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
493 TYPE(params_3c), POINTER :: params
494
495 NULLIFY (params, p_deriv)
496 ALLOCATE (params)
497
498 permute_1 = [4, 5, 6]
499 permute_2 = [7, 8, 9]
500
501 dr_ab = 0.0_dp
502 dr_bc = 0.0_dp
503 dr_ac = 0.0_dp
504
505 op = potential_parameter%potential_type
506
507 IF (op == do_potential_truncated .OR. op == do_potential_short &
508 .OR. op == do_potential_mix_cl_trunc) THEN
509 dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
510 dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
511 ELSE IF (op == do_potential_coulomb) THEN
512 dr_bc = 1000000.0_dp
513 dr_ac = 1000000.0_dp
514 END IF
515
516 do_ext = .false.
517 IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .true.
518 der_abc_1_ext_prv = 0.0_dp
519 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 DO ipgf = 1, npgfa
527 zeti = zeta(ipgf)
528 a_start = (ipgf - 1)*ncoset(la_max)
529
530 DO jpgf = 1, npgfb
531
532 ! screening
533 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
534
535 zetj = zetb(jpgf)
536 b_start = (jpgf - 1)*ncoset(lb_max)
537
538 DO kpgf = 1, npgfc
539
540 ! screening
541 IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) cycle
542 IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) cycle
543
544 zetk = zetc(kpgf)
545 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 potential_parameter=potential_parameter, params_out=params)
550
551 DO li = la_min, la_max
552 a_offset = a_start + ncoset(li - 1)
553 ncoa = nco(li)
554 DO lj = max(li, lb_min), lb_max
555 b_offset = b_start + ncoset(lj - 1)
556 ncob = nco(lj)
557 DO lk = lc_min, lc_max
558 c_offset = c_start + ncoset(lk - 1)
559 ncoc = nco(lk)
560
561 a_mysize(1) = ncoa*ncob*ncoc
562
563 CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
564
565 IF (do_ext) THEN
566 DO i_deriv = 1, 3
567 DO k = 1, ncoc
568 p1 = (k - 1)*ncob
569 DO j = 1, ncob
570 p2 = (p1 + j - 1)*ncoa
571 DO i = 1, ncoa
572 p3 = p2 + i
573
574 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
575 p_deriv(p3, permute_2(i_deriv))
576 der_abc_1_ext_prv(i_deriv) = max(der_abc_1_ext_prv(i_deriv), &
577 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 p_deriv(p3, permute_1(i_deriv))
581 der_abc_2_ext_prv(i_deriv) = max(der_abc_2_ext_prv(i_deriv), &
582 abs(p_deriv(p3, permute_1(i_deriv))))
583
584 END DO
585 END DO
586 END DO
587 END DO
588 ELSE
589 DO i_deriv = 1, 3
590 DO k = 1, ncoc
591 p1 = (k - 1)*ncob
592 DO j = 1, ncob
593 p2 = (p1 + j - 1)*ncoa
594 DO i = 1, ncoa
595 p3 = p2 + i
596
597 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
598 p_deriv(p3, permute_2(i_deriv))
599
600 der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
601 p_deriv(p3, permute_1(i_deriv))
602 END DO
603 END DO
604 END DO
605 END DO
606 END IF
607
608 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 CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
615
616 DO lj = lb_min, lb_max
617 b_offset = b_start + ncoset(lj - 1)
618 ncob = nco(lj)
619 DO li = max(lj + 1, la_min), la_max
620 a_offset = a_start + ncoset(li - 1)
621 ncoa = nco(li)
622 DO lk = lc_min, lc_max
623 c_offset = c_start + ncoset(lk - 1)
624 ncoc = nco(lk)
625
626 a_mysize(1) = ncoa*ncob*ncoc
627 CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
628
629 IF (do_ext) THEN
630 DO i_deriv = 1, 3
631 DO k = 1, ncoc
632 p1 = (k - 1)*ncoa
633 DO i = 1, ncoa
634 p2 = (p1 + i - 1)*ncob
635 DO j = 1, ncob
636 p3 = p2 + j
637
638 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
639 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 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 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 abs(p_deriv(p3, permute_2(i_deriv))))
649 END DO
650 END DO
651 END DO
652 END DO
653 ELSE
654 DO i_deriv = 1, 3
655 DO k = 1, ncoc
656 p1 = (k - 1)*ncoa
657 DO i = 1, ncoa
658 p2 = (p1 + i - 1)*ncob
659 DO j = 1, ncob
660 p3 = p2 + j
661
662 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
663 p_deriv(p3, permute_1(i_deriv))
664
665 der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
666 p_deriv(p3, permute_2(i_deriv))
667 END DO
668 END DO
669 END DO
670 END DO
671 END IF
672
673 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 IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
683 IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
684
685 DEALLOCATE (params)
686
687 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:) :: fm
723 TYPE(params_3c), POINTER :: params
724
725 IF (PRESENT(params_in)) THEN
726 params => params_in
727
728 ELSE
729 params => params_out
730
731 params%m_max = li_max + lj_max + lk_max + 1
732 gammaq = zeti + zetj
733 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
734 params%ZetapEtaInv = 1._dp/(zetk + gammaq)
735
736 params%Q = (zeti*ri + zetj*rj)*params%EtaInv
737 params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
738 params%Rho = zetk*gammaq/(zetk + gammaq)
739
740 params%Fm = 0.0_dp
741 SELECT CASE (potential_parameter%potential_type)
742 CASE (do_potential_coulomb)
743 t = params%Rho*sum((params%Q - rk)**2)
744 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
745 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
746
747 CALL fgamma(params%m_max, t, params%Fm)
748 params%Fm = prefac*params%Fm
749 CASE (do_potential_truncated)
750 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
751 t = params%Rho*sum((params%Q - rk)**2)
752 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
753 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
754
755 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
756 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
757 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
758 params%Fm = prefac*params%Fm
759 CASE (do_potential_short)
760 t = params%Rho*sum((params%Q - rk)**2)
761 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
762 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
763
764 CALL fgamma(params%m_max, t, params%Fm)
765
766 omega2 = potential_parameter%omega**2
767 omega_corr2 = omega2/(omega2 + params%Rho)
768 omega_corr = sqrt(omega_corr2)
769 t = t*omega_corr2
770 ALLOCATE (fm(prim_data_f_size))
771
772 CALL fgamma(params%m_max, t, fm)
773 tmp = -omega_corr
774 DO l = 1, params%m_max + 1
775 params%Fm(l) = params%Fm(l) + fm(l)*tmp
776 tmp = tmp*omega_corr2
777 END DO
778 params%Fm = prefac*params%Fm
779 CASE (do_potential_mix_cl_trunc)
780 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
781 t = params%Rho*sum((params%Q - rk)**2)
782 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
783 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
784
785 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
786 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
787 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
788
789 ALLOCATE (fm(prim_data_f_size))
790 CALL fgamma(params%m_max, t, fm)
791 DO l = 1, params%m_max + 1
792 params%Fm(l) = params%Fm(l) &
793 *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
794 - fm(l)*potential_parameter%scale_longrange
795 END DO
796 DEALLOCATE (fm)
797
798 omega2 = potential_parameter%omega**2
799 omega_corr2 = omega2/(omega2 + params%Rho)
800 omega_corr = sqrt(omega_corr2)
801 t = t*omega_corr2
802
803 ALLOCATE (fm(prim_data_f_size))
804 CALL fgamma(params%m_max, t, fm)
805 tmp = omega_corr
806 DO l = 1, params%m_max + 1
807 params%Fm(l) = params%Fm(l) + fm(l)*tmp*potential_parameter%scale_longrange
808 tmp = tmp*omega_corr2
809 END DO
810 params%Fm = prefac*params%Fm
811 CASE (do_potential_id)
812 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2) &
813 - gammaq*zetk*params%ZetapEtaInv*sum((params%Q - rk)**2))
814 prefac = sqrt((pi*params%ZetapEtaInv)**3)*s1234
815
816 params%Fm(:) = prefac
817 CASE DEFAULT
818 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 params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
826
827 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 SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
853 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 REAL(dp), DIMENSION(:), POINTER :: p_work
872
873 NULLIFY (p_work)
874
875 dr_ab = 0.0_dp
876
877 IF (potential_parameter%potential_type == do_potential_truncated .OR. &
878 potential_parameter%potential_type == do_potential_short .OR. &
879 potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
880 dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
881 ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
882 dr_ab = 1000000.0_dp
883 END IF
884
885 !Looping over the pgfs
886 DO ipgf = 1, npgfa
887 zeti = zeta(ipgf)
888 a_start = (ipgf - 1)*ncoset(la_max)
889
890 DO jpgf = 1, npgfb
891 zetj = zetb(jpgf)
892 b_start = (jpgf - 1)*ncoset(lb_max)
893
894 !screening
895 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
896
897 CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
898
899 DO li = la_min, la_max
900 a_offset = a_start + ncoset(li - 1)
901 ncoa = nco(li)
902 DO lj = lb_min, lb_max
903 b_offset = b_start + ncoset(lj - 1)
904 ncob = nco(lj)
905
906 a_mysize(1) = ncoa*ncob
907 CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
908
909 DO j = 1, ncob
910 p1 = (j - 1)*ncoa
911 DO i = 1, ncoa
912 p2 = p1 + i
913 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 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 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 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 op = potential_parameter%potential_type
961 params%m_max = lj_max + lk_max
962 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
963 params%ZetapEtaInv = 1._dp/(zetk + zetj)
964
965 params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
966 params%Rho = zetk*zetj/(zetk + zetj)
967
968 params%Fm = 0.0_dp
969 SELECT CASE (op)
970 CASE (do_potential_coulomb)
971 t = params%Rho*sum((rj - rk)**2)
972 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
973 CALL fgamma(params%m_max, t, params%Fm)
974 params%Fm = prefac*params%Fm
975 CASE (do_potential_truncated)
976 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
977 t = params%Rho*sum((rj - rk)**2)
978 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
979
980 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
981 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
982 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
983 params%Fm = prefac*params%Fm
984 CASE (do_potential_short)
985 t = params%Rho*sum((rj - rk)**2)
986 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
987
988 CALL fgamma(params%m_max, t, params%Fm)
989
990 omega2 = potential_parameter%omega**2
991 omega_corr2 = omega2/(omega2 + params%Rho)
992 omega_corr = sqrt(omega_corr2)
993 t = t*omega_corr2
994 ALLOCATE (fm(prim_data_f_size))
995
996 CALL fgamma(params%m_max, t, fm)
997 tmp = -omega_corr
998 DO l = 1, params%m_max + 1
999 params%Fm(l) = params%Fm(l) + fm(l)*tmp
1000 tmp = tmp*omega_corr2
1001 END DO
1002 params%Fm = prefac*params%Fm
1003 CASE (do_potential_mix_cl_trunc)
1004 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
1005 t = params%Rho*sum((rj - rk)**2)
1006 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1007
1008 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
1009 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
1010 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
1011
1012 ALLOCATE (fm(prim_data_f_size))
1013 CALL fgamma(params%m_max, t, fm)
1014 DO l = 1, params%m_max + 1
1015 params%Fm(l) = params%Fm(l) &
1016 *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
1017 - fm(l)*potential_parameter%scale_longrange
1018 END DO
1019 DEALLOCATE (fm)
1020
1021 omega2 = potential_parameter%omega**2
1022 omega_corr2 = omega2/(omega2 + params%Rho)
1023 omega_corr = sqrt(omega_corr2)
1024 t = t*omega_corr2
1025
1026 ALLOCATE (fm(prim_data_f_size))
1027 CALL fgamma(params%m_max, t, fm)
1028 tmp = omega_corr
1029 DO l = 1, params%m_max + 1
1030 params%Fm(l) = params%Fm(l) + fm(l)*tmp*potential_parameter%scale_longrange
1031 tmp = tmp*omega_corr2
1032 END DO
1033 params%Fm = prefac*params%Fm
1034 CASE (do_potential_id)
1035
1036 prefac = sqrt((pi*params%ZetapEtaInv)**3)*exp(-zetj*zetk*params%ZetapEtaInv*sum((rk - rj)**2))
1037 params%Fm(:) = prefac
1038 CASE DEFAULT
1039 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 params%m_max, params%Fm)
1045
1046 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 PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
1055 TYPE(coulomb_operator_type), INTENT(IN) :: potential1, potential2
1056 LOGICAL :: equals
1057
1058 IF (potential1%potential_type /= potential2%potential_type) THEN
1059 equals = .false.
1060 ELSE
1061 equals = .true.
1062 SELECT CASE (potential1%potential_type)
1063 CASE (do_potential_short, do_potential_long)
1064 IF (potential1%omega /= potential2%omega) equals = .false.
1065 CASE (do_potential_truncated)
1066 IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .false.
1067 CASE (do_potential_mix_cl_trunc)
1068 IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .false.
1069 IF (potential1%omega /= potential2%omega) equals = .false.
1070 IF (potential1%scale_coulomb /= potential2%scale_coulomb) equals = .false.
1071 IF (potential1%scale_longrange /= potential2%scale_longrange) equals = .false.
1072 END SELECT
1073 END IF
1074
1075 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 SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
1100 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 REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
1120
1121 NULLIFY (p_deriv)
1122
1123 permute = [4, 5, 6]
1124
1125 dr_ab = 0.0_dp
1126
1127 IF (potential_parameter%potential_type == do_potential_truncated .OR. &
1128 potential_parameter%potential_type == do_potential_short .OR. &
1129 potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
1130 dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
1131 ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
1132 dr_ab = 1000000.0_dp
1133 END IF
1134
1135 !Looping over the pgfs
1136 DO ipgf = 1, npgfa
1137 zeti = zeta(ipgf)
1138 a_start = (ipgf - 1)*ncoset(la_max)
1139
1140 DO jpgf = 1, npgfb
1141 zetj = zetb(jpgf)
1142 b_start = (jpgf - 1)*ncoset(lb_max)
1143
1144 !screening
1145 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
1146
1147 CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
1148
1149 DO li = la_min, la_max
1150 a_offset = a_start + ncoset(li - 1)
1151 ncoa = nco(li)
1152 DO lj = lb_min, lb_max
1153 b_offset = b_start + ncoset(lj - 1)
1154 ncob = nco(lj)
1155
1156 a_mysize(1) = ncoa*ncob
1157 CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
1158
1159 DO i_deriv = 1, 3
1160 DO j = 1, ncob
1161 p1 = (j - 1)*ncoa
1162 DO i = 1, ncoa
1163 p2 = p1 + i
1164 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 DEALLOCATE (p_deriv)
1170 END DO
1171 END DO
1172
1173 END DO
1174 END DO
1175
1176 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 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 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 op = potential_parameter%potential_type
1214 params%m_max = lj_max + lk_max + 1
1215 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
1216 params%ZetapEtaInv = 1._dp/(zetk + zetj)
1217
1218 params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
1219 params%Rho = zetk*zetj/(zetk + zetj)
1220
1221 params%Fm = 0.0_dp
1222 SELECT CASE (op)
1223 CASE (do_potential_coulomb)
1224 t = params%Rho*sum((rj - rk)**2)
1225 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1226 CALL fgamma(params%m_max, t, params%Fm)
1227 params%Fm = prefac*params%Fm
1228 CASE (do_potential_truncated)
1229 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
1230 t = params%Rho*sum((rj - rk)**2)
1231 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1232
1233 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
1234 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
1235 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
1236 params%Fm = prefac*params%Fm
1237 CASE (do_potential_short)
1238 t = params%Rho*sum((rj - rk)**2)
1239 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1240
1241 CALL fgamma(params%m_max, t, params%Fm)
1242
1243 omega2 = potential_parameter%omega**2
1244 omega_corr2 = omega2/(omega2 + params%Rho)
1245 omega_corr = sqrt(omega_corr2)
1246 t = t*omega_corr2
1247 ALLOCATE (fm(prim_data_f_size))
1248
1249 CALL fgamma(params%m_max, t, fm)
1250 tmp = -omega_corr
1251 DO l = 1, params%m_max + 1
1252 params%Fm(l) = params%Fm(l) + fm(l)*tmp
1253 tmp = tmp*omega_corr2
1254 END DO
1255 params%Fm = prefac*params%Fm
1256 CASE (do_potential_mix_cl_trunc)
1257 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
1258 t = params%Rho*sum((rj - rk)**2)
1259 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1260
1261 cpassert(get_lmax_init() >= params%m_max) !check if truncated coulomb init correctly
1262 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
1263 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
1264
1265 ALLOCATE (fm(prim_data_f_size))
1266 CALL fgamma(params%m_max, t, fm)
1267 DO l = 1, params%m_max + 1
1268 params%Fm(l) = params%Fm(l) &
1269 *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
1270 - fm(l)*potential_parameter%scale_longrange
1271 END DO
1272 DEALLOCATE (fm)
1273
1274 omega2 = potential_parameter%omega**2
1275 omega_corr2 = omega2/(omega2 + params%Rho)
1276 omega_corr = sqrt(omega_corr2)
1277 t = t*omega_corr2
1278
1279 ALLOCATE (fm(prim_data_f_size))
1280 CALL fgamma(params%m_max, t, fm)
1281 tmp = omega_corr
1282 DO l = 1, params%m_max + 1
1283 params%Fm(l) = params%Fm(l) + fm(l)*tmp*potential_parameter%scale_longrange
1284 tmp = tmp*omega_corr2
1285 END DO
1286 params%Fm = prefac*params%Fm
1287 CASE (do_potential_id)
1288
1289 prefac = sqrt((pi*params%ZetapEtaInv)**3)*exp(-zetj*zetk*params%ZetapEtaInv*sum((rk - rj)**2))
1290 params%Fm(:) = prefac
1291 CASE DEFAULT
1292 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 params%m_max, params%Fm)
1299
1300 END SUBROUTINE set_params_2c_deriv
1301
1302END MODULE libint_2c_3c
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
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:839
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:839
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:839
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:839
input_constants::do_potential_mix_cl_trunc
integer, parameter, public do_potential_mix_cl_trunc
Definition input_constants.F:839
input_constants::do_potential_long
integer, parameter, public do_potential_long
Definition input_constants.F:839
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::dp
integer, parameter, public dp
Definition kinds.F:34
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_2c_3c::eri_2center
subroutine, public eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, dab, lib, potential_parameter)
Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian gaussian orbita...
Definition libint_2c_3c.F:855
libint_2c_3c::compare_potential_types
pure logical function, public compare_potential_types(potential1, potential2)
Helper function to compare Coulomb operator types.
Definition libint_2c_3c.F:1055
libint_2c_3c::eri_3center
subroutine, public eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, int_abc_ext)
Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian gaussian orbit...
Definition libint_2c_3c.F:108
libint_2c_3c::cutoff_screen_factor
real(kind=dp), parameter, public cutoff_screen_factor
Definition libint_2c_3c.F:50
libint_2c_3c::eri_2center_derivs
subroutine, public eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, dab, lib, potential_parameter)
Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given set of cartes...
Definition libint_2c_3c.F:1102
libint_2c_3c::eri_3center_derivs
subroutine, public eri_3center_derivs(der_abc_1, der_abc_2, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, der_abc_1_ext, der_abc_2_ext)
Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given set of carte...
Definition libint_2c_3c.F:470
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
libint_wrapper::cp_libint_set_params_eri_deriv
subroutine, public cp_libint_set_params_eri_deriv(libint, a, b, c, d, p, q, w, zeta_a, zeta_b, zeta_c, zeta_d, zetainv, etainv, zetapetainv, rho, m_max, f)
Definition libint_wrapper.F:368
libint_wrapper::cp_libint_get_2eri_derivs
subroutine, public cp_libint_get_2eri_derivs(n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:1063
libint_wrapper::cp_libint_set_params_eri
subroutine, public cp_libint_set_params_eri(libint, a, b, c, d, zetainv, etainv, zetapetainv, rho, p, q, w, m_max, f)
Definition libint_wrapper.F:649
libint_wrapper::cp_libint_get_3eris
subroutine, public cp_libint_get_3eris(n_c, n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:977
libint_wrapper::cp_libint_get_2eris
subroutine, public cp_libint_get_2eris(n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:1105
libint_wrapper::prim_data_f_size
integer, parameter, public prim_data_f_size
Definition libint_wrapper.F:65
libint_wrapper::cp_libint_get_3eri_derivs
subroutine, public cp_libint_get_3eri_derivs(n_c, n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:1020
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
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::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:58
t_c_g0::t_c_g0_n
subroutine, public t_c_g0_n(res, use_gamma, r, t, nderiv)
...
Definition t_c_g0.F:88
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition t_c_g0.F:1468
integral_library_types::coulomb_operator_type
Definition integral_library_types.F:28
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:71