LCOV - code coverage report
Current view: top level - src - libint_2c_3c.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 408 442 92.3 %
Date: 2024-03-29 07:50:05 Functions: 9 12 75.0 %

          Line data    Source code
       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     2689329 :    SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
     110     2689329 :                           lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
     111     2689329 :                           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     2689329 :       REAL(dp), DIMENSION(:), POINTER                    :: p_work
     134             :       TYPE(params_3c), POINTER                           :: params
     135             : 
     136     2689329 :       NULLIFY (params, p_work)
     137    80679870 :       ALLOCATE (params)
     138             : 
     139     2689329 :       dr_ab = 0.0_dp
     140     2689329 :       dr_bc = 0.0_dp
     141     2689329 :       dr_ac = 0.0_dp
     142             : 
     143     2689329 :       op = potential_parameter%potential_type
     144             : 
     145     2689329 :       IF (op == do_potential_truncated .OR. op == do_potential_short) THEN
     146     1891992 :          dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
     147     1891992 :          dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
     148      797337 :       ELSEIF (op == do_potential_coulomb) THEN
     149       16371 :          dr_bc = 1000000.0_dp
     150       16371 :          dr_ac = 1000000.0_dp
     151             :       END IF
     152             : 
     153     2689329 :       IF (PRESENT(int_abc_ext)) THEN
     154     2658957 :          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     9137657 :       DO ipgf = 1, npgfa
     163     6448328 :          zeti = zeta(ipgf)
     164     6448328 :          a_start = (ipgf - 1)*ncoset(la_max)
     165             : 
     166    27405254 :          DO jpgf = 1, npgfb
     167             : 
     168             :             ! screening
     169    18267597 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
     170             : 
     171     9087856 :             zetj = zetb(jpgf)
     172     9087856 :             b_start = (jpgf - 1)*ncoset(lb_max)
     173             : 
     174    34044342 :             DO kpgf = 1, npgfc
     175             : 
     176             :                ! screening
     177    18508158 :                IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
     178    12371187 :                IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
     179             : 
     180    10223144 :                zetk = zetc(kpgf)
     181    10223144 :                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    10223144 :                                   potential_parameter=potential_parameter, params_out=params)
     186             : 
     187    24598294 :                DO li = la_min, la_max
     188    14375150 :                   a_offset = a_start + ncoset(li - 1)
     189    14375150 :                   ncoa = nco(li)
     190    40306514 :                   DO lj = MAX(li, lb_min), lb_max
     191    15708220 :                      b_offset = b_start + ncoset(lj - 1)
     192    15708220 :                      ncob = nco(lj)
     193    54097306 :                      DO lk = lc_min, lc_max
     194    24013936 :                         c_offset = c_start + ncoset(lk - 1)
     195    24013936 :                         ncoc = nco(lk)
     196             : 
     197    24013936 :                         a_mysize(1) = ncoa*ncob*ncoc
     198    24013936 :                         CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
     199             : 
     200    39722156 :                         IF (PRESENT(int_abc_ext)) THEN
     201    88249854 :                            DO k = 1, ncoc
     202    64868976 :                               p1 = (k - 1)*ncob
     203   234958580 :                               DO j = 1, ncob
     204   146708726 :                                  p2 = (p1 + j - 1)*ncoa
     205   455341648 :                                  DO i = 1, ncoa
     206   243763946 :                                     p3 = p2 + i
     207   243763946 :                                     int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
     208   390472672 :                                     int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
     209             :                                  END DO
     210             :                               END DO
     211             :                            END DO
     212             :                         ELSE
     213     2750604 :                            DO k = 1, ncoc
     214     2117546 :                               p1 = (k - 1)*ncob
     215     6941472 :                               DO j = 1, ncob
     216     4190868 :                                  p2 = (p1 + j - 1)*ncoa
     217    12548432 :                                  DO i = 1, ncoa
     218     6240018 :                                     p3 = p2 + i
     219    10430886 :                                     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    10223144 :                CALL set_params_3c(lib, rb, ra, rc, params_in=params)
     231             : 
     232    42890675 :                DO lj = lb_min, lb_max
     233    14399934 :                   b_offset = b_start + ncoset(lj - 1)
     234    14399934 :                   ncob = nco(lj)
     235    37925370 :                   DO li = MAX(lj + 1, la_min), la_max
     236     5017278 :                      a_offset = a_start + ncoset(li - 1)
     237     5017278 :                      ncoa = nco(li)
     238    27337213 :                      DO lk = lc_min, lc_max
     239     7920001 :                         c_offset = c_start + ncoset(lk - 1)
     240     7920001 :                         ncoc = nco(lk)
     241             : 
     242     7920001 :                         a_mysize(1) = ncoa*ncob*ncoc
     243     7920001 :                         CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
     244             : 
     245    12937279 :                         IF (PRESENT(int_abc_ext)) THEN
     246    29671975 :                            DO k = 1, ncoc
     247    21931700 :                               p1 = (k - 1)*ncoa
     248   110694785 :                               DO i = 1, ncoa
     249    81022810 :                                  p2 = (p1 + i - 1)*ncob
     250   204637780 :                                  DO j = 1, ncob
     251   101683270 :                                     p3 = p2 + j
     252   101683270 :                                     int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
     253   182706080 :                                     int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
     254             :                                  END DO
     255             :                               END DO
     256             :                            END DO
     257             :                         ELSE
     258      785845 :                            DO k = 1, ncoc
     259      606119 :                               p1 = (k - 1)*ncoa
     260     2922451 :                               DO i = 1, ncoa
     261     2136606 :                                  p2 = (p1 + i - 1)*ncob
     262     5275235 :                                  DO j = 1, ncob
     263     2532510 :                                     p3 = p2 + j
     264     4669116 :                                     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     2689329 :       DEALLOCATE (params)
     279             : 
     280     2689329 :    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    20446288 :    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    20446288 :       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    20446288 :       IF (PRESENT(params_in)) THEN
     328    10223144 :          params => params_in
     329             : 
     330             :          !If no external parameters to use, compute them
     331             :       ELSE
     332    10223144 :          params => params_out
     333             : 
     334             :          !Note: some variable of 4-center integrals simplify with a dummy center:
     335             :          !      P -> rk, gammap -> zetk
     336    10223144 :          params%m_max = li_max + lj_max + lk_max
     337    10223144 :          gammaq = zeti + zetj
     338    10223144 :          params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
     339    10223144 :          params%ZetapEtaInv = 1._dp/(zetk + gammaq)
     340             : 
     341    40892576 :          params%Q = (zeti*ri + zetj*rj)*params%EtaInv
     342    40892576 :          params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
     343    10223144 :          params%Rho = zetk*gammaq/(zetk + gammaq)
     344             : 
     345   224909168 :          params%Fm = 0.0_dp
     346    10223144 :          SELECT CASE (potential_parameter%potential_type)
     347             :          CASE (do_potential_coulomb)
     348     1455148 :             T = params%Rho*SUM((params%Q - rk)**2)
     349     1455148 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     350      363787 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     351             : 
     352      363787 :             CALL fgamma(params%m_max, T, params%Fm)
     353     8003314 :             params%Fm = prefac*params%Fm
     354             :          CASE (do_potential_truncated)
     355     6313163 :             R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     356    25252652 :             T = params%Rho*SUM((params%Q - rk)**2)
     357    25252652 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     358     6313163 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     359             : 
     360     6313163 :             CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     361     6313163 :             CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     362     6313163 :             IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     363   138889586 :             params%Fm = prefac*params%Fm
     364             :          CASE (do_potential_short)
     365     8439944 :             T = params%Rho*SUM((params%Q - rk)**2)
     366     8439944 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     367     2109986 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     368             : 
     369     2109986 :             CALL fgamma(params%m_max, T, params%Fm)
     370             : 
     371     2109986 :             omega2 = potential_parameter%omega**2
     372     2109986 :             omega_corr2 = omega2/(omega2 + params%Rho)
     373     2109986 :             omega_corr = SQRT(omega_corr2)
     374     2109986 :             T = T*omega_corr2
     375     2109986 :             ALLOCATE (Fm(prim_data_f_size))
     376             : 
     377     2109986 :             CALL fgamma(params%m_max, T, Fm)
     378     2109986 :             tmp = -omega_corr
     379    11707924 :             DO l = 1, params%m_max + 1
     380     9597938 :                params%Fm(l) = params%Fm(l) + Fm(l)*tmp
     381    11707924 :                tmp = tmp*omega_corr2
     382             :             END DO
     383    46419692 :             params%Fm = prefac*params%Fm
     384             :          CASE (do_potential_id)
     385             :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
     386    10053456 :                         - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
     387     1436208 :             prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
     388             : 
     389    31596576 :             params%Fm(:) = prefac
     390             :          CASE DEFAULT
     391    10223144 :             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    20446288 :                                     params%m_max, params%Fm)
     399             : 
     400    20446288 :    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      337191 :    SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
     438      337191 :                                  la_min, la_max, npgfa, zeta, rpgfa, ra, &
     439      337191 :                                  lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
     440      337191 :                                  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      337191 :       REAL(dp), DIMENSION(:, :), POINTER                 :: p_deriv
     466             :       TYPE(params_3c), POINTER                           :: params
     467             : 
     468      337191 :       NULLIFY (params, p_deriv)
     469    10115730 :       ALLOCATE (params)
     470             : 
     471      337191 :       permute_1 = [4, 5, 6]
     472      337191 :       permute_2 = [7, 8, 9]
     473             : 
     474      337191 :       dr_ab = 0.0_dp
     475      337191 :       dr_bc = 0.0_dp
     476      337191 :       dr_ac = 0.0_dp
     477             : 
     478      337191 :       op = potential_parameter%potential_type
     479             : 
     480      337191 :       IF (op == do_potential_truncated .OR. op == do_potential_short) THEN
     481      119807 :          dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
     482      119807 :          dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
     483      217384 :       ELSEIF (op == do_potential_coulomb) THEN
     484        9988 :          dr_bc = 1000000.0_dp
     485        9988 :          dr_ac = 1000000.0_dp
     486             :       END IF
     487             : 
     488      337191 :       do_ext = .FALSE.
     489      337191 :       IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
     490      337191 :       der_abc_1_ext_prv = 0.0_dp
     491      337191 :       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     1168188 :       DO ipgf = 1, npgfa
     499      830997 :          zeti = zeta(ipgf)
     500      830997 :          a_start = (ipgf - 1)*ncoset(la_max)
     501             : 
     502     4723948 :          DO jpgf = 1, npgfb
     503             : 
     504             :             ! screening
     505     3555760 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
     506             : 
     507     1121639 :             zetj = zetb(jpgf)
     508     1121639 :             b_start = (jpgf - 1)*ncoset(lb_max)
     509             : 
     510     8248366 :             DO kpgf = 1, npgfc
     511             : 
     512             :                ! screening
     513     6295730 :                IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
     514     2636184 :                IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
     515             : 
     516     1847000 :                zetk = zetc(kpgf)
     517     1847000 :                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     1847000 :                                         potential_parameter=potential_parameter, params_out=params)
     522             : 
     523     4256937 :                DO li = la_min, la_max
     524     2409937 :                   a_offset = a_start + ncoset(li - 1)
     525     2409937 :                   ncoa = nco(li)
     526     6831831 :                   DO lj = MAX(li, lb_min), lb_max
     527     2574894 :                      b_offset = b_start + ncoset(lj - 1)
     528     2574894 :                      ncob = nco(lj)
     529     8719891 :                      DO lk = lc_min, lc_max
     530     3735060 :                         c_offset = c_start + ncoset(lk - 1)
     531     3735060 :                         ncoc = nco(lk)
     532             : 
     533     3735060 :                         a_mysize(1) = ncoa*ncob*ncoc
     534             : 
     535     3735060 :                         CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
     536             : 
     537     3735060 :                         IF (do_ext) THEN
     538    14940240 :                            DO i_deriv = 1, 3
     539    42339804 :                               DO k = 1, ncoc
     540    27399564 :                                  p1 = (k - 1)*ncob
     541    90160152 :                                  DO j = 1, ncob
     542    51555408 :                                     p2 = (p1 + j - 1)*ncoa
     543   156637704 :                                     DO i = 1, ncoa
     544    77682732 :                                        p3 = p2 + i
     545             : 
     546             :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     547    77682732 :                                           p_deriv(p3, permute_2(i_deriv))
     548             :                                        der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
     549    77682732 :                                                                         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    77682732 :                                           p_deriv(p3, permute_1(i_deriv))
     553             :                                        der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
     554   129238140 :                                                                         ABS(p_deriv(p3, permute_1(i_deriv))))
     555             : 
     556             :                                     END DO
     557             :                                  END DO
     558             :                               END DO
     559             :                            END DO
     560             :                         ELSE
     561           0 :                            DO i_deriv = 1, 3
     562           0 :                               DO k = 1, ncoc
     563           0 :                                  p1 = (k - 1)*ncob
     564           0 :                                  DO j = 1, ncob
     565           0 :                                     p2 = (p1 + j - 1)*ncoa
     566           0 :                                     DO i = 1, ncoa
     567           0 :                                        p3 = p2 + i
     568             : 
     569             :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     570           0 :                                           p_deriv(p3, permute_2(i_deriv))
     571             : 
     572             :                                        der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     573           0 :                                           p_deriv(p3, permute_1(i_deriv))
     574             :                                     END DO
     575             :                                  END DO
     576             :                               END DO
     577             :                            END DO
     578             :                         END IF
     579             : 
     580     6309954 :                         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     1847000 :                CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
     587             : 
     588     7821822 :                DO lj = lb_min, lb_max
     589     2419062 :                   b_offset = b_start + ncoset(lj - 1)
     590     2419062 :                   ncob = nco(lj)
     591     9372505 :                   DO li = MAX(lj + 1, la_min), la_max
     592      657713 :                      a_offset = a_start + ncoset(li - 1)
     593      657713 :                      ncoa = nco(li)
     594     4043246 :                      DO lk = lc_min, lc_max
     595      966471 :                         c_offset = c_start + ncoset(lk - 1)
     596      966471 :                         ncoc = nco(lk)
     597             : 
     598      966471 :                         a_mysize(1) = ncoa*ncob*ncoc
     599      966471 :                         CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
     600             : 
     601      966471 :                         IF (do_ext) THEN
     602     3865884 :                            DO i_deriv = 1, 3
     603    11172267 :                               DO k = 1, ncoc
     604     7306383 :                                  p1 = (k - 1)*ncoa
     605    34168179 :                                  DO i = 1, ncoa
     606    23962383 :                                     p2 = (p1 + i - 1)*ncob
     607    58639773 :                                     DO j = 1, ncob
     608    27371007 :                                        p3 = p2 + j
     609             : 
     610             :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     611    27371007 :                                           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    27371007 :                                                                         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    27371007 :                                           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    51333390 :                                                                         ABS(p_deriv(p3, permute_2(i_deriv))))
     621             :                                     END DO
     622             :                                  END DO
     623             :                               END DO
     624             :                            END DO
     625             :                         ELSE
     626           0 :                            DO i_deriv = 1, 3
     627           0 :                               DO k = 1, ncoc
     628           0 :                                  p1 = (k - 1)*ncoa
     629           0 :                                  DO i = 1, ncoa
     630           0 :                                     p2 = (p1 + i - 1)*ncob
     631           0 :                                     DO j = 1, ncob
     632           0 :                                        p3 = p2 + j
     633             : 
     634             :                                        der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     635           0 :                                           p_deriv(p3, permute_1(i_deriv))
     636             : 
     637             :                                        der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
     638           0 :                                           p_deriv(p3, permute_2(i_deriv))
     639             :                                     END DO
     640             :                                  END DO
     641             :                               END DO
     642             :                            END DO
     643             :                         END IF
     644             : 
     645     1624184 :                         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      337191 :       IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
     655      337191 :       IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
     656             : 
     657      337191 :       DEALLOCATE (params)
     658             : 
     659      337191 :    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     3694000 :    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     3694000 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: Fm
     695             :       TYPE(params_3c), POINTER                           :: params
     696             : 
     697     3694000 :       IF (PRESENT(params_in)) THEN
     698     1847000 :          params => params_in
     699             : 
     700             :       ELSE
     701     1847000 :          params => params_out
     702             : 
     703     1847000 :          params%m_max = li_max + lj_max + lk_max + 1
     704     1847000 :          gammaq = zeti + zetj
     705     1847000 :          params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
     706     1847000 :          params%ZetapEtaInv = 1._dp/(zetk + gammaq)
     707             : 
     708     7388000 :          params%Q = (zeti*ri + zetj*rj)*params%EtaInv
     709     7388000 :          params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
     710     1847000 :          params%Rho = zetk*gammaq/(zetk + gammaq)
     711             : 
     712    40634000 :          params%Fm = 0.0_dp
     713     1847000 :          SELECT CASE (potential_parameter%potential_type)
     714             :          CASE (do_potential_coulomb)
     715      416424 :             T = params%Rho*SUM((params%Q - rk)**2)
     716      416424 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     717      104106 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     718             : 
     719      104106 :             CALL fgamma(params%m_max, T, params%Fm)
     720     2290332 :             params%Fm = prefac*params%Fm
     721             :          CASE (do_potential_truncated)
     722     1244984 :             R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     723     4979936 :             T = params%Rho*SUM((params%Q - rk)**2)
     724     4979936 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     725     1244984 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     726             : 
     727     1244984 :             CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     728     1244984 :             CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     729     1244984 :             IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     730    27389648 :             params%Fm = prefac*params%Fm
     731             :          CASE (do_potential_short)
     732           0 :             T = params%Rho*SUM((params%Q - rk)**2)
     733           0 :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
     734           0 :             prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
     735             : 
     736           0 :             CALL fgamma(params%m_max, T, params%Fm)
     737             : 
     738           0 :             omega2 = potential_parameter%omega**2
     739           0 :             omega_corr2 = omega2/(omega2 + params%Rho)
     740           0 :             omega_corr = SQRT(omega_corr2)
     741           0 :             T = T*omega_corr2
     742           0 :             ALLOCATE (Fm(prim_data_f_size))
     743             : 
     744           0 :             CALL fgamma(params%m_max, T, Fm)
     745           0 :             tmp = -omega_corr
     746           0 :             DO l = 1, params%m_max + 1
     747           0 :                params%Fm(l) = params%Fm(l) + Fm(l)*tmp
     748           0 :                tmp = tmp*omega_corr2
     749             :             END DO
     750           0 :             params%Fm = prefac*params%Fm
     751             :          CASE (do_potential_id)
     752             :             S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
     753     3485370 :                         - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
     754      497910 :             prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
     755             : 
     756    10954020 :             params%Fm(:) = prefac
     757             :          CASE DEFAULT
     758     1847000 :             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     3694000 :                                           params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
     766             : 
     767     3694000 :    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      387677 :    SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
     793      387677 :                           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      387677 :       REAL(dp), DIMENSION(:), POINTER                    :: p_work
     812             : 
     813      387677 :       NULLIFY (p_work)
     814             : 
     815      387677 :       dr_ab = 0.0_dp
     816             : 
     817      387677 :       IF (potential_parameter%potential_type == do_potential_truncated .OR. &
     818             :           potential_parameter%potential_type == do_potential_short) THEN
     819      223607 :          dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
     820      164070 :       ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
     821         984 :          dr_ab = 1000000.0_dp
     822             :       END IF
     823             : 
     824             :       !Looping over the pgfs
     825     1229208 :       DO ipgf = 1, npgfa
     826      841531 :          zeti = zeta(ipgf)
     827      841531 :          a_start = (ipgf - 1)*ncoset(la_max)
     828             : 
     829     5439004 :          DO jpgf = 1, npgfb
     830     4209796 :             zetj = zetb(jpgf)
     831     4209796 :             b_start = (jpgf - 1)*ncoset(lb_max)
     832             : 
     833             :             !screening
     834     4209796 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
     835             : 
     836     2791512 :             CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
     837             : 
     838     8829385 :             DO li = la_min, la_max
     839     5196342 :                a_offset = a_start + ncoset(li - 1)
     840     5196342 :                ncoa = nco(li)
     841    19725008 :                DO lj = lb_min, lb_max
     842    10318870 :                   b_offset = b_start + ncoset(lj - 1)
     843    10318870 :                   ncob = nco(lj)
     844             : 
     845    10318870 :                   a_mysize(1) = ncoa*ncob
     846    10318870 :                   CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
     847             : 
     848    47046018 :                   DO j = 1, ncob
     849    31530806 :                      p1 = (j - 1)*ncoa
     850   139972368 :                      DO i = 1, ncoa
     851    98122692 :                         p2 = p1 + i
     852   129653498 :                         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      387677 :    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     2791512 :    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     2791512 :       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     2791512 :       op = potential_parameter%potential_type
     900     2791512 :       params%m_max = lj_max + lk_max
     901     2791512 :       params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
     902     2791512 :       params%ZetapEtaInv = 1._dp/(zetk + zetj)
     903             : 
     904    11166048 :       params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
     905     2791512 :       params%Rho = zetk*zetj/(zetk + zetj)
     906             : 
     907    61413264 :       params%Fm = 0.0_dp
     908             :       SELECT CASE (op)
     909             :       CASE (do_potential_coulomb)
     910       94768 :          T = params%Rho*SUM((rj - rk)**2)
     911       23692 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
     912       23692 :          CALL fgamma(params%m_max, T, params%Fm)
     913      521224 :          params%Fm = prefac*params%Fm
     914             :       CASE (do_potential_truncated)
     915      173001 :          R = potential_parameter%cutoff_radius*SQRT(params%Rho)
     916      692004 :          T = params%Rho*SUM((rj - rk)**2)
     917      173001 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
     918             : 
     919      173001 :          CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
     920      173001 :          CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
     921      173001 :          IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
     922     3806022 :          params%Fm = prefac*params%Fm
     923             :       CASE (do_potential_short)
     924    10017192 :          T = params%Rho*SUM((rj - rk)**2)
     925     2504298 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
     926             : 
     927     2504298 :          CALL fgamma(params%m_max, T, params%Fm)
     928             : 
     929     2504298 :          omega2 = potential_parameter%omega**2
     930     2504298 :          omega_corr2 = omega2/(omega2 + params%Rho)
     931     2504298 :          omega_corr = SQRT(omega_corr2)
     932     2504298 :          T = T*omega_corr2
     933     2504298 :          ALLOCATE (Fm(prim_data_f_size))
     934             : 
     935     2504298 :          CALL fgamma(params%m_max, T, Fm)
     936     2504298 :          tmp = -omega_corr
     937    12305810 :          DO l = 1, params%m_max + 1
     938     9801512 :             params%Fm(l) = params%Fm(l) + Fm(l)*tmp
     939    12305810 :             tmp = tmp*omega_corr2
     940             :          END DO
     941    55094556 :          params%Fm = prefac*params%Fm
     942             :       CASE (do_potential_id)
     943             : 
     944      362084 :          prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
     945     1991462 :          params%Fm(:) = prefac
     946             :       CASE DEFAULT
     947     2791512 :          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     2791512 :                                     params%m_max, params%Fm)
     953             : 
     954    72579312 :    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       13552 :    PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
     963             :       TYPE(libint_potential_type), INTENT(IN)            :: potential1, potential2
     964             :       LOGICAL                                            :: equals
     965             : 
     966       13552 :       IF (potential1%potential_type /= potential2%potential_type) THEN
     967             :          equals = .FALSE.
     968             :       ELSE
     969       12776 :          equals = .TRUE.
     970         734 :          SELECT CASE (potential1%potential_type)
     971             :          CASE (do_potential_short, do_potential_long)
     972         734 :             IF (potential1%omega /= potential2%omega) equals = .FALSE.
     973             :          CASE (do_potential_truncated)
     974       12776 :             IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
     975             :          END SELECT
     976             :       END IF
     977             : 
     978       13552 :    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      137567 :    SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
    1003      137567 :                                  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      137567 :       REAL(dp), DIMENSION(:, :), POINTER                 :: p_deriv
    1023             : 
    1024      137567 :       NULLIFY (p_deriv)
    1025             : 
    1026      137567 :       permute = [4, 5, 6]
    1027             : 
    1028      137567 :       dr_ab = 0.0_dp
    1029             : 
    1030      137567 :       IF (potential_parameter%potential_type == do_potential_truncated .OR. &
    1031             :           potential_parameter%potential_type == do_potential_short) THEN
    1032       61518 :          dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
    1033       76049 :       ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
    1034        9831 :          dr_ab = 1000000.0_dp
    1035             :       END IF
    1036             : 
    1037             :       !Looping over the pgfs
    1038      620990 :       DO ipgf = 1, npgfa
    1039      483423 :          zeti = zeta(ipgf)
    1040      483423 :          a_start = (ipgf - 1)*ncoset(la_max)
    1041             : 
    1042     3811929 :          DO jpgf = 1, npgfb
    1043     3190939 :             zetj = zetb(jpgf)
    1044     3190939 :             b_start = (jpgf - 1)*ncoset(lb_max)
    1045             : 
    1046             :             !screening
    1047     3190939 :             IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
    1048             : 
    1049     2166114 :             CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
    1050             : 
    1051     6287730 :             DO li = la_min, la_max
    1052     3638193 :                a_offset = a_start + ncoset(li - 1)
    1053     3638193 :                ncoa = nco(li)
    1054    12986663 :                DO lj = lb_min, lb_max
    1055     6157531 :                   b_offset = b_start + ncoset(lj - 1)
    1056     6157531 :                   ncob = nco(lj)
    1057             : 
    1058     6157531 :                   a_mysize(1) = ncoa*ncob
    1059     6157531 :                   CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
    1060             : 
    1061    24630124 :                   DO i_deriv = 1, 3
    1062    75239905 :                      DO j = 1, ncob
    1063    50609781 :                         p1 = (j - 1)*ncoa
    1064   208172862 :                         DO i = 1, ncoa
    1065   139090488 :                            p2 = p1 + i
    1066   189700269 :                            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     9795724 :                   DEALLOCATE (p_deriv)
    1072             :                END DO
    1073             :             END DO
    1074             : 
    1075             :          END DO
    1076             :       END DO
    1077             : 
    1078      137567 :    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     2166114 :    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     2166114 :       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     2166114 :       op = potential_parameter%potential_type
    1116     2166114 :       params%m_max = lj_max + lk_max + 1
    1117     2166114 :       params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
    1118     2166114 :       params%ZetapEtaInv = 1._dp/(zetk + zetj)
    1119             : 
    1120     8664456 :       params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
    1121     2166114 :       params%Rho = zetk*zetj/(zetk + zetj)
    1122             : 
    1123    47654508 :       params%Fm = 0.0_dp
    1124             :       SELECT CASE (op)
    1125             :       CASE (do_potential_coulomb)
    1126      102916 :          T = params%Rho*SUM((rj - rk)**2)
    1127       25729 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1128       25729 :          CALL fgamma(params%m_max, T, params%Fm)
    1129      566038 :          params%Fm = prefac*params%Fm
    1130             :       CASE (do_potential_truncated)
    1131       60843 :          R = potential_parameter%cutoff_radius*SQRT(params%Rho)
    1132      243372 :          T = params%Rho*SUM((rj - rk)**2)
    1133       60843 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1134             : 
    1135       60843 :          CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
    1136       60843 :          CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
    1137       60843 :          IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
    1138     1338546 :          params%Fm = prefac*params%Fm
    1139             :       CASE (do_potential_short)
    1140     8096600 :          T = params%Rho*SUM((rj - rk)**2)
    1141     2024150 :          prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
    1142             : 
    1143     2024150 :          CALL fgamma(params%m_max, T, params%Fm)
    1144             : 
    1145     2024150 :          omega2 = potential_parameter%omega**2
    1146     2024150 :          omega_corr2 = omega2/(omega2 + params%Rho)
    1147     2024150 :          omega_corr = SQRT(omega_corr2)
    1148     2024150 :          T = T*omega_corr2
    1149     2024150 :          ALLOCATE (Fm(prim_data_f_size))
    1150             : 
    1151     2024150 :          CALL fgamma(params%m_max, T, Fm)
    1152     2024150 :          tmp = -omega_corr
    1153    11369176 :          DO l = 1, params%m_max + 1
    1154     9345026 :             params%Fm(l) = params%Fm(l) + Fm(l)*tmp
    1155    11369176 :             tmp = tmp*omega_corr2
    1156             :          END DO
    1157    44531300 :          params%Fm = prefac*params%Fm
    1158             :       CASE (do_potential_id)
    1159             : 
    1160      221568 :          prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
    1161     1218624 :          params%Fm(:) = prefac
    1162             :       CASE DEFAULT
    1163     2166114 :          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     2166114 :                                           params%m_max, params%Fm)
    1170             : 
    1171    56318964 :    END SUBROUTINE set_params_2c_deriv
    1172             : 
    1173           0 : END MODULE libint_2c_3c
    1174             : 

Generated by: LCOV version 1.15