LCOV - code coverage report
Current view: top level - src - libint_2c_3c.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 89.4 % 536 479
Test Date: 2025-07-25 12:55:17 Functions: 75.0 % 12 9

            Line data    Source code
       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              : 
      14              : MODULE libint_2c_3c
      15              :    USE gamma,                           ONLY: fgamma => fgamma_0
      16              :    USE input_constants,                 ONLY: do_potential_coulomb,&
      17              :                                               do_potential_id,&
      18              :                                               do_potential_long,&
      19              :                                               do_potential_mix_cl_trunc,&
      20              :                                               do_potential_short,&
      21              :                                               do_potential_truncated
      22              :    USE 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              : 
      76              : CONTAINS
      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      2163603 :    SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
     111      2163603 :                           lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
     112      2163603 :                           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      2163603 :       REAL(dp), DIMENSION(:), POINTER                    :: p_work
     135              :       TYPE(params_3c), POINTER                           :: params
     136              : 
     137      2163603 :       NULLIFY (params, p_work)
     138     64908090 :       ALLOCATE (params)
     139              : 
     140      2163603 :       dr_ab = 0.0_dp
     141      2163603 :       dr_bc = 0.0_dp
     142      2163603 :       dr_ac = 0.0_dp
     143              : 
     144      2163603 :       op = potential_parameter%potential_type
     145              : 
     146              :       IF (op == do_potential_truncated .OR. op == do_potential_short &
     147      2163603 :           .OR. op == do_potential_mix_cl_trunc) THEN
     148      1347753 :          dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
     149      1347753 :          dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
     150       815850 :       ELSEIF (op == do_potential_coulomb) THEN
     151        19291 :          dr_bc = 1000000.0_dp
     152        19291 :          dr_ac = 1000000.0_dp
     153              :       END IF
     154              : 
     155      2163603 :       IF (PRESENT(int_abc_ext)) THEN
     156      2127951 :          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      6980806 :       DO ipgf = 1, npgfa
     165      4817203 :          zeti = zeta(ipgf)
     166      4817203 :          a_start = (ipgf - 1)*ncoset(la_max)
     167              : 
     168     20249252 :          DO jpgf = 1, npgfb
     169              : 
     170              :             ! screening
     171     13268446 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
     172              : 
     173      6292252 :             zetj = zetb(jpgf)
     174      6292252 :             b_start = (jpgf - 1)*ncoset(lb_max)
     175              : 
     176     30315371 :             DO kpgf = 1, npgfc
     177              : 
     178              :                ! screening
     179     19205916 :                IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
     180     11921910 :                IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
     181              : 
     182      9427201 :                zetk = zetc(kpgf)
     183      9427201 :                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      9427201 :                                   potential_parameter=potential_parameter, params_out=params)
     188              : 
     189     21517482 :                DO li = la_min, la_max
     190     12090281 :                   a_offset = a_start + ncoset(li - 1)
     191     12090281 :                   ncoa = nco(li)
     192     33828431 :                   DO lj = MAX(li, lb_min), lb_max
     193     12310949 :                      b_offset = b_start + ncoset(lj - 1)
     194     12310949 :                      ncob = nco(lj)
     195     46255097 :                      DO lk = lc_min, lc_max
     196     21853867 :                         c_offset = c_start + ncoset(lk - 1)
     197     21853867 :                         ncoc = nco(lk)
     198              : 
     199     21853867 :                         a_mysize(1) = ncoa*ncob*ncoc
     200     21853867 :                         CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
     201              : 
     202     34164816 :                         IF (PRESENT(int_abc_ext)) THEN
     203     86109397 :                            DO k = 1, ncoc
     204     64989989 :                               p1 = (k - 1)*ncob
     205    230224999 :                               DO j = 1, ncob
     206    144115602 :                                  p2 = (p1 + j - 1)*ncoa
     207    444475001 :                                  DO i = 1, ncoa
     208    235369410 :                                     p3 = p2 + i
     209    235369410 :                                     int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
     210    379485012 :                                     int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
     211              :                                  END DO
     212              :                               END DO
     213              :                            END DO
     214              :                         ELSE
     215      3219446 :                            DO k = 1, ncoc
     216      2484987 :                               p1 = (k - 1)*ncob
     217      8203101 :                               DO j = 1, ncob
     218      4983655 :                                  p2 = (p1 + j - 1)*ncoa
     219     14964443 :                                  DO i = 1, ncoa
     220      7495801 :                                     p3 = p2 + i
     221     12479456 :                                     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      9427201 :                CALL set_params_3c(lib, rb, ra, rc, params_in=params)
     233              : 
     234     34808502 :                DO lj = lb_min, lb_max
     235     12112855 :                   b_offset = b_start + ncoset(lj - 1)
     236     12112855 :                   ncob = nco(lj)
     237     34935133 :                   DO li = MAX(lj + 1, la_min), la_max
     238      3616362 :                      a_offset = a_start + ncoset(li - 1)
     239      3616362 :                      ncoa = nco(li)
     240     22633106 :                      DO lk = lc_min, lc_max
     241      6903889 :                         c_offset = c_start + ncoset(lk - 1)
     242      6903889 :                         ncoc = nco(lk)
     243              : 
     244      6903889 :                         a_mysize(1) = ncoa*ncob*ncoc
     245      6903889 :                         CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
     246              : 
     247     10520251 :                         IF (PRESENT(int_abc_ext)) THEN
     248     27954368 :                            DO k = 1, ncoc
     249     21261697 :                               p1 = (k - 1)*ncoa
     250    108278160 :                               DO i = 1, ncoa
     251     80323792 :                                  p2 = (p1 + i - 1)*ncob
     252    204374325 :                                  DO j = 1, ncob
     253    102788836 :                                     p3 = p2 + j
     254    102788836 :                                     int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
     255    183112628 :                                     int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
     256              :                                  END DO
     257              :                               END DO
     258              :                            END DO
     259              :                         ELSE
     260       933928 :                            DO k = 1, ncoc
     261       722710 :                               p1 = (k - 1)*ncoa
     262      3520291 :                               DO i = 1, ncoa
     263      2586363 :                                  p2 = (p1 + i - 1)*ncob
     264      6422524 :                                  DO j = 1, ncob
     265      3113451 :                                     p3 = p2 + j
     266      5699814 :                                     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      2163603 :       DEALLOCATE (params)
     281              : 
     282      2163603 :    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     18854402 :    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     18854402 :       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     18854402 :       IF (PRESENT(params_in)) THEN
     330      9427201 :          params => params_in
     331              : 
     332              :          !If no external parameters to use, compute them
     333              :       ELSE
     334      9427201 :          params => params_out
     335              : 
     336              :          !Note: some variable of 4-center integrals simplify with a dummy center:
     337              :          !      P -> rk, gammap -> zetk
     338      9427201 :          params%m_max = li_max + lj_max + lk_max
     339      9427201 :          gammaq = zeti + zetj
     340      9427201 :          params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
     341      9427201 :          params%ZetapEtaInv = 1._dp/(zetk + gammaq)
     342              : 
     343     37708804 :          params%Q = (zeti*ri + zetj*rj)*params%EtaInv
     344     37708804 :          params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
     345      9427201 :          params%Rho = zetk*gammaq/(zetk + gammaq)
     346              : 
     347    207398422 :          params%Fm = 0.0_dp
     348      9427201 :          SELECT CASE (potential_parameter%potential_type)
     349              :          CASE (do_potential_coulomb)
     350      1632720 :             T = params%Rho*SUM((params%Q - rk)**2)
     351      1632720 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     352       408180 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     353              : 
     354       408180 :             CALL fgamma(params%m_max, T, params%Fm)
     355      8979960 :             params%Fm = prefac*params%Fm
     356              :          CASE (do_potential_truncated)
     357      4002001 :             R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     358     16008004 :             T = params%Rho*SUM((params%Q - rk)**2)
     359     16008004 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     360      4002001 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     361              : 
     362      4002001 :             CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     363      4002001 :             CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     364      4002001 :             IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     365     88044022 :             params%Fm = prefac*params%Fm
     366              :          CASE (do_potential_short)
     367      8439944 :             T = params%Rho*SUM((params%Q - rk)**2)
     368      8439944 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     369      2109986 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     370              : 
     371      2109986 :             CALL fgamma(params%m_max, T, params%Fm)
     372              : 
     373      2109986 :             omega2 = potential_parameter%omega**2
     374      2109986 :             omega_corr2 = omega2/(omega2 + params%Rho)
     375      2109986 :             omega_corr = SQRT(omega_corr2)
     376      2109986 :             T = T*omega_corr2
     377      2109986 :             ALLOCATE (Fm(prim_data_f_size))
     378              : 
     379      2109986 :             CALL fgamma(params%m_max, T, Fm)
     380      2109986 :             tmp = -omega_corr
     381     11707924 :             DO l = 1, params%m_max + 1
     382      9597938 :                params%Fm(l) = params%Fm(l) + Fm(l)*tmp
     383     11707924 :                tmp = tmp*omega_corr2
     384              :             END DO
     385     46419692 :             params%Fm = prefac*params%Fm
     386              :          CASE (do_potential_mix_cl_trunc)
     387      1399165 :             R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     388      5596660 :             T = params%Rho*SUM((params%Q - rk)**2)
     389      5596660 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     390      1399165 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     391              : 
     392      1399165 :             CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     393      1399165 :             CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     394      1399165 :             IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     395              : 
     396      1399165 :             ALLOCATE (Fm(prim_data_f_size))
     397      1399165 :             CALL fgamma(params%m_max, T, Fm)
     398      5377685 :             DO l = 1, params%m_max + 1
     399              :                params%Fm(l) = params%Fm(l) &
     400              :                               *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
     401      5377685 :                               - Fm(l)*potential_parameter%scale_longrange
     402              :             END DO
     403      1399165 :             DEALLOCATE (Fm)
     404              : 
     405      1399165 :             omega2 = potential_parameter%omega**2
     406      1399165 :             omega_corr2 = omega2/(omega2 + params%Rho)
     407      1399165 :             omega_corr = SQRT(omega_corr2)
     408      1399165 :             T = T*omega_corr2
     409              : 
     410      1399165 :             ALLOCATE (Fm(prim_data_f_size))
     411      1399165 :             CALL fgamma(params%m_max, T, Fm)
     412      1399165 :             tmp = omega_corr
     413      5377685 :             DO l = 1, params%m_max + 1
     414      3978520 :                params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
     415      5377685 :                tmp = tmp*omega_corr2
     416              :             END DO
     417     30781630 :             params%Fm = prefac*params%Fm
     418              :          CASE (do_potential_id)
     419              :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
     420     10555083 :                         - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
     421      1507869 :             prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
     422              : 
     423     33173118 :             params%Fm(:) = prefac
     424              :          CASE DEFAULT
     425      9427201 :             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     18854402 :                                     params%m_max, params%Fm)
     433              : 
     434     18854402 :    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       342310 :    SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
     472       342310 :                                  la_min, la_max, npgfa, zeta, rpgfa, ra, &
     473       342310 :                                  lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
     474       342310 :                                  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       342310 :       REAL(dp), DIMENSION(:, :), POINTER                 :: p_deriv
     500              :       TYPE(params_3c), POINTER                           :: params
     501              : 
     502       342310 :       NULLIFY (params, p_deriv)
     503     10269300 :       ALLOCATE (params)
     504              : 
     505       342310 :       permute_1 = [4, 5, 6]
     506       342310 :       permute_2 = [7, 8, 9]
     507              : 
     508       342310 :       dr_ab = 0.0_dp
     509       342310 :       dr_bc = 0.0_dp
     510       342310 :       dr_ac = 0.0_dp
     511              : 
     512       342310 :       op = potential_parameter%potential_type
     513              : 
     514              :       IF (op == do_potential_truncated .OR. op == do_potential_short &
     515       342310 :           .OR. op == do_potential_mix_cl_trunc) THEN
     516       119807 :          dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
     517       119807 :          dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
     518       222503 :       ELSEIF (op == do_potential_coulomb) THEN
     519         9988 :          dr_bc = 1000000.0_dp
     520         9988 :          dr_ac = 1000000.0_dp
     521              :       END IF
     522              : 
     523       342310 :       do_ext = .FALSE.
     524       342310 :       IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
     525       342310 :       der_abc_1_ext_prv = 0.0_dp
     526       342310 :       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      1180747 :       DO ipgf = 1, npgfa
     534       838437 :          zeti = zeta(ipgf)
     535       838437 :          a_start = (ipgf - 1)*ncoset(la_max)
     536              : 
     537      4746758 :          DO jpgf = 1, npgfb
     538              : 
     539              :             ! screening
     540      3566011 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
     541              : 
     542      1129693 :             zetj = zetb(jpgf)
     543      1129693 :             b_start = (jpgf - 1)*ncoset(lb_max)
     544              : 
     545      8285048 :             DO kpgf = 1, npgfc
     546              : 
     547              :                ! screening
     548      6316918 :                IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
     549      2649783 :                IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
     550              : 
     551      1857820 :                zetk = zetc(kpgf)
     552      1857820 :                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      1857820 :                                         potential_parameter=potential_parameter, params_out=params)
     557              : 
     558      4281757 :                DO li = la_min, la_max
     559      2423937 :                   a_offset = a_start + ncoset(li - 1)
     560      2423937 :                   ncoa = nco(li)
     561      6872693 :                   DO lj = MAX(li, lb_min), lb_max
     562      2590936 :                      b_offset = b_start + ncoset(lj - 1)
     563      2590936 :                      ncob = nco(lj)
     564      8765957 :                      DO lk = lc_min, lc_max
     565      3751084 :                         c_offset = c_start + ncoset(lk - 1)
     566      3751084 :                         ncoc = nco(lk)
     567              : 
     568      3751084 :                         a_mysize(1) = ncoa*ncob*ncoc
     569              : 
     570      3751084 :                         CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
     571              : 
     572      3751084 :                         IF (do_ext) THEN
     573     15004336 :                            DO i_deriv = 1, 3
     574     42533851 :                               DO k = 1, ncoc
     575     27529515 :                                  p1 = (k - 1)*ncob
     576     90660276 :                                  DO j = 1, ncob
     577     51877509 :                                     p2 = (p1 + j - 1)*ncoa
     578    157726497 :                                     DO i = 1, ncoa
     579     78319473 :                                        p3 = p2 + i
     580              : 
     581              :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     582     78319473 :                                           p_deriv(p3, permute_2(i_deriv))
     583              :                                        der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
     584     78319473 :                                                                         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     78319473 :                                           p_deriv(p3, permute_1(i_deriv))
     588              :                                        der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
     589    130196982 :                                                                         ABS(p_deriv(p3, permute_1(i_deriv))))
     590              : 
     591              :                                     END DO
     592              :                                  END DO
     593              :                               END DO
     594              :                            END DO
     595              :                         ELSE
     596            0 :                            DO i_deriv = 1, 3
     597            0 :                               DO k = 1, ncoc
     598            0 :                                  p1 = (k - 1)*ncob
     599            0 :                                  DO j = 1, ncob
     600            0 :                                     p2 = (p1 + j - 1)*ncoa
     601            0 :                                     DO i = 1, ncoa
     602            0 :                                        p3 = p2 + i
     603              : 
     604              :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     605            0 :                                           p_deriv(p3, permute_2(i_deriv))
     606              : 
     607              :                                        der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     608            0 :                                           p_deriv(p3, permute_1(i_deriv))
     609              :                                     END DO
     610              :                                  END DO
     611              :                               END DO
     612              :                            END DO
     613              :                         END IF
     614              : 
     615      6342020 :                         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      1857820 :                CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
     622              : 
     623      7856965 :                DO lj = lb_min, lb_max
     624      2433134 :                   b_offset = b_start + ncoset(lj - 1)
     625      2433134 :                   ncob = nco(lj)
     626      9411651 :                   DO li = MAX(lj + 1, la_min), la_max
     627       661599 :                      a_offset = a_start + ncoset(li - 1)
     628       661599 :                      ncoa = nco(li)
     629      4065070 :                      DO lk = lc_min, lc_max
     630       970337 :                         c_offset = c_start + ncoset(lk - 1)
     631       970337 :                         ncoc = nco(lk)
     632              : 
     633       970337 :                         a_mysize(1) = ncoa*ncob*ncoc
     634       970337 :                         CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
     635              : 
     636       970337 :                         IF (do_ext) THEN
     637      3881348 :                            DO i_deriv = 1, 3
     638     11231387 :                               DO k = 1, ncoc
     639      7350039 :                                  p1 = (k - 1)*ncoa
     640     34394001 :                                  DO i = 1, ncoa
     641     24132951 :                                     p2 = (p1 + i - 1)*ncob
     642     59103765 :                                     DO j = 1, ncob
     643     27620775 :                                        p3 = p2 + j
     644              : 
     645              :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     646     27620775 :                                           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     27620775 :                                                                         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     27620775 :                                           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     51753726 :                                                                         ABS(p_deriv(p3, permute_2(i_deriv))))
     656              :                                     END DO
     657              :                                  END DO
     658              :                               END DO
     659              :                            END DO
     660              :                         ELSE
     661            0 :                            DO i_deriv = 1, 3
     662            0 :                               DO k = 1, ncoc
     663            0 :                                  p1 = (k - 1)*ncoa
     664            0 :                                  DO i = 1, ncoa
     665            0 :                                     p2 = (p1 + i - 1)*ncob
     666            0 :                                     DO j = 1, ncob
     667            0 :                                        p3 = p2 + j
     668              : 
     669              :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     670            0 :                                           p_deriv(p3, permute_1(i_deriv))
     671              : 
     672              :                                        der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     673            0 :                                           p_deriv(p3, permute_2(i_deriv))
     674              :                                     END DO
     675              :                                  END DO
     676              :                               END DO
     677              :                            END DO
     678              :                         END IF
     679              : 
     680      1631936 :                         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       342310 :       IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
     690       342310 :       IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
     691              : 
     692       342310 :       DEALLOCATE (params)
     693              : 
     694       342310 :    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      3715640 :    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      3715640 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: Fm
     730              :       TYPE(params_3c), POINTER                           :: params
     731              : 
     732      3715640 :       IF (PRESENT(params_in)) THEN
     733      1857820 :          params => params_in
     734              : 
     735              :       ELSE
     736      1857820 :          params => params_out
     737              : 
     738      1857820 :          params%m_max = li_max + lj_max + lk_max + 1
     739      1857820 :          gammaq = zeti + zetj
     740      1857820 :          params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
     741      1857820 :          params%ZetapEtaInv = 1._dp/(zetk + gammaq)
     742              : 
     743      7431280 :          params%Q = (zeti*ri + zetj*rj)*params%EtaInv
     744      7431280 :          params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
     745      1857820 :          params%Rho = zetk*gammaq/(zetk + gammaq)
     746              : 
     747     40872040 :          params%Fm = 0.0_dp
     748      1857820 :          SELECT CASE (potential_parameter%potential_type)
     749              :          CASE (do_potential_coulomb)
     750       416424 :             T = params%Rho*SUM((params%Q - rk)**2)
     751       416424 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     752       104106 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     753              : 
     754       104106 :             CALL fgamma(params%m_max, T, params%Fm)
     755      2290332 :             params%Fm = prefac*params%Fm
     756              :          CASE (do_potential_truncated)
     757      1244984 :             R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     758      4979936 :             T = params%Rho*SUM((params%Q - rk)**2)
     759      4979936 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     760      1244984 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     761              : 
     762      1244984 :             CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     763      1244984 :             CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     764      1244984 :             IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     765     27389648 :             params%Fm = prefac*params%Fm
     766              :          CASE (do_potential_short)
     767            0 :             T = params%Rho*SUM((params%Q - rk)**2)
     768            0 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     769            0 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     770              : 
     771            0 :             CALL fgamma(params%m_max, T, params%Fm)
     772              : 
     773            0 :             omega2 = potential_parameter%omega**2
     774            0 :             omega_corr2 = omega2/(omega2 + params%Rho)
     775            0 :             omega_corr = SQRT(omega_corr2)
     776            0 :             T = T*omega_corr2
     777            0 :             ALLOCATE (Fm(prim_data_f_size))
     778              : 
     779            0 :             CALL fgamma(params%m_max, T, Fm)
     780            0 :             tmp = -omega_corr
     781            0 :             DO l = 1, params%m_max + 1
     782            0 :                params%Fm(l) = params%Fm(l) + Fm(l)*tmp
     783            0 :                tmp = tmp*omega_corr2
     784              :             END DO
     785            0 :             params%Fm = prefac*params%Fm
     786              :          CASE (do_potential_mix_cl_trunc)
     787            0 :             R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     788            0 :             T = params%Rho*SUM((params%Q - rk)**2)
     789            0 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     790            0 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     791              : 
     792            0 :             CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     793            0 :             CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     794            0 :             IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     795              : 
     796            0 :             ALLOCATE (Fm(prim_data_f_size))
     797            0 :             CALL fgamma(params%m_max, T, Fm)
     798            0 :             DO l = 1, params%m_max + 1
     799              :                params%Fm(l) = params%Fm(l) &
     800              :                               *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
     801            0 :                               - Fm(l)*potential_parameter%scale_longrange
     802              :             END DO
     803            0 :             DEALLOCATE (Fm)
     804              : 
     805            0 :             omega2 = potential_parameter%omega**2
     806            0 :             omega_corr2 = omega2/(omega2 + params%Rho)
     807            0 :             omega_corr = SQRT(omega_corr2)
     808            0 :             T = T*omega_corr2
     809              : 
     810            0 :             ALLOCATE (Fm(prim_data_f_size))
     811            0 :             CALL fgamma(params%m_max, T, Fm)
     812            0 :             tmp = omega_corr
     813            0 :             DO l = 1, params%m_max + 1
     814            0 :                params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
     815            0 :                tmp = tmp*omega_corr2
     816              :             END DO
     817            0 :             params%Fm = prefac*params%Fm
     818              :          CASE (do_potential_id)
     819              :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
     820      3561110 :                         - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
     821       508730 :             prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
     822              : 
     823     11192060 :             params%Fm(:) = prefac
     824              :          CASE DEFAULT
     825      1857820 :             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      3715640 :                                           params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
     833              : 
     834      3715640 :    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       398419 :    SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
     860       398419 :                           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       398419 :       REAL(dp), DIMENSION(:), POINTER                    :: p_work
     879              : 
     880       398419 :       NULLIFY (p_work)
     881              : 
     882       398419 :       dr_ab = 0.0_dp
     883              : 
     884              :       IF (potential_parameter%potential_type == do_potential_truncated .OR. &
     885       180884 :           potential_parameter%potential_type == do_potential_short .OR. &
     886              :           potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
     887       230171 :          dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
     888       168248 :       ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
     889         1616 :          dr_ab = 1000000.0_dp
     890              :       END IF
     891              : 
     892              :       !Looping over the pgfs
     893      1292204 :       DO ipgf = 1, npgfa
     894       893785 :          zeti = zeta(ipgf)
     895       893785 :          a_start = (ipgf - 1)*ncoset(la_max)
     896              : 
     897      5721224 :          DO jpgf = 1, npgfb
     898      4429020 :             zetj = zetb(jpgf)
     899      4429020 :             b_start = (jpgf - 1)*ncoset(lb_max)
     900              : 
     901              :             !screening
     902      4429020 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
     903              : 
     904      2856168 :             CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
     905              : 
     906      9051913 :             DO li = la_min, la_max
     907      5301960 :                a_offset = a_start + ncoset(li - 1)
     908      5301960 :                ncoa = nco(li)
     909     20231143 :                DO lj = lb_min, lb_max
     910     10500163 :                   b_offset = b_start + ncoset(lj - 1)
     911     10500163 :                   ncob = nco(lj)
     912              : 
     913     10500163 :                   a_mysize(1) = ncoa*ncob
     914     10500163 :                   CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
     915              : 
     916     47859327 :                   DO j = 1, ncob
     917     32057204 :                      p1 = (j - 1)*ncoa
     918    142254741 :                      DO i = 1, ncoa
     919     99697374 :                         p2 = p1 + i
     920    131754578 :                         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       398419 :    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      2856168 :    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      2856168 :       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      2856168 :       op = potential_parameter%potential_type
     968      2856168 :       params%m_max = lj_max + lk_max
     969      2856168 :       params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
     970      2856168 :       params%ZetapEtaInv = 1._dp/(zetk + zetj)
     971              : 
     972     11424672 :       params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
     973      2856168 :       params%Rho = zetk*zetj/(zetk + zetj)
     974              : 
     975     62835696 :       params%Fm = 0.0_dp
     976              :       SELECT CASE (op)
     977              :       CASE (do_potential_coulomb)
     978       106992 :          T = params%Rho*SUM((rj - rk)**2)
     979        26748 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
     980        26748 :          CALL fgamma(params%m_max, T, params%Fm)
     981       588456 :          params%Fm = prefac*params%Fm
     982              :       CASE (do_potential_truncated)
     983       168755 :          R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     984       675020 :          T = params%Rho*SUM((rj - rk)**2)
     985       168755 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
     986              : 
     987       168755 :          CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     988       168755 :          CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     989       168755 :          IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     990      3712610 :          params%Fm = prefac*params%Fm
     991              :       CASE (do_potential_short)
     992     10017192 :          T = params%Rho*SUM((rj - rk)**2)
     993      2504298 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
     994              : 
     995      2504298 :          CALL fgamma(params%m_max, T, params%Fm)
     996              : 
     997      2504298 :          omega2 = potential_parameter%omega**2
     998      2504298 :          omega_corr2 = omega2/(omega2 + params%Rho)
     999      2504298 :          omega_corr = SQRT(omega_corr2)
    1000      2504298 :          T = T*omega_corr2
    1001      2504298 :          ALLOCATE (Fm(prim_data_f_size))
    1002              : 
    1003      2504298 :          CALL fgamma(params%m_max, T, Fm)
    1004      2504298 :          tmp = -omega_corr
    1005     12305810 :          DO l = 1, params%m_max + 1
    1006      9801512 :             params%Fm(l) = params%Fm(l) + Fm(l)*tmp
    1007     12305810 :             tmp = tmp*omega_corr2
    1008              :          END DO
    1009     55094556 :          params%Fm = prefac*params%Fm
    1010              :       CASE (do_potential_mix_cl_trunc)
    1011        61464 :          R = potential_parameter%cutoff_radius*SQRT(params%Rho)
    1012       245856 :          T = params%Rho*SUM((rj - rk)**2)
    1013        61464 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1014              : 
    1015        61464 :          CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
    1016        61464 :          CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
    1017        61464 :          IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
    1018              : 
    1019        61464 :          ALLOCATE (Fm(prim_data_f_size))
    1020        61464 :          CALL fgamma(params%m_max, T, Fm)
    1021       223681 :          DO l = 1, params%m_max + 1
    1022              :             params%Fm(l) = params%Fm(l) &
    1023              :                            *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
    1024       223681 :                            - Fm(l)*potential_parameter%scale_longrange
    1025              :          END DO
    1026        61464 :          DEALLOCATE (Fm)
    1027              : 
    1028        61464 :          omega2 = potential_parameter%omega**2
    1029        61464 :          omega_corr2 = omega2/(omega2 + params%Rho)
    1030        61464 :          omega_corr = SQRT(omega_corr2)
    1031        61464 :          T = T*omega_corr2
    1032              : 
    1033        61464 :          ALLOCATE (Fm(prim_data_f_size))
    1034        61464 :          CALL fgamma(params%m_max, T, Fm)
    1035        61464 :          tmp = omega_corr
    1036       223681 :          DO l = 1, params%m_max + 1
    1037       162217 :             params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
    1038       223681 :             tmp = tmp*omega_corr2
    1039              :          END DO
    1040      1352208 :          params%Fm = prefac*params%Fm
    1041              :       CASE (do_potential_id)
    1042              : 
    1043       379612 :          prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
    1044      2087866 :          params%Fm(:) = prefac
    1045              :       CASE DEFAULT
    1046      2856168 :          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      2856168 :                                     params%m_max, params%Fm)
    1052              : 
    1053     74260368 :    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        14370 :    PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
    1062              :       TYPE(libint_potential_type), INTENT(IN)            :: potential1, potential2
    1063              :       LOGICAL                                            :: equals
    1064              : 
    1065        14370 :       IF (potential1%potential_type /= potential2%potential_type) THEN
    1066              :          equals = .FALSE.
    1067              :       ELSE
    1068        13552 :          equals = .TRUE.
    1069          736 :          SELECT CASE (potential1%potential_type)
    1070              :          CASE (do_potential_short, do_potential_long)
    1071          736 :             IF (potential1%omega /= potential2%omega) equals = .FALSE.
    1072              :          CASE (do_potential_truncated)
    1073           14 :             IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
    1074              :          CASE (do_potential_mix_cl_trunc)
    1075            2 :             IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
    1076            2 :             IF (potential1%omega /= potential2%omega) equals = .FALSE.
    1077            2 :             IF (potential1%scale_coulomb /= potential2%scale_coulomb) equals = .FALSE.
    1078        13554 :             IF (potential1%scale_longrange /= potential2%scale_longrange) equals = .FALSE.
    1079              :          END SELECT
    1080              :       END IF
    1081              : 
    1082        14370 :    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       143159 :    SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
    1107       143159 :                                  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       143159 :       REAL(dp), DIMENSION(:, :), POINTER                 :: p_deriv
    1127              : 
    1128       143159 :       NULLIFY (p_deriv)
    1129              : 
    1130       143159 :       permute = [4, 5, 6]
    1131              : 
    1132       143159 :       dr_ab = 0.0_dp
    1133              : 
    1134              :       IF (potential_parameter%potential_type == do_potential_truncated .OR. &
    1135        79391 :           potential_parameter%potential_type == do_potential_short .OR. &
    1136              :           potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
    1137        64460 :          dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
    1138        78699 :       ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
    1139         9831 :          dr_ab = 1000000.0_dp
    1140              :       END IF
    1141              : 
    1142              :       !Looping over the pgfs
    1143       635114 :       DO ipgf = 1, npgfa
    1144       491955 :          zeti = zeta(ipgf)
    1145       491955 :          a_start = (ipgf - 1)*ncoset(la_max)
    1146              : 
    1147      3845505 :          DO jpgf = 1, npgfb
    1148      3210391 :             zetj = zetb(jpgf)
    1149      3210391 :             b_start = (jpgf - 1)*ncoset(lb_max)
    1150              : 
    1151              :             !screening
    1152      3210391 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
    1153              : 
    1154      2171616 :             CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
    1155              : 
    1156      6307276 :             DO li = la_min, la_max
    1157      3643705 :                a_offset = a_start + ncoset(li - 1)
    1158      3643705 :                ncoa = nco(li)
    1159     13017159 :                DO lj = lb_min, lb_max
    1160      6163063 :                   b_offset = b_start + ncoset(lj - 1)
    1161      6163063 :                   ncob = nco(lj)
    1162              : 
    1163      6163063 :                   a_mysize(1) = ncoa*ncob
    1164      6163063 :                   CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
    1165              : 
    1166     24652252 :                   DO i_deriv = 1, 3
    1167     75297271 :                      DO j = 1, ncob
    1168     50645019 :                         p1 = (j - 1)*ncoa
    1169    208323612 :                         DO i = 1, ncoa
    1170    139189404 :                            p2 = p1 + i
    1171    189834423 :                            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      9806768 :                   DEALLOCATE (p_deriv)
    1177              :                END DO
    1178              :             END DO
    1179              : 
    1180              :          END DO
    1181              :       END DO
    1182              : 
    1183       143159 :    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      2171616 :    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      2171616 :       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      2171616 :       op = potential_parameter%potential_type
    1221      2171616 :       params%m_max = lj_max + lk_max + 1
    1222      2171616 :       params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
    1223      2171616 :       params%ZetapEtaInv = 1._dp/(zetk + zetj)
    1224              : 
    1225      8686464 :       params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
    1226      2171616 :       params%Rho = zetk*zetj/(zetk + zetj)
    1227              : 
    1228     47775552 :       params%Fm = 0.0_dp
    1229              :       SELECT CASE (op)
    1230              :       CASE (do_potential_coulomb)
    1231       102916 :          T = params%Rho*SUM((rj - rk)**2)
    1232        25729 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1233        25729 :          CALL fgamma(params%m_max, T, params%Fm)
    1234       566038 :          params%Fm = prefac*params%Fm
    1235              :       CASE (do_potential_truncated)
    1236        61941 :          R = potential_parameter%cutoff_radius*SQRT(params%Rho)
    1237       247764 :          T = params%Rho*SUM((rj - rk)**2)
    1238        61941 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1239              : 
    1240        61941 :          CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
    1241        61941 :          CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
    1242        61941 :          IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
    1243      1362702 :          params%Fm = prefac*params%Fm
    1244              :       CASE (do_potential_short)
    1245      8096600 :          T = params%Rho*SUM((rj - rk)**2)
    1246      2024150 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1247              : 
    1248      2024150 :          CALL fgamma(params%m_max, T, params%Fm)
    1249              : 
    1250      2024150 :          omega2 = potential_parameter%omega**2
    1251      2024150 :          omega_corr2 = omega2/(omega2 + params%Rho)
    1252      2024150 :          omega_corr = SQRT(omega_corr2)
    1253      2024150 :          T = T*omega_corr2
    1254      2024150 :          ALLOCATE (Fm(prim_data_f_size))
    1255              : 
    1256      2024150 :          CALL fgamma(params%m_max, T, Fm)
    1257      2024150 :          tmp = -omega_corr
    1258     11369176 :          DO l = 1, params%m_max + 1
    1259      9345026 :             params%Fm(l) = params%Fm(l) + Fm(l)*tmp
    1260     11369176 :             tmp = tmp*omega_corr2
    1261              :          END DO
    1262     44531300 :          params%Fm = prefac*params%Fm
    1263              :       CASE (do_potential_mix_cl_trunc)
    1264         2578 :          R = potential_parameter%cutoff_radius*SQRT(params%Rho)
    1265        10312 :          T = params%Rho*SUM((rj - rk)**2)
    1266         2578 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1267              : 
    1268         2578 :          CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
    1269         2578 :          CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
    1270         2578 :          IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
    1271              : 
    1272         2578 :          ALLOCATE (Fm(prim_data_f_size))
    1273         2578 :          CALL fgamma(params%m_max, T, Fm)
    1274         8234 :          DO l = 1, params%m_max + 1
    1275              :             params%Fm(l) = params%Fm(l) &
    1276              :                            *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
    1277         8234 :                            - Fm(l)*potential_parameter%scale_longrange
    1278              :          END DO
    1279         2578 :          DEALLOCATE (Fm)
    1280              : 
    1281         2578 :          omega2 = potential_parameter%omega**2
    1282         2578 :          omega_corr2 = omega2/(omega2 + params%Rho)
    1283         2578 :          omega_corr = SQRT(omega_corr2)
    1284         2578 :          T = T*omega_corr2
    1285              : 
    1286         2578 :          ALLOCATE (Fm(prim_data_f_size))
    1287         2578 :          CALL fgamma(params%m_max, T, Fm)
    1288         2578 :          tmp = omega_corr
    1289         8234 :          DO l = 1, params%m_max + 1
    1290         5656 :             params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
    1291         8234 :             tmp = tmp*omega_corr2
    1292              :          END DO
    1293        56716 :          params%Fm = prefac*params%Fm
    1294              :       CASE (do_potential_id)
    1295              : 
    1296       228872 :          prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
    1297      1258796 :          params%Fm(:) = prefac
    1298              :       CASE DEFAULT
    1299      2171616 :          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      2171616 :                                           params%m_max, params%Fm)
    1306              : 
    1307     56462016 :    END SUBROUTINE set_params_2c_deriv
    1308              : 
    1309            0 : END MODULE libint_2c_3c
    1310              : 
        

Generated by: LCOV version 2.0-1