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