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