LCOV - code coverage report
Current view: top level - src - hfx_libint_interface.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:70636b1) Lines: 96.7 % 983 951
Test Date: 2026-02-11 07:00:35 Functions: 100.0 % 6 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Interface to the Libint-Library
      10              : !> \par History
      11              : !>      11.2006 created [Manuel Guidon]
      12              : !>      11.2019 Fixed potential_id initial values (A. Bussy)
      13              : !> \author Manuel Guidon
      14              : ! **************************************************************************************************
      15              : MODULE hfx_libint_interface
      16              : 
      17              :    USE cell_types,                      ONLY: cell_type,&
      18              :                                               real_to_scaled
      19              :    USE gamma,                           ONLY: fgamma => fgamma_0
      20              :    USE hfx_contraction_methods,         ONLY: contract
      21              :    USE hfx_types,                       ONLY: hfx_pgf_product_list,&
      22              :                                               hfx_potential_type
      23              :    USE input_constants,                 ONLY: &
      24              :         do_potential_coulomb, do_potential_gaussian, do_potential_id, do_potential_long, &
      25              :         do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, &
      26              :         do_potential_truncated
      27              :    USE kinds,                           ONLY: dp,&
      28              :                                               int_8
      29              :    USE libint_wrapper,                  ONLY: cp_libint_get_derivs,&
      30              :                                               cp_libint_get_eris,&
      31              :                                               cp_libint_set_params_eri,&
      32              :                                               cp_libint_set_params_eri_deriv,&
      33              :                                               cp_libint_set_params_eri_screen,&
      34              :                                               cp_libint_t,&
      35              :                                               get_ssss_f_val,&
      36              :                                               prim_data_f_size
      37              :    USE mathconstants,                   ONLY: pi
      38              :    USE orbital_pointers,                ONLY: nco
      39              :    USE t_c_g0,                          ONLY: t_c_g0_n
      40              : #include "./base/base_uses.f90"
      41              : 
      42              :    IMPLICIT NONE
      43              :    PRIVATE
      44              :    PUBLIC :: evaluate_eri, &
      45              :              evaluate_deriv_eri, &
      46              :              evaluate_eri_screen
      47              : 
      48              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm1 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
      49              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm2 = [4, 5, 6, 1, 2, 3, 7, 8, 9, 10, 11, 12]
      50              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm3 = [1, 2, 3, 4, 5, 6, 10, 11, 12, 7, 8, 9]
      51              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm4 = [4, 5, 6, 1, 2, 3, 10, 11, 12, 7, 8, 9]
      52              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm5 = [7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6]
      53              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm6 = [7, 8, 9, 10, 11, 12, 4, 5, 6, 1, 2, 3]
      54              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm7 = [10, 11, 12, 7, 8, 9, 1, 2, 3, 4, 5, 6]
      55              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm8 = [10, 11, 12, 7, 8, 9, 4, 5, 6, 1, 2, 3]
      56              : 
      57              : !***
      58              : 
      59              : CONTAINS
      60              : 
      61              : ! **************************************************************************************************
      62              : !> \brief Fill data structure used in libint
      63              : !> \param A ...
      64              : !> \param B ...
      65              : !> \param C ...
      66              : !> \param D ...
      67              : !> \param Zeta_A ...
      68              : !> \param Zeta_B ...
      69              : !> \param Zeta_C ...
      70              : !> \param Zeta_D ...
      71              : !> \param m_max ...
      72              : !> \param potential_parameter ...
      73              : !> \param libint ...
      74              : !> \param R11 ...
      75              : !> \param R22 ...
      76              : !> \par History
      77              : !>      03.2007 created [Manuel Guidon]
      78              : !> \author Manuel Guidon
      79              : ! **************************************************************************************************
      80     87276120 :    SUBROUTINE build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, &
      81              :                                         potential_parameter, libint, R11, R22)
      82              : 
      83              :       REAL(KIND=dp)                                      :: A(3), B(3), C(3), D(3)
      84              :       REAL(KIND=dp), INTENT(IN)                          :: Zeta_A, Zeta_B, Zeta_C, Zeta_D
      85              :       INTEGER, INTENT(IN)                                :: m_max
      86              :       TYPE(hfx_potential_type)                           :: potential_parameter
      87              :       TYPE(cp_libint_t)                                  :: libint
      88              :       REAL(dp)                                           :: R11, R22
      89              : 
      90              :       INTEGER                                            :: i
      91              :       LOGICAL                                            :: use_gamma
      92              :       REAL(KIND=dp) :: AB(3), AB2, CD(3), CD2, den, Eta, EtaInv, factor, G(3), num, omega2, &
      93              :          omega_corr, omega_corr2, P(3), PQ(3), PQ2, Q(3), R, R1, R2, Rho, RhoInv, S1234, ssss, T, &
      94              :          tmp, W(3), Zeta, ZetaInv, ZetapEtaInv
      95              :       REAL(KIND=dp), DIMENSION(prim_data_f_size)         :: F, Fm
      96              : 
      97     87276120 :       Zeta = Zeta_A + Zeta_B
      98     87276120 :       ZetaInv = 1.0_dp/Zeta
      99     87276120 :       Eta = Zeta_C + Zeta_D
     100     87276120 :       EtaInv = 1.0_dp/Eta
     101     87276120 :       ZetapEtaInv = Zeta + Eta
     102     87276120 :       ZetapEtaInv = 1.0_dp/ZetapEtaInv
     103     87276120 :       Rho = Zeta*Eta*ZetapEtaInv
     104     87276120 :       RhoInv = 1.0_dp/Rho
     105              : 
     106    349104480 :       DO i = 1, 3
     107    261828360 :          P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv
     108    261828360 :          Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv
     109    261828360 :          AB(i) = A(i) - B(i)
     110    261828360 :          CD(i) = C(i) - D(i)
     111    261828360 :          PQ(i) = P(i) - Q(i)
     112    349104480 :          W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv
     113              :       END DO
     114              : 
     115    349104480 :       AB2 = DOT_PRODUCT(AB, AB)
     116    349104480 :       CD2 = DOT_PRODUCT(CD, CD)
     117    349104480 :       PQ2 = DOT_PRODUCT(PQ, PQ)
     118              : 
     119     87276120 :       S1234 = EXP((-Zeta_A*Zeta_B*ZetaInv*AB2) + (-Zeta_C*Zeta_D*EtaInv*CD2))
     120     87276120 :       T = Rho*PQ2
     121              : 
     122    100482072 :       SELECT CASE (potential_parameter%potential_type)
     123              :       CASE (do_potential_truncated)
     124     13205952 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     125     13205952 :          R1 = R11
     126     13205952 :          R2 = R22
     127     13205952 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     128            0 :             RETURN
     129              :          END IF
     130     13205952 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     131     13205952 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     132     13205952 :          factor = 2.0_dp*Pi*RhoInv
     133              :       CASE (do_potential_coulomb)
     134     66967848 :          CALL fgamma(m_max, T, F(1))
     135     66967848 :          factor = 2.0_dp*Pi*RhoInv
     136              :       CASE (do_potential_short)
     137      3091812 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     138      3091812 :          R1 = R11
     139      3091812 :          R2 = R22
     140      3091812 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     141              :             RETURN
     142              :          END IF
     143      3091812 :          CALL fgamma(m_max, T, F(1))
     144      3091812 :          omega2 = potential_parameter%omega**2
     145      3091812 :          omega_corr2 = omega2/(omega2 + Rho)
     146      3091812 :          omega_corr = SQRT(omega_corr2)
     147      3091812 :          T = T*omega_corr2
     148      3091812 :          CALL fgamma(m_max, T, Fm)
     149      3091812 :          tmp = -omega_corr
     150     11625504 :          DO i = 1, m_max + 1
     151      8533692 :             F(i) = F(i) + Fm(i)*tmp
     152     11625504 :             tmp = tmp*omega_corr2
     153              :          END DO
     154      3091812 :          factor = 2.0_dp*Pi*RhoInv
     155              :       CASE (do_potential_long)
     156       978084 :          omega2 = potential_parameter%omega**2
     157       978084 :          omega_corr2 = omega2/(omega2 + Rho)
     158       978084 :          omega_corr = SQRT(omega_corr2)
     159       978084 :          T = T*omega_corr2
     160       978084 :          CALL fgamma(m_max, T, F(1))
     161       978084 :          tmp = omega_corr
     162      3565704 :          DO i = 1, m_max + 1
     163      2587620 :             F(i) = F(i)*tmp
     164      3565704 :             tmp = tmp*omega_corr2
     165              :          END DO
     166       978084 :          factor = 2.0_dp*Pi*RhoInv
     167              :       CASE (do_potential_mix_cl)
     168      1558026 :          CALL fgamma(m_max, T, F(1))
     169      1558026 :          omega2 = potential_parameter%omega**2
     170      1558026 :          omega_corr2 = omega2/(omega2 + Rho)
     171      1558026 :          omega_corr = SQRT(omega_corr2)
     172      1558026 :          T = T*omega_corr2
     173      1558026 :          CALL fgamma(m_max, T, Fm)
     174      1558026 :          tmp = omega_corr
     175      5416428 :          DO i = 1, m_max + 1
     176      3858402 :             F(i) = F(i)*potential_parameter%scale_coulomb + Fm(i)*tmp*potential_parameter%scale_longrange
     177      5416428 :             tmp = tmp*omega_corr2
     178              :          END DO
     179      1558026 :          factor = 2.0_dp*Pi*RhoInv
     180              :       CASE (do_potential_mix_cl_trunc)
     181              : 
     182              :          ! truncated
     183       959298 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     184       959298 :          R1 = R11
     185       959298 :          R2 = R22
     186       959298 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     187              :             RETURN
     188              :          END IF
     189       959298 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     190       959298 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     191              : 
     192              :          ! Coulomb
     193       959298 :          CALL fgamma(m_max, T, Fm)
     194              : 
     195      3419052 :          DO i = 1, m_max + 1
     196              :             F(i) = F(i)*(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
     197      3419052 :                    Fm(i)*potential_parameter%scale_longrange
     198              :          END DO
     199              : 
     200              :          ! longrange
     201       959298 :          omega2 = potential_parameter%omega**2
     202       959298 :          omega_corr2 = omega2/(omega2 + Rho)
     203       959298 :          omega_corr = SQRT(omega_corr2)
     204       959298 :          T = T*omega_corr2
     205       959298 :          CALL fgamma(m_max, T, Fm)
     206       959298 :          tmp = omega_corr
     207      3419052 :          DO i = 1, m_max + 1
     208      2459754 :             F(i) = F(i) + Fm(i)*tmp*potential_parameter%scale_longrange
     209      3419052 :             tmp = tmp*omega_corr2
     210              :          END DO
     211       959298 :          factor = 2.0_dp*Pi*RhoInv
     212              : 
     213              :       CASE (do_potential_gaussian)
     214            0 :          omega2 = potential_parameter%omega**2
     215            0 :          T = -omega2*T/(Rho + omega2)
     216            0 :          tmp = 1.0_dp
     217            0 :          DO i = 1, m_max + 1
     218            0 :             F(i) = EXP(T)*tmp
     219            0 :             tmp = tmp*omega2/(Rho + omega2)
     220              :          END DO
     221            0 :          factor = (Pi/(Rho + omega2))**(1.5_dp)
     222              :       CASE (do_potential_mix_lg)
     223       272700 :          omega2 = potential_parameter%omega**2
     224       272700 :          omega_corr2 = omega2/(omega2 + Rho)
     225       272700 :          omega_corr = SQRT(omega_corr2)
     226       272700 :          T = T*omega_corr2
     227       272700 :          CALL fgamma(m_max, T, Fm)
     228       272700 :          tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange
     229       981720 :          DO i = 1, m_max + 1
     230       709020 :             Fm(i) = Fm(i)*tmp
     231       981720 :             tmp = tmp*omega_corr2
     232              :          END DO
     233              :          T = Rho*PQ2
     234       272700 :          T = -omega2*T/(Rho + omega2)
     235       272700 :          tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
     236       981720 :          DO i = 1, m_max + 1
     237       709020 :             F(i) = EXP(T)*tmp + Fm(i)
     238       981720 :             tmp = tmp*omega2/(Rho + omega2)
     239              :          END DO
     240       242400 :          factor = 1.0_dp
     241              :       CASE (do_potential_id)
     242       242400 :          ssss = -Zeta_A*Zeta_B*ZetaInv*AB2
     243              : 
     244       242400 :          num = (Zeta_A + Zeta_B)*Zeta_C
     245       242400 :          den = Zeta_A + Zeta_B + Zeta_C
     246       969600 :          ssss = ssss - num/den*SUM((P - C)**2)
     247              : 
     248       969600 :          G(:) = (Zeta_A*A(:) + Zeta_B*B(:) + Zeta_C*C(:))/den
     249       242400 :          num = den*Zeta_D
     250       242400 :          den = den + Zeta_D
     251       969600 :          ssss = ssss - num/den*SUM((G - D)**2)
     252              : 
     253      5332800 :          F(:) = EXP(ssss)
     254       242400 :          factor = 1.0_dp
     255    101683770 :          IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
     256              :       END SELECT
     257              : 
     258     87276120 :       tmp = (Pi*ZetapEtaInv)**3
     259     87276120 :       factor = factor*S1234*SQRT(tmp)
     260              : 
     261    337418376 :       DO i = 1, m_max + 1
     262    337418376 :          F(i) = F(i)*factor
     263              :       END DO
     264              : 
     265     87276120 :       CALL cp_libint_set_params_eri_screen(libint, A, B, C, D, P, Q, W, ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F)
     266              : 
     267              :    END SUBROUTINE build_quartet_data_screen
     268              : 
     269              : ! **************************************************************************************************
     270              : !> \brief Fill data structure used in libderiv
     271              : !> \param libint ...
     272              : !> \param A ...
     273              : !> \param B ...
     274              : !> \param C ...
     275              : !> \param D ...
     276              : !> \param Zeta_A ...
     277              : !> \param Zeta_B ...
     278              : !> \param Zeta_C ...
     279              : !> \param Zeta_D ...
     280              : !> \param ZetaInv ...
     281              : !> \param EtaInv ...
     282              : !> \param ZetapEtaInv ...
     283              : !> \param Rho ...
     284              : !> \param RhoInv ...
     285              : !> \param P ...
     286              : !> \param Q ...
     287              : !> \param W ...
     288              : !> \param m_max ...
     289              : !> \param F ...
     290              : !> \par History
     291              : !>      03.2007 created [Manuel Guidon]
     292              : !> \author Manuel Guidon
     293              : ! **************************************************************************************************
     294    101110749 :    SUBROUTINE build_deriv_data(libint, A, B, C, D, &
     295              :                                Zeta_A, Zeta_B, Zeta_C, Zeta_D, &
     296              :                                ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, &
     297    101110749 :                                P, Q, W, m_max, F)
     298              : 
     299              :       TYPE(cp_libint_t)                                  :: libint
     300              :       REAL(KIND=dp), INTENT(IN)                          :: A(3), B(3), C(3), D(3)
     301              :       REAL(dp), INTENT(IN)                               :: Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, &
     302              :                                                             EtaInv, ZetapEtaInv, Rho, RhoInv, &
     303              :                                                             P(3), Q(3), W(3)
     304              :       INTEGER, INTENT(IN)                                :: m_max
     305              :       REAL(KIND=dp), DIMENSION(:)                        :: F
     306              : 
     307              :       MARK_USED(D) ! todo: fix
     308              :       MARK_USED(Zeta_C)
     309              :       MARK_USED(RhoInv)
     310              : 
     311              :       CALL cp_libint_set_params_eri_deriv(libint, A, B, C, D, P, Q, W, zeta_A, zeta_B, zeta_C, zeta_D, &
     312    101110749 :                                           ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F)
     313              : 
     314    101110749 :    END SUBROUTINE build_deriv_data
     315              : 
     316              : ! **************************************************************************************************
     317              : !> \brief Evaluates derivatives of electron repulsion integrals for a primitive quartet
     318              : !> \param deriv ...
     319              : !> \param nproducts ...
     320              : !> \param pgf_product_list ...
     321              : !> \param n_a ...
     322              : !> \param n_b ...
     323              : !> \param n_c ...
     324              : !> \param n_d ...
     325              : !> \param ncoa ...
     326              : !> \param ncob ...
     327              : !> \param ncoc ...
     328              : !> \param ncod ...
     329              : !> \param nsgfa ...
     330              : !> \param nsgfb ...
     331              : !> \param nsgfc ...
     332              : !> \param nsgfd ...
     333              : !> \param primitives ...
     334              : !> \param max_contraction ...
     335              : !> \param tmp_max_all ...
     336              : !> \param eps_schwarz ...
     337              : !> \param neris ...
     338              : !> \param Zeta_A ...
     339              : !> \param Zeta_B ...
     340              : !> \param Zeta_C ...
     341              : !> \param Zeta_D ...
     342              : !> \param ZetaInv ...
     343              : !> \param EtaInv ...
     344              : !> \param s_offset_a ...
     345              : !> \param s_offset_b ...
     346              : !> \param s_offset_c ...
     347              : !> \param s_offset_d ...
     348              : !> \param nl_a ...
     349              : !> \param nl_b ...
     350              : !> \param nl_c ...
     351              : !> \param nl_d ...
     352              : !> \param nsoa ...
     353              : !> \param nsob ...
     354              : !> \param nsoc ...
     355              : !> \param nsod ...
     356              : !> \param sphi_a ...
     357              : !> \param sphi_b ...
     358              : !> \param sphi_c ...
     359              : !> \param sphi_d ...
     360              : !> \param work ...
     361              : !> \param work2 ...
     362              : !> \param work_forces ...
     363              : !> \param buffer1 ...
     364              : !> \param buffer2 ...
     365              : !> \param primitives_tmp ...
     366              : !> \param use_virial ...
     367              : !> \param work_virial ...
     368              : !> \param work2_virial ...
     369              : !> \param primitives_tmp_virial ...
     370              : !> \param primitives_virial ...
     371              : !> \param cell ...
     372              : !> \param tmp_max_all_virial ...
     373              : !> \par History
     374              : !>      03.2007 created [Manuel Guidon]
     375              : !>      08.2007 refactured permutation part [Manuel Guidon]
     376              : !> \author Manuel Guidon
     377              : ! **************************************************************************************************
     378     34308234 :    SUBROUTINE evaluate_deriv_eri(deriv, nproducts, pgf_product_list, &
     379              :                                  n_a, n_b, n_c, n_d, &
     380              :                                  ncoa, ncob, ncoc, ncod, &
     381              :                                  nsgfa, nsgfb, nsgfc, nsgfd, &
     382     34308234 :                                  primitives, max_contraction, tmp_max_all, &
     383              :                                  eps_schwarz, neris, &
     384              :                                  Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, EtaInv, &
     385              :                                  s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
     386              :                                  nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, &
     387     34308234 :                                  sphi_a, sphi_b, sphi_c, sphi_d, &
     388     34308234 :                                  work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
     389     34308234 :                                  use_virial, work_virial, work2_virial, primitives_tmp_virial, &
     390     34308234 :                                  primitives_virial, cell, tmp_max_all_virial)
     391              : 
     392              :       TYPE(cp_libint_t)                                  :: deriv
     393              :       INTEGER, INTENT(IN)                                :: nproducts
     394              :       TYPE(hfx_pgf_product_list), DIMENSION(nproducts)   :: pgf_product_list
     395              :       INTEGER, INTENT(IN)                                :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, &
     396              :                                                             ncod, nsgfa, nsgfb, nsgfc, nsgfd
     397              :       REAL(dp), &
     398              :          DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12)       :: primitives
     399              :       REAL(dp)                                           :: max_contraction, tmp_max_all, eps_schwarz
     400              :       INTEGER(int_8)                                     :: neris
     401              :       REAL(dp), INTENT(IN)                               :: Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, &
     402              :                                                             EtaInv
     403              :       INTEGER                                            :: s_offset_a, s_offset_b, s_offset_c, &
     404              :                                                             s_offset_d, nl_a, nl_b, nl_c, nl_d, &
     405              :                                                             nsoa, nsob, nsoc, nsod
     406              :       REAL(dp), DIMENSION(ncoa, nsoa*nl_a)               :: sphi_a
     407              :       REAL(dp), DIMENSION(ncob, nsob*nl_b)               :: sphi_b
     408              :       REAL(dp), DIMENSION(ncoc, nsoc*nl_c)               :: sphi_c
     409              :       REAL(dp), DIMENSION(ncod, nsod*nl_d)               :: sphi_d
     410              :       REAL(dp) :: work(ncoa*ncob*ncoc*ncod, 12), work2(ncoa, ncob, ncoc, ncod, 12), &
     411              :          work_forces(ncoa*ncob*ncoc*ncod, 12)
     412              :       REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod)           :: buffer1, buffer2
     413              :       REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*&
     414              :          nl_c, nsod*nl_d)                                :: primitives_tmp
     415              :       LOGICAL, INTENT(IN)                                :: use_virial
     416              :       REAL(dp) :: work_virial(ncoa*ncob*ncoc*ncod, 12, 3), &
     417              :          work2_virial(ncoa, ncob, ncoc, ncod, 12, 3)
     418              :       REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*&
     419              :          nl_c, nsod*nl_d)                                :: primitives_tmp_virial
     420              :       REAL(dp), &
     421              :          DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12, 3)    :: primitives_virial
     422              :       TYPE(cell_type), POINTER                           :: cell
     423              :       REAL(dp)                                           :: tmp_max_all_virial
     424              : 
     425              :       INTEGER                                            :: a_mysize(1), i, j, k, l, m, m_max, &
     426              :                                                             mysize, n, p1, p2, p3, p4, perm_case
     427              :       REAL(dp)                                           :: A(3), AB(3), B(3), C(3), CD(3), D(3), &
     428              :                                                             P(3), Q(3), Rho, RhoInv, scoord(12), &
     429              :                                                             tmp_max, tmp_max_virial, W(3), &
     430              :                                                             ZetapEtaInv
     431              : 
     432     34308234 :       m_max = n_a + n_b + n_c + n_d
     433     34308234 :       m_max = m_max + 1
     434     34308234 :       mysize = ncoa*ncob*ncoc*ncod
     435     68616468 :       a_mysize = mysize
     436              : 
     437   4812837426 :       work = 0.0_dp
     438   9149440434 :       work2 = 0.0_dp
     439              : 
     440     34308234 :       IF (use_virial) THEN
     441    850401248 :          work_virial = 0.0_dp
     442   1465200824 :          work2_virial = 0.0_dp
     443              :       END IF
     444              : 
     445     34308234 :       perm_case = 1
     446     34308234 :       IF (n_a < n_b) THEN
     447      3319354 :          perm_case = perm_case + 1
     448              :       END IF
     449     34308234 :       IF (n_c < n_d) THEN
     450      5631476 :          perm_case = perm_case + 2
     451              :       END IF
     452     34308234 :       IF (n_a + n_b > n_c + n_d) THEN
     453      6346968 :          perm_case = perm_case + 4
     454              :       END IF
     455              : 
     456              :       SELECT CASE (perm_case)
     457              :       CASE (1)
     458     83151984 :          DO i = 1, nproducts
     459    245619496 :             A = pgf_product_list(i)%ra
     460    245619496 :             B = pgf_product_list(i)%rb
     461    245619496 :             C = pgf_product_list(i)%rc
     462    245619496 :             D = pgf_product_list(i)%rd
     463     61404874 :             Rho = pgf_product_list(i)%Rho
     464     61404874 :             RhoInv = pgf_product_list(i)%RhoInv
     465    245619496 :             P = pgf_product_list(i)%P
     466    245619496 :             Q = pgf_product_list(i)%Q
     467    245619496 :             W = pgf_product_list(i)%W
     468    245619496 :             AB = pgf_product_list(i)%AB
     469    245619496 :             CD = pgf_product_list(i)%CD
     470     61404874 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
     471              : 
     472              :             CALL build_deriv_data(deriv, A, B, C, D, &
     473              :                                   Zeta_A, Zeta_B, Zeta_C, Zeta_D, &
     474              :                                   ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, &
     475     61404874 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     476     61404874 :             CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize)
     477    245619496 :             DO k = 4, 6
     478   1920047443 :                DO j = 1, mysize
     479              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     480              :                                                work_forces(j, k + 3) + &
     481   1858642569 :                                                work_forces(j, k + 6))
     482              :                END DO
     483              :             END DO
     484    798263362 :             DO k = 1, 12
     485   7495975150 :                DO j = 1, mysize
     486   7434570276 :                   work(j, k) = work(j, k) + work_forces(j, k)
     487              :                END DO
     488              :             END DO
     489     61404874 :             neris = neris + 12*mysize
     490     83151984 :             IF (use_virial) THEN
     491     10723386 :                CALL real_to_scaled(scoord(1:3), A, cell)
     492     10723386 :                CALL real_to_scaled(scoord(4:6), B, cell)
     493     10723386 :                CALL real_to_scaled(scoord(7:9), C, cell)
     494     10723386 :                CALL real_to_scaled(scoord(10:12), D, cell)
     495    139404018 :                DO k = 1, 12
     496   1164286866 :                   DO j = 1, mysize
     497   4228212024 :                      DO m = 1, 3
     498   4099531392 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     499              :                      END DO
     500              :                   END DO
     501              :                END DO
     502              :             END IF
     503              :          END DO
     504              : 
     505    282712430 :          DO n = 1, 12
     506              :             tmp_max = 0.0_dp
     507   2345558820 :             DO i = 1, mysize
     508   2345558820 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     509              :             END DO
     510    260965320 :             tmp_max = tmp_max*max_contraction
     511    260965320 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     512              : 
     513    648126902 :             DO i = 1, ncoa
     514    365414472 :                p1 = (i - 1)*ncob
     515   1031352504 :                DO j = 1, ncob
     516    404972712 :                   p2 = (p1 + j - 1)*ncoc
     517   1880542692 :                   DO k = 1, ncoc
     518   1110155508 :                      p3 = (p2 + k - 1)*ncod
     519   3599721720 :                      DO l = 1, ncod
     520   2084593500 :                         p4 = p3 + l
     521   3194749008 :                         work2(i, j, k, l, full_perm1(n)) = work(p4, n)
     522              :                      END DO
     523              :                   END DO
     524              :                END DO
     525              :             END DO
     526              :          END DO
     527     21747110 :          IF (use_virial) THEN
     528      6233656 :             DO n = 1, 12
     529              :                tmp_max_virial = 0.0_dp
     530    112731564 :                DO i = 1, mysize
     531              :                   tmp_max_virial = MAX(tmp_max_virial, &
     532    112731564 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     533              :                END DO
     534      5754144 :                tmp_max_virial = tmp_max_virial*max_contraction
     535      5754144 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     536              : 
     537     16527736 :                DO i = 1, ncoa
     538     10294080 :                   p1 = (i - 1)*ncob
     539     29965200 :                   DO j = 1, ncob
     540     13916976 :                      p2 = (p1 + j - 1)*ncoc
     541     69604716 :                      DO k = 1, ncoc
     542     45393660 :                         p3 = (p2 + k - 1)*ncod
     543    166288056 :                         DO l = 1, ncod
     544    106977420 :                            p4 = p3 + l
     545    473303340 :                            work2_virial(i, j, k, l, full_perm1(n), 1:3) = work_virial(p4, n, 1:3)
     546              :                         END DO
     547              :                      END DO
     548              :                   END DO
     549              :                END DO
     550              :             END DO
     551              :          END IF
     552              :       CASE (2)
     553      6029727 :          DO i = 1, nproducts
     554     18808960 :             A = pgf_product_list(i)%ra
     555     18808960 :             B = pgf_product_list(i)%rb
     556     18808960 :             C = pgf_product_list(i)%rc
     557     18808960 :             D = pgf_product_list(i)%rd
     558      4702240 :             Rho = pgf_product_list(i)%Rho
     559      4702240 :             RhoInv = pgf_product_list(i)%RhoInv
     560     18808960 :             P = pgf_product_list(i)%P
     561     18808960 :             Q = pgf_product_list(i)%Q
     562     18808960 :             W = pgf_product_list(i)%W
     563     18808960 :             AB = pgf_product_list(i)%AB
     564     18808960 :             CD = pgf_product_list(i)%CD
     565      4702240 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
     566              : 
     567              :             CALL build_deriv_data(deriv, B, A, C, D, &
     568              :                                   Zeta_B, Zeta_A, Zeta_C, Zeta_D, &
     569              :                                   ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, &
     570      4702240 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     571      4702240 :             CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize)
     572     18808960 :             DO k = 4, 6
     573    339547657 :                DO j = 1, mysize
     574              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     575              :                                                work_forces(j, k + 3) + &
     576    334845417 :                                                work_forces(j, k + 6))
     577              :                END DO
     578              :             END DO
     579     61129120 :             DO k = 1, 12
     580   1344083908 :                DO j = 1, mysize
     581   1339381668 :                   work(j, k) = work(j, k) + work_forces(j, k)
     582              :                END DO
     583              :             END DO
     584      4702240 :             neris = neris + 12*mysize
     585      6029727 :             IF (use_virial) THEN
     586       674510 :                CALL real_to_scaled(scoord(1:3), B, cell)
     587       674510 :                CALL real_to_scaled(scoord(4:6), A, cell)
     588       674510 :                CALL real_to_scaled(scoord(7:9), C, cell)
     589       674510 :                CALL real_to_scaled(scoord(10:12), D, cell)
     590      8768630 :                DO k = 1, 12
     591    180406874 :                   DO j = 1, mysize
     592    694647096 :                      DO m = 1, 3
     593    686552976 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     594              :                      END DO
     595              :                   END DO
     596              :                END DO
     597              :             END IF
     598              : 
     599              :          END DO
     600     17257331 :          DO n = 1, 12
     601              :             tmp_max = 0.0_dp
     602    430611612 :             DO i = 1, mysize
     603    430611612 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     604              :             END DO
     605     15929844 :             tmp_max = tmp_max*max_contraction
     606     15929844 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     607              : 
     608     68122451 :             DO j = 1, ncob
     609     50865120 :                p1 = (j - 1)*ncoa
     610    119533524 :                DO i = 1, ncoa
     611     52738560 :                   p2 = (p1 + i - 1)*ncoc
     612    295448400 :                   DO k = 1, ncoc
     613    191844720 :                      p3 = (p2 + k - 1)*ncod
     614    659265048 :                      DO l = 1, ncod
     615    414681768 :                         p4 = p3 + l
     616    606526488 :                         work2(i, j, k, l, full_perm2(n)) = work(p4, n)
     617              :                      END DO
     618              :                   END DO
     619              :                END DO
     620              :             END DO
     621              :          END DO
     622      1327487 :          IF (use_virial) THEN
     623      1273922 :             DO n = 1, 12
     624              :                tmp_max_virial = 0.0_dp
     625     34681632 :                DO i = 1, mysize
     626              :                   tmp_max_virial = MAX(tmp_max_virial, &
     627     34681632 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     628              :                END DO
     629      1175928 :                tmp_max_virial = tmp_max_virial*max_contraction
     630      1175928 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     631              : 
     632      5078186 :                DO j = 1, ncob
     633      3804264 :                   p1 = (j - 1)*ncoa
     634      8968776 :                   DO i = 1, ncoa
     635      3988584 :                      p2 = (p1 + i - 1)*ncoc
     636     22271544 :                      DO k = 1, ncoc
     637     14478696 :                         p3 = (p2 + k - 1)*ncod
     638     51972984 :                         DO l = 1, ncod
     639     33505704 :                            p4 = p3 + l
     640    148501512 :                            work2_virial(i, j, k, l, full_perm2(n), 1:3) = work_virial(p4, n, 1:3)
     641              :                         END DO
     642              :                      END DO
     643              :                   END DO
     644              :                END DO
     645              :             END DO
     646              :          END IF
     647              :       CASE (3)
     648     17354355 :          DO i = 1, nproducts
     649     52781268 :             A = pgf_product_list(i)%ra
     650     52781268 :             B = pgf_product_list(i)%rb
     651     52781268 :             C = pgf_product_list(i)%rc
     652     52781268 :             D = pgf_product_list(i)%rd
     653     13195317 :             Rho = pgf_product_list(i)%Rho
     654     13195317 :             RhoInv = pgf_product_list(i)%RhoInv
     655     52781268 :             P = pgf_product_list(i)%P
     656     52781268 :             Q = pgf_product_list(i)%Q
     657     52781268 :             W = pgf_product_list(i)%W
     658     52781268 :             AB = pgf_product_list(i)%AB
     659     52781268 :             CD = pgf_product_list(i)%CD
     660     13195317 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
     661              : 
     662              :             CALL build_deriv_data(deriv, A, B, D, C, &
     663              :                                   Zeta_A, Zeta_B, Zeta_D, Zeta_C, &
     664              :                                   ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, &
     665     13195317 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     666     13195317 :             CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize)
     667     52781268 :             DO k = 4, 6
     668    475669875 :                DO j = 1, mysize
     669              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     670              :                                                work_forces(j, k + 3) + &
     671    462474558 :                                                work_forces(j, k + 6))
     672              :                END DO
     673              :             END DO
     674    171539121 :             DO k = 1, 12
     675   1863093549 :                DO j = 1, mysize
     676   1849898232 :                   work(j, k) = work(j, k) + work_forces(j, k)
     677              :                END DO
     678              :             END DO
     679     13195317 :             neris = neris + 12*mysize
     680     17354355 :             IF (use_virial) THEN
     681      1825532 :                CALL real_to_scaled(scoord(1:3), A, cell)
     682      1825532 :                CALL real_to_scaled(scoord(4:6), B, cell)
     683      1825532 :                CALL real_to_scaled(scoord(7:9), D, cell)
     684      1825532 :                CALL real_to_scaled(scoord(10:12), C, cell)
     685     23731916 :                DO k = 1, 12
     686    205909556 :                   DO j = 1, mysize
     687    750616944 :                      DO m = 1, 3
     688    728710560 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     689              :                      END DO
     690              :                   END DO
     691              :                END DO
     692              :             END IF
     693              : 
     694              :          END DO
     695     54067494 :          DO n = 1, 12
     696              :             tmp_max = 0.0_dp
     697    612619524 :             DO i = 1, mysize
     698    612619524 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     699              :             END DO
     700     49908456 :             tmp_max = tmp_max*max_contraction
     701     49908456 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     702              : 
     703    141838302 :             DO i = 1, ncoa
     704     87770808 :                p1 = (i - 1)*ncob
     705    234475776 :                DO j = 1, ncob
     706     96796512 :                   p2 = (p1 + j - 1)*ncod
     707    564673428 :                   DO l = 1, ncod
     708    380106108 :                      p3 = (p2 + l - 1)*ncoc
     709   1039613688 :                      DO k = 1, ncoc
     710    562711068 :                         p4 = p3 + k
     711    942817176 :                         work2(i, j, k, l, full_perm3(n)) = work(p4, n)
     712              :                      END DO
     713              :                   END DO
     714              :                END DO
     715              :             END DO
     716              :          END DO
     717      4159038 :          IF (use_virial) THEN
     718      1829516 :             DO n = 1, 12
     719              :                tmp_max_virial = 0.0_dp
     720     32066520 :                DO i = 1, mysize
     721              :                   tmp_max_virial = MAX(tmp_max_virial, &
     722     32066520 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     723              :                END DO
     724      1688784 :                tmp_max_virial = tmp_max_virial*max_contraction
     725      1688784 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     726              : 
     727      5285348 :                DO i = 1, ncoa
     728      3455832 :                   p1 = (i - 1)*ncob
     729      9372288 :                   DO j = 1, ncob
     730      4227672 :                      p2 = (p1 + j - 1)*ncod
     731     25806840 :                      DO l = 1, ncod
     732     18123336 :                         p3 = (p2 + l - 1)*ncoc
     733     52728744 :                         DO k = 1, ncoc
     734     30377736 :                            p4 = p3 + k
     735    139634280 :                            work2_virial(i, j, k, l, full_perm3(n), 1:3) = work_virial(p4, n, 1:3)
     736              :                         END DO
     737              :                      END DO
     738              :                   END DO
     739              :                END DO
     740              :             END DO
     741              :          END IF
     742              :       CASE (4)
     743      3050715 :          DO i = 1, nproducts
     744      9292336 :             A = pgf_product_list(i)%ra
     745      9292336 :             B = pgf_product_list(i)%rb
     746      9292336 :             C = pgf_product_list(i)%rc
     747      9292336 :             D = pgf_product_list(i)%rd
     748      2323084 :             Rho = pgf_product_list(i)%Rho
     749      2323084 :             RhoInv = pgf_product_list(i)%RhoInv
     750      9292336 :             P = pgf_product_list(i)%P
     751      9292336 :             Q = pgf_product_list(i)%Q
     752      9292336 :             W = pgf_product_list(i)%W
     753      9292336 :             AB = pgf_product_list(i)%AB
     754      9292336 :             CD = pgf_product_list(i)%CD
     755      2323084 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
     756              :             CALL build_deriv_data(deriv, B, A, D, C, &
     757              :                                   Zeta_B, Zeta_A, Zeta_D, Zeta_C, &
     758              :                                   ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, &
     759      2323084 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     760      2323084 :             CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize)
     761      9292336 :             DO k = 4, 6
     762    129888376 :                DO j = 1, mysize
     763              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     764              :                                                work_forces(j, k + 3) + &
     765    127565292 :                                                work_forces(j, k + 6))
     766              :                END DO
     767              :             END DO
     768     30200092 :             DO k = 1, 12
     769    512584252 :                DO j = 1, mysize
     770    510261168 :                   work(j, k) = work(j, k) + work_forces(j, k)
     771              :                END DO
     772              :             END DO
     773      2323084 :             neris = neris + 12*mysize
     774      3050715 :             IF (use_virial) THEN
     775       327543 :                CALL real_to_scaled(scoord(1:3), B, cell)
     776       327543 :                CALL real_to_scaled(scoord(4:6), A, cell)
     777       327543 :                CALL real_to_scaled(scoord(7:9), D, cell)
     778       327543 :                CALL real_to_scaled(scoord(10:12), C, cell)
     779      4258059 :                DO k = 1, 12
     780     59544879 :                   DO j = 1, mysize
     781    225077796 :                      DO m = 1, 3
     782    221147280 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     783              :                      END DO
     784              :                   END DO
     785              :                END DO
     786              :             END IF
     787              : 
     788              :          END DO
     789      9459203 :          DO n = 1, 12
     790              :             tmp_max = 0.0_dp
     791    186700344 :             DO i = 1, mysize
     792    186700344 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     793              :             END DO
     794      8731572 :             tmp_max = tmp_max*max_contraction
     795      8731572 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     796              : 
     797     36990383 :             DO j = 1, ncob
     798     27531180 :                p1 = (j - 1)*ncoa
     799     65337900 :                DO i = 1, ncoa
     800     29075148 :                   p2 = (p1 + i - 1)*ncod
     801    170839980 :                   DO l = 1, ncod
     802    114233652 :                      p3 = (p2 + l - 1)*ncoc
     803    321277572 :                      DO k = 1, ncoc
     804    177968772 :                         p4 = p3 + k
     805    292202424 :                         work2(i, j, k, l, full_perm4(n)) = work(p4, n)
     806              :                      END DO
     807              :                   END DO
     808              :                END DO
     809              :             END DO
     810              :          END DO
     811       727631 :          IF (use_virial) THEN
     812       662467 :             DO n = 1, 12
     813              :                tmp_max_virial = 0.0_dp
     814     14437128 :                DO i = 1, mysize
     815              :                   tmp_max_virial = MAX(tmp_max_virial, &
     816     14437128 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     817              :                END DO
     818       611508 :                tmp_max_virial = tmp_max_virial*max_contraction
     819       611508 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     820              : 
     821      2607583 :                DO j = 1, ncob
     822      1945116 :                   p1 = (j - 1)*ncoa
     823      4649196 :                   DO i = 1, ncoa
     824      2092572 :                      p2 = (p1 + i - 1)*ncod
     825     12389004 :                      DO l = 1, ncod
     826      8351316 :                         p3 = (p2 + l - 1)*ncoc
     827     24269508 :                         DO k = 1, ncoc
     828     13825620 :                            p4 = p3 + k
     829     63653796 :                            work2_virial(i, j, k, l, full_perm4(n), 1:3) = work_virial(p4, n, 1:3)
     830              :                         END DO
     831              :                      END DO
     832              :                   END DO
     833              :                END DO
     834              :             END DO
     835              :          END IF
     836              :       CASE (5)
     837     18046390 :          DO i = 1, nproducts
     838     54398524 :             A = pgf_product_list(i)%ra
     839     54398524 :             B = pgf_product_list(i)%rb
     840     54398524 :             C = pgf_product_list(i)%rc
     841     54398524 :             D = pgf_product_list(i)%rd
     842     13599631 :             Rho = pgf_product_list(i)%Rho
     843     13599631 :             RhoInv = pgf_product_list(i)%RhoInv
     844     54398524 :             P = pgf_product_list(i)%P
     845     54398524 :             Q = pgf_product_list(i)%Q
     846     54398524 :             W = pgf_product_list(i)%W
     847     54398524 :             AB = pgf_product_list(i)%AB
     848     54398524 :             CD = pgf_product_list(i)%CD
     849     13599631 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
     850              :             CALL build_deriv_data(deriv, C, D, A, B, &
     851              :                                   Zeta_C, Zeta_D, Zeta_A, Zeta_B, &
     852              :                                   EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, &
     853     13599631 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     854     13599631 :             CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize)
     855     54398524 :             DO k = 4, 6
     856    530288815 :                DO j = 1, mysize
     857              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     858              :                                                work_forces(j, k + 3) + &
     859    516689184 :                                                work_forces(j, k + 6))
     860              :                END DO
     861              :             END DO
     862    176795203 :             DO k = 1, 12
     863   2080356367 :                DO j = 1, mysize
     864   2066756736 :                   work(j, k) = work(j, k) + work_forces(j, k)
     865              :                END DO
     866              :             END DO
     867     13599631 :             neris = neris + 12*mysize
     868     18046390 :             IF (use_virial) THEN
     869      1995603 :                CALL real_to_scaled(scoord(1:3), C, cell)
     870      1995603 :                CALL real_to_scaled(scoord(4:6), D, cell)
     871      1995603 :                CALL real_to_scaled(scoord(7:9), A, cell)
     872      1995603 :                CALL real_to_scaled(scoord(10:12), B, cell)
     873     25942839 :                DO k = 1, 12
     874    267747675 :                   DO j = 1, mysize
     875    991166580 :                      DO m = 1, 3
     876    967219344 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     877              :                      END DO
     878              :                   END DO
     879              :                END DO
     880              :             END IF
     881              : 
     882              :          END DO
     883     57807867 :          DO n = 1, 12
     884              :             tmp_max = 0.0_dp
     885    678170232 :             DO i = 1, mysize
     886    678170232 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     887              :             END DO
     888     53361108 :             tmp_max = tmp_max*max_contraction
     889     53361108 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     890              : 
     891    130774083 :             DO k = 1, ncoc
     892     72966216 :                p1 = (k - 1)*ncod
     893    203930556 :                DO l = 1, ncod
     894     77603232 :                   p2 = (p1 + l - 1)*ncoa
     895    447100044 :                   DO i = 1, ncoa
     896    296530596 :                      p3 = (p2 + i - 1)*ncob
     897    998942952 :                      DO j = 1, ncob
     898    624809124 :                         p4 = p3 + j
     899    921339720 :                         work2(i, j, k, l, full_perm5(n)) = work(p4, n)
     900              :                      END DO
     901              :                   END DO
     902              :                END DO
     903              :             END DO
     904              :          END DO
     905      4446759 :          IF (use_virial) THEN
     906      2252003 :             DO n = 1, 12
     907              :                tmp_max_virial = 0.0_dp
     908     44899908 :                DO i = 1, mysize
     909              :                   tmp_max_virial = MAX(tmp_max_virial, &
     910     44899908 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     911              :                END DO
     912      2078772 :                tmp_max_virial = tmp_max_virial*max_contraction
     913      2078772 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     914              : 
     915      5747327 :                DO k = 1, ncoc
     916      3495324 :                   p1 = (k - 1)*ncod
     917      9490044 :                   DO l = 1, ncod
     918      3915948 :                      p2 = (p1 + l - 1)*ncoa
     919     22978104 :                      DO i = 1, ncoa
     920     15566832 :                         p3 = (p2 + i - 1)*ncob
     921     62303916 :                         DO j = 1, ncob
     922     42821136 :                            p4 = p3 + j
     923    186851376 :                            work2_virial(i, j, k, l, full_perm5(n), 1:3) = work_virial(p4, n, 1:3)
     924              :                         END DO
     925              :                      END DO
     926              :                   END DO
     927              :                END DO
     928              :             END DO
     929              :          END IF
     930              :       CASE (6)
     931      4667905 :          DO i = 1, nproducts
     932     14050012 :             A = pgf_product_list(i)%ra
     933     14050012 :             B = pgf_product_list(i)%rb
     934     14050012 :             C = pgf_product_list(i)%rc
     935     14050012 :             D = pgf_product_list(i)%rd
     936      3512503 :             Rho = pgf_product_list(i)%Rho
     937      3512503 :             RhoInv = pgf_product_list(i)%RhoInv
     938     14050012 :             P = pgf_product_list(i)%P
     939     14050012 :             Q = pgf_product_list(i)%Q
     940     14050012 :             W = pgf_product_list(i)%W
     941     14050012 :             AB = pgf_product_list(i)%AB
     942     14050012 :             CD = pgf_product_list(i)%CD
     943      3512503 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
     944              : 
     945              :             CALL build_deriv_data(deriv, C, D, B, A, &
     946              :                                   Zeta_C, Zeta_D, Zeta_B, Zeta_A, &
     947              :                                   EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, &
     948      3512503 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     949      3512503 :             CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize)
     950     14050012 :             DO k = 4, 6
     951    136657552 :                DO j = 1, mysize
     952              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     953              :                                                work_forces(j, k + 3) + &
     954    133145049 :                                                work_forces(j, k + 6))
     955              :                END DO
     956              :             END DO
     957     45662539 :             DO k = 1, 12
     958    536092699 :                DO j = 1, mysize
     959    532580196 :                   work(j, k) = work(j, k) + work_forces(j, k)
     960              :                END DO
     961              :             END DO
     962      3512503 :             neris = neris + 12*mysize
     963      4667905 :             IF (use_virial) THEN
     964       381872 :                CALL real_to_scaled(scoord(1:3), C, cell)
     965       381872 :                CALL real_to_scaled(scoord(4:6), D, cell)
     966       381872 :                CALL real_to_scaled(scoord(7:9), B, cell)
     967       381872 :                CALL real_to_scaled(scoord(10:12), A, cell)
     968      4964336 :                DO k = 1, 12
     969     51070436 :                   DO j = 1, mysize
     970    189006864 :                      DO m = 1, 3
     971    184424400 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     972              :                      END DO
     973              :                   END DO
     974              :                END DO
     975              :             END IF
     976              : 
     977              :          END DO
     978     15020226 :          DO n = 1, 12
     979              :             tmp_max = 0.0_dp
     980    195212556 :             DO i = 1, mysize
     981    195212556 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     982              :             END DO
     983     13864824 :             tmp_max = tmp_max*max_contraction
     984     13864824 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     985              : 
     986     32866194 :             DO k = 1, ncoc
     987     17845968 :                p1 = (k - 1)*ncod
     988     52452456 :                DO l = 1, ncod
     989     20741664 :                   p2 = (p1 + l - 1)*ncob
     990    131403588 :                   DO j = 1, ncob
     991     92815956 :                      p3 = (p2 + j - 1)*ncoa
     992    294905352 :                      DO i = 1, ncoa
     993    181347732 :                         p4 = p3 + i
     994    274163688 :                         work2(i, j, k, l, full_perm6(n)) = work(p4, n)
     995              :                      END DO
     996              :                   END DO
     997              :                END DO
     998              :             END DO
     999              :          END DO
    1000      1155402 :          IF (use_virial) THEN
    1001       836784 :             DO n = 1, 12
    1002              :                tmp_max_virial = 0.0_dp
    1003     16323840 :                DO i = 1, mysize
    1004              :                   tmp_max_virial = MAX(tmp_max_virial, &
    1005     16323840 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1006              :                END DO
    1007       772416 :                tmp_max_virial = tmp_max_virial*max_contraction
    1008       772416 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1009              : 
    1010      1965552 :                DO k = 1, ncoc
    1011      1128768 :                   p1 = (k - 1)*ncod
    1012      3324864 :                   DO l = 1, ncod
    1013      1423680 :                      p2 = (p1 + l - 1)*ncob
    1014      9551424 :                      DO j = 1, ncob
    1015      6998976 :                         p3 = (p2 + j - 1)*ncoa
    1016     23974080 :                         DO i = 1, ncoa
    1017     15551424 :                            p4 = p3 + i
    1018     69204672 :                            work2_virial(i, j, k, l, full_perm6(n), 1:3) = work_virial(p4, n, 1:3)
    1019              :                         END DO
    1020              :                      END DO
    1021              :                   END DO
    1022              :                END DO
    1023              :             END DO
    1024              :          END IF
    1025              :       CASE (7)
    1026      2744720 :          DO i = 1, nproducts
    1027      8434988 :             A = pgf_product_list(i)%ra
    1028      8434988 :             B = pgf_product_list(i)%rb
    1029      8434988 :             C = pgf_product_list(i)%rc
    1030      8434988 :             D = pgf_product_list(i)%rd
    1031      2108747 :             Rho = pgf_product_list(i)%Rho
    1032      2108747 :             RhoInv = pgf_product_list(i)%RhoInv
    1033      8434988 :             P = pgf_product_list(i)%P
    1034      8434988 :             Q = pgf_product_list(i)%Q
    1035      8434988 :             W = pgf_product_list(i)%W
    1036      8434988 :             AB = pgf_product_list(i)%AB
    1037      8434988 :             CD = pgf_product_list(i)%CD
    1038      2108747 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1039              : 
    1040              :             CALL build_deriv_data(deriv, D, C, A, B, &
    1041              :                                   Zeta_D, Zeta_C, Zeta_A, Zeta_B, &
    1042              :                                   EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, &
    1043      2108747 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1044      2108747 :             CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize)
    1045      8434988 :             DO k = 4, 6
    1046    206547515 :                DO j = 1, mysize
    1047              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1048              :                                                work_forces(j, k + 3) + &
    1049    204438768 :                                                work_forces(j, k + 6))
    1050              :                END DO
    1051              :             END DO
    1052     27413711 :             DO k = 1, 12
    1053    819863819 :                DO j = 1, mysize
    1054    817755072 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1055              :                END DO
    1056              :             END DO
    1057      2108747 :             neris = neris + 12*mysize
    1058      2744720 :             IF (use_virial) THEN
    1059       325942 :                CALL real_to_scaled(scoord(1:3), D, cell)
    1060       325942 :                CALL real_to_scaled(scoord(4:6), C, cell)
    1061       325942 :                CALL real_to_scaled(scoord(7:9), A, cell)
    1062       325942 :                CALL real_to_scaled(scoord(10:12), B, cell)
    1063      4237246 :                DO k = 1, 12
    1064    119756746 :                   DO j = 1, mysize
    1065    465989304 :                      DO m = 1, 3
    1066    462078000 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
    1067              :                      END DO
    1068              :                   END DO
    1069              :                END DO
    1070              :             END IF
    1071              : 
    1072              :          END DO
    1073      8267649 :          DO n = 1, 12
    1074              :             tmp_max = 0.0_dp
    1075    273131856 :             DO i = 1, mysize
    1076    273131856 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1077              :             END DO
    1078      7631676 :             tmp_max = tmp_max*max_contraction
    1079      7631676 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1080              : 
    1081     31880121 :             DO l = 1, ncod
    1082     23612472 :                p1 = (l - 1)*ncoc
    1083     55298124 :                DO k = 1, ncoc
    1084     24053976 :                   p2 = (p1 + k - 1)*ncoa
    1085    145878300 :                   DO i = 1, ncoa
    1086     98211852 :                      p3 = (p2 + i - 1)*ncob
    1087    387766008 :                      DO j = 1, ncob
    1088    265500180 :                         p4 = p3 + j
    1089    363712032 :                         work2(i, j, k, l, full_perm7(n)) = work(p4, n)
    1090              :                      END DO
    1091              :                   END DO
    1092              :                END DO
    1093              :             END DO
    1094              :          END DO
    1095       635973 :          IF (use_virial) THEN
    1096       640536 :             DO n = 1, 12
    1097              :                tmp_max_virial = 0.0_dp
    1098     21929472 :                DO i = 1, mysize
    1099              :                   tmp_max_virial = MAX(tmp_max_virial, &
    1100     21929472 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1101              :                END DO
    1102       591264 :                tmp_max_virial = tmp_max_virial*max_contraction
    1103       591264 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1104              : 
    1105      2470920 :                DO l = 1, ncod
    1106      1830384 :                   p1 = (l - 1)*ncoc
    1107      4288896 :                   DO k = 1, ncoc
    1108      1867248 :                      p2 = (p1 + k - 1)*ncoa
    1109     10867968 :                      DO i = 1, ncoa
    1110      7170336 :                         p3 = (p2 + i - 1)*ncob
    1111     30375792 :                         DO j = 1, ncob
    1112     21338208 :                            p4 = p3 + j
    1113     92523168 :                            work2_virial(i, j, k, l, full_perm7(n), 1:3) = work_virial(p4, n, 1:3)
    1114              :                         END DO
    1115              :                      END DO
    1116              :                   END DO
    1117              :                END DO
    1118              :             END DO
    1119              :          END IF
    1120              :       CASE (8)
    1121       373187 :          DO i = 1, nproducts
    1122      1057412 :             A = pgf_product_list(i)%ra
    1123      1057412 :             B = pgf_product_list(i)%rb
    1124      1057412 :             C = pgf_product_list(i)%rc
    1125      1057412 :             D = pgf_product_list(i)%rd
    1126       264353 :             Rho = pgf_product_list(i)%Rho
    1127       264353 :             RhoInv = pgf_product_list(i)%RhoInv
    1128      1057412 :             P = pgf_product_list(i)%P
    1129      1057412 :             Q = pgf_product_list(i)%Q
    1130      1057412 :             W = pgf_product_list(i)%W
    1131      1057412 :             AB = pgf_product_list(i)%AB
    1132      1057412 :             CD = pgf_product_list(i)%CD
    1133       264353 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1134              : 
    1135              :             CALL build_deriv_data(deriv, D, C, B, A, &
    1136              :                                   Zeta_D, Zeta_C, Zeta_B, Zeta_A, &
    1137              :                                   EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, &
    1138       264353 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1139       264353 :             CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize)
    1140      1057412 :             DO k = 4, 6
    1141     33923648 :                DO j = 1, mysize
    1142              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1143              :                                                work_forces(j, k + 3) + &
    1144     33659295 :                                                work_forces(j, k + 6))
    1145              :                END DO
    1146              :             END DO
    1147      3436589 :             DO k = 1, 12
    1148    134901533 :                DO j = 1, mysize
    1149    134637180 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1150              :                END DO
    1151              :             END DO
    1152       264353 :             neris = neris + 12*mysize
    1153       373187 :             IF (use_virial) THEN
    1154        22456 :                CALL real_to_scaled(scoord(1:3), D, cell)
    1155        22456 :                CALL real_to_scaled(scoord(4:6), C, cell)
    1156        22456 :                CALL real_to_scaled(scoord(7:9), B, cell)
    1157        22456 :                CALL real_to_scaled(scoord(10:12), A, cell)
    1158       291928 :                DO k = 1, 12
    1159     11660440 :                   DO j = 1, mysize
    1160     45743520 :                      DO m = 1, 3
    1161     45474048 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
    1162              :                      END DO
    1163              :                   END DO
    1164              :                END DO
    1165              :             END IF
    1166              : 
    1167              :          END DO
    1168      1414842 :          DO n = 1, 12
    1169              :             tmp_max = 0.0_dp
    1170     56524248 :             DO i = 1, mysize
    1171     56524248 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1172              :             END DO
    1173      1306008 :             tmp_max = tmp_max*max_contraction
    1174      1306008 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1175              : 
    1176      5757306 :             DO l = 1, ncod
    1177      4342464 :                p1 = (l - 1)*ncoc
    1178      9990936 :                DO k = 1, ncoc
    1179      4342464 :                   p2 = (p1 + k - 1)*ncob
    1180     34739712 :                   DO j = 1, ncob
    1181     26054784 :                      p3 = (p2 + j - 1)*ncoa
    1182     85615488 :                      DO i = 1, ncoa
    1183     55218240 :                         p4 = p3 + i
    1184     81273024 :                         work2(i, j, k, l, full_perm8(n)) = work(p4, n)
    1185              :                      END DO
    1186              :                   END DO
    1187              :                END DO
    1188              :             END DO
    1189              :          END DO
    1190     34417068 :          IF (use_virial) THEN
    1191       119808 :             DO n = 1, 12
    1192              :                tmp_max_virial = 0.0_dp
    1193      4976640 :                DO i = 1, mysize
    1194              :                   tmp_max_virial = MAX(tmp_max_virial, &
    1195      4976640 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1196              :                END DO
    1197       110592 :                tmp_max_virial = tmp_max_virial*max_contraction
    1198       110592 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1199              : 
    1200       488448 :                DO l = 1, ncod
    1201       368640 :                   p1 = (l - 1)*ncoc
    1202       847872 :                   DO k = 1, ncoc
    1203       368640 :                      p2 = (p1 + k - 1)*ncob
    1204      2949120 :                      DO j = 1, ncob
    1205      2211840 :                         p3 = (p2 + j - 1)*ncoa
    1206      7446528 :                         DO i = 1, ncoa
    1207      4866048 :                            p4 = p3 + i
    1208     21676032 :                            work2_virial(i, j, k, l, full_perm8(n), 1:3) = work_virial(p4, n, 1:3)
    1209              :                         END DO
    1210              :                      END DO
    1211              :                   END DO
    1212              :                END DO
    1213              :             END DO
    1214              :          END IF
    1215              :       END SELECT
    1216              : 
    1217     34308234 :       IF (.NOT. use_virial) THEN
    1218     33242950 :          IF (tmp_max_all < eps_schwarz) RETURN
    1219              :       END IF
    1220              : 
    1221     28626339 :       IF (tmp_max_all >= eps_schwarz) THEN
    1222    369983718 :          DO i = 1, 12
    1223  20256805776 :             primitives_tmp(:, :, :, :) = 0.0_dp
    1224              :             CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
    1225              :                           n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2(1, 1, 1, 1, i), &
    1226              :                           sphi_a, &
    1227              :                           sphi_b, &
    1228              :                           sphi_c, &
    1229              :                           sphi_d, &
    1230              :                           primitives_tmp(1, 1, 1, 1), &
    1231    341523432 :                           buffer1, buffer2)
    1232              : 
    1233              :             primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1234              :                        s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1235              :                        s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1236              :                        s_offset_d + 1:s_offset_d + nsod*nl_d, i) = &
    1237              :                primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1238              :                           s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1239              :                           s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1240  20285266062 :                           s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
    1241              :          END DO
    1242              :       END IF
    1243              : 
    1244     28626339 :       IF (use_virial .AND. tmp_max_all_virial >= eps_schwarz) THEN
    1245     11813685 :          DO i = 1, 12
    1246     44528505 :             DO m = 1, 3
    1247   4937303412 :                primitives_tmp_virial(:, :, :, :) = 0.0_dp
    1248              :                CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
    1249              :                              n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2_virial(1, 1, 1, 1, i, m), &
    1250              :                              sphi_a, &
    1251              :                              sphi_b, &
    1252              :                              sphi_c, &
    1253              :                              sphi_d, &
    1254              :                              primitives_tmp_virial(1, 1, 1, 1), &
    1255     32714820 :                              buffer1, buffer2)
    1256              : 
    1257              :                primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1258              :                                  s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1259              :                                  s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1260              :                                  s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) = &
    1261              :                   primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1262              :                                     s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1263              :                                     s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1264   4948208352 :                                     s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) + primitives_tmp_virial(:, :, :, :)
    1265              :             END DO
    1266              :          END DO
    1267              :       END IF
    1268              : 
    1269              :    END SUBROUTINE evaluate_deriv_eri
    1270              : 
    1271              : ! **************************************************************************************************
    1272              : !> \brief Evaluates max-abs values of  electron repulsion integrals for a primitive quartet
    1273              : !> \param libint ...
    1274              : !> \param A ...
    1275              : !> \param B ...
    1276              : !> \param C ...
    1277              : !> \param D ...
    1278              : !> \param Zeta_A ...
    1279              : !> \param Zeta_B ...
    1280              : !> \param Zeta_C ...
    1281              : !> \param Zeta_D ...
    1282              : !> \param n_a ...
    1283              : !> \param n_b ...
    1284              : !> \param n_c ...
    1285              : !> \param n_d ...
    1286              : !> \param max_val ...
    1287              : !> \param potential_parameter ...
    1288              : !> \param R1 ...
    1289              : !> \param R2 ...
    1290              : !> \param p_work ...
    1291              : !> \par History
    1292              : !>      03.2007 created [Manuel Guidon]
    1293              : !>      08.2007 refactured permutation part [Manuel Guidon]
    1294              : !> \author Manuel Guidon
    1295              : ! **************************************************************************************************
    1296     87276120 :    SUBROUTINE evaluate_eri_screen(libint, A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, &
    1297              :                                   n_a, n_b, n_c, n_d, &
    1298              :                                   max_val, potential_parameter, R1, R2, &
    1299              :                                   p_work)
    1300              : 
    1301              :       TYPE(cp_libint_t)                                  :: libint
    1302              :       REAL(dp), INTENT(IN)                               :: A(3), B(3), C(3), D(3), Zeta_A, Zeta_B, &
    1303              :                                                             Zeta_C, Zeta_D
    1304              :       INTEGER, INTENT(IN)                                :: n_a, n_b, n_c, n_d
    1305              :       REAL(dp), INTENT(INOUT)                            :: max_val
    1306              :       TYPE(hfx_potential_type)                           :: potential_parameter
    1307              :       REAL(dp)                                           :: R1, R2
    1308              :       REAL(dp), DIMENSION(:), POINTER                    :: p_work
    1309              : 
    1310              :       INTEGER                                            :: a_mysize(1), i, m_max, mysize, perm_case
    1311              : 
    1312     87276120 :       m_max = n_a + n_b + n_c + n_d
    1313     87276120 :       mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
    1314    174552240 :       a_mysize = mysize
    1315              : 
    1316     87276120 :       IF (m_max /= 0) THEN
    1317     54836334 :          perm_case = 1
    1318     54836334 :          IF (n_a < n_b) THEN
    1319     22260198 :             perm_case = perm_case + 1
    1320              :          END IF
    1321     54836334 :          IF (n_c < n_d) THEN
    1322     22260198 :             perm_case = perm_case + 2
    1323              :          END IF
    1324     54836334 :          IF (n_a + n_b > n_c + n_d) THEN
    1325            0 :             perm_case = perm_case + 4
    1326              :          END IF
    1327              : 
    1328     32576136 :          SELECT CASE (perm_case)
    1329              :          CASE (1)
    1330              :             CALL build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, &
    1331     32576136 :                                            potential_parameter, libint, R1, R2)
    1332              : 
    1333     32576136 :             CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1334   2796479112 :             DO i = 1, mysize
    1335   2796479112 :                max_val = MAX(max_val, ABS(p_work(i)))
    1336              :             END DO
    1337              :          CASE (2)
    1338              :             CALL build_quartet_data_screen(B, A, C, D, Zeta_B, Zeta_A, Zeta_C, Zeta_D, m_max, &
    1339            0 :                                            potential_parameter, libint, R1, R2)
    1340              : 
    1341            0 :             CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
    1342            0 :             DO i = 1, mysize
    1343            0 :                max_val = MAX(max_val, ABS(p_work(i)))
    1344              :             END DO
    1345              :          CASE (3)
    1346              :             CALL build_quartet_data_screen(A, B, D, C, Zeta_A, Zeta_B, Zeta_D, Zeta_C, m_max, &
    1347            0 :                                            potential_parameter, libint, R1, R2)
    1348              : 
    1349            0 :             CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
    1350            0 :             DO i = 1, mysize
    1351            0 :                max_val = MAX(max_val, ABS(p_work(i)))
    1352              :             END DO
    1353              :          CASE (4)
    1354              :             CALL build_quartet_data_screen(B, A, D, C, Zeta_B, Zeta_A, Zeta_D, Zeta_C, m_max, &
    1355     22260198 :                                            potential_parameter, libint, R1, R2)
    1356              : 
    1357     22260198 :             CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1358   1005563070 :             DO i = 1, mysize
    1359   1005563070 :                max_val = MAX(max_val, ABS(p_work(i)))
    1360              :             END DO
    1361              :          CASE (5)
    1362              :             CALL build_quartet_data_screen(C, D, A, B, Zeta_C, Zeta_D, Zeta_A, Zeta_B, m_max, &
    1363            0 :                                            potential_parameter, libint, R1, R2)
    1364              : 
    1365            0 :             CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
    1366            0 :             DO i = 1, mysize
    1367            0 :                max_val = MAX(max_val, ABS(p_work(i)))
    1368              :             END DO
    1369              :          CASE (6)
    1370              :             CALL build_quartet_data_screen(C, D, B, A, Zeta_C, Zeta_D, Zeta_B, Zeta_A, m_max, &
    1371            0 :                                            potential_parameter, libint, R1, R2)
    1372              : 
    1373            0 :             CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
    1374            0 :             DO i = 1, mysize
    1375            0 :                max_val = MAX(max_val, ABS(p_work(i)))
    1376              :             END DO
    1377              :          CASE (7)
    1378              :             CALL build_quartet_data_screen(D, C, A, B, Zeta_D, Zeta_C, Zeta_A, Zeta_B, m_max, &
    1379            0 :                                            potential_parameter, libint, R1, R2)
    1380              : 
    1381            0 :             CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
    1382            0 :             DO i = 1, mysize
    1383            0 :                max_val = MAX(max_val, ABS(p_work(i)))
    1384              :             END DO
    1385              :          CASE (8)
    1386              :             CALL build_quartet_data_screen(D, C, B, A, Zeta_D, Zeta_C, Zeta_B, Zeta_A, m_max, &
    1387            0 :                                            potential_parameter, libint, R1, R2)
    1388              : 
    1389            0 :             CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
    1390     54836334 :             DO i = 1, mysize
    1391            0 :                max_val = MAX(max_val, ABS(p_work(i)))
    1392              :             END DO
    1393              :          END SELECT
    1394              :       ELSE
    1395              :          CALL build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, &
    1396     32439786 :                                         potential_parameter, libint, R1, R2)
    1397     32439786 :          max_val = ABS(get_ssss_f_val(libint))
    1398              :       END IF
    1399              : 
    1400     87276120 :    END SUBROUTINE evaluate_eri_screen
    1401              : 
    1402              : ! **************************************************************************************************
    1403              : !> \brief Evaluate electron repulsion integrals for a primitive quartet
    1404              : !> \param libint ...
    1405              : !> \param nproducts ...
    1406              : !> \param pgf_product_list ...
    1407              : !> \param n_a ...
    1408              : !> \param n_b ...
    1409              : !> \param n_c ...
    1410              : !> \param n_d ...
    1411              : !> \param ncoa ...
    1412              : !> \param ncob ...
    1413              : !> \param ncoc ...
    1414              : !> \param ncod ...
    1415              : !> \param nsgfa ...
    1416              : !> \param nsgfb ...
    1417              : !> \param nsgfc ...
    1418              : !> \param nsgfd ...
    1419              : !> \param primitives ...
    1420              : !> \param max_contraction ...
    1421              : !> \param tmp_max ...
    1422              : !> \param eps_schwarz ...
    1423              : !> \param neris ...
    1424              : !> \param ZetaInv ...
    1425              : !> \param EtaInv ...
    1426              : !> \param s_offset_a ...
    1427              : !> \param s_offset_b ...
    1428              : !> \param s_offset_c ...
    1429              : !> \param s_offset_d ...
    1430              : !> \param nl_a ...
    1431              : !> \param nl_b ...
    1432              : !> \param nl_c ...
    1433              : !> \param nl_d ...
    1434              : !> \param nsoa ...
    1435              : !> \param nsob ...
    1436              : !> \param nsoc ...
    1437              : !> \param nsod ...
    1438              : !> \param sphi_a ...
    1439              : !> \param sphi_b ...
    1440              : !> \param sphi_c ...
    1441              : !> \param sphi_d ...
    1442              : !> \param work ...
    1443              : !> \param work2 ...
    1444              : !> \param buffer1 ...
    1445              : !> \param buffer2 ...
    1446              : !> \param primitives_tmp ...
    1447              : !> \param p_work ...
    1448              : !> \par History
    1449              : !>      11.2006 created [Manuel Guidon]
    1450              : !>      08.2007 refactured permutation part [Manuel Guidon]
    1451              : !> \author Manuel Guidon
    1452              : ! **************************************************************************************************
    1453    176174052 :    SUBROUTINE evaluate_eri(libint, nproducts, pgf_product_list, &
    1454              :                            n_a, n_b, n_c, n_d, &
    1455              :                            ncoa, ncob, ncoc, ncod, &
    1456              :                            nsgfa, nsgfb, nsgfc, nsgfd, &
    1457    176174052 :                            primitives, max_contraction, tmp_max, &
    1458              :                            eps_schwarz, neris, &
    1459              :                            ZetaInv, EtaInv, &
    1460              :                            s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
    1461              :                            nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, &
    1462    176174052 :                            sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
    1463    176174052 :                            primitives_tmp, p_work)
    1464              : 
    1465              :       TYPE(cp_libint_t)                                  :: libint
    1466              :       INTEGER, INTENT(IN)                                :: nproducts
    1467              :       TYPE(hfx_pgf_product_list), DIMENSION(nproducts)   :: pgf_product_list
    1468              :       INTEGER, INTENT(IN)                                :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, &
    1469              :                                                             ncod, nsgfa, nsgfb, nsgfc, nsgfd
    1470              :       REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd)    :: primitives
    1471              :       REAL(dp)                                           :: max_contraction, tmp_max, eps_schwarz
    1472              :       INTEGER(int_8)                                     :: neris
    1473              :       REAL(dp), INTENT(IN)                               :: ZetaInv, EtaInv
    1474              :       INTEGER                                            :: s_offset_a, s_offset_b, s_offset_c, &
    1475              :                                                             s_offset_d, nl_a, nl_b, nl_c, nl_d, &
    1476              :                                                             nsoa, nsob, nsoc, nsod
    1477              :       REAL(dp), DIMENSION(ncoa, nsoa*nl_a), INTENT(IN)   :: sphi_a
    1478              :       REAL(dp), DIMENSION(ncob, nsob*nl_b), INTENT(IN)   :: sphi_b
    1479              :       REAL(dp), DIMENSION(ncoc, nsoc*nl_c), INTENT(IN)   :: sphi_c
    1480              :       REAL(dp), DIMENSION(ncod, nsod*nl_d), INTENT(IN)   :: sphi_d
    1481              :       REAL(dp)                                           :: work(ncoa*ncob*ncoc*ncod), &
    1482              :                                                             work2(ncoa, ncob, ncoc, ncod)
    1483              :       REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod)           :: buffer1, buffer2
    1484              :       REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*&
    1485              :          nl_c, nsod*nl_d)                                :: primitives_tmp
    1486              :       REAL(dp), DIMENSION(:), POINTER                    :: p_work
    1487              : 
    1488              :       INTEGER                                            :: a_mysize(1), i, j, k, l, m_max, mysize, &
    1489              :                                                             p1, p2, p3, p4, perm_case
    1490              :       REAL(dp)                                           :: A(3), AB(3), B(3), C(3), CD(3), D(3), &
    1491              :                                                             P(3), Q(3), Rho, RhoInv, W(3), &
    1492              :                                                             ZetapEtaInv
    1493              :       REAL(KIND=dp), DIMENSION(prim_data_f_size)         :: F
    1494              : 
    1495    176174052 :       m_max = n_a + n_b + n_c + n_d
    1496    176174052 :       mysize = ncoa*ncob*ncoc*ncod
    1497    352348104 :       a_mysize = mysize
    1498              : 
    1499   3614769530 :       work = 0.0_dp
    1500    176174052 :       IF (m_max /= 0) THEN
    1501    139385347 :          perm_case = 1
    1502    139385347 :          IF (n_a < n_b) THEN
    1503     33885765 :             perm_case = perm_case + 1
    1504              :          END IF
    1505    139385347 :          IF (n_c < n_d) THEN
    1506     40561352 :             perm_case = perm_case + 2
    1507              :          END IF
    1508    139385347 :          IF (n_a + n_b > n_c + n_d) THEN
    1509     49222448 :             perm_case = perm_case + 4
    1510              :          END IF
    1511              :          SELECT CASE (perm_case)
    1512              :          CASE (1)
    1513    200172830 :             DO i = 1, nproducts
    1514    611773916 :                A = pgf_product_list(i)%ra
    1515    611773916 :                B = pgf_product_list(i)%rb
    1516    611773916 :                C = pgf_product_list(i)%rc
    1517    611773916 :                D = pgf_product_list(i)%rd
    1518    152943479 :                Rho = pgf_product_list(i)%Rho
    1519    152943479 :                RhoInv = pgf_product_list(i)%RhoInv
    1520    611773916 :                P = pgf_product_list(i)%P
    1521    611773916 :                Q = pgf_product_list(i)%Q
    1522    611773916 :                W = pgf_product_list(i)%W
    1523    611773916 :                AB = pgf_product_list(i)%AB
    1524    611773916 :                CD = pgf_product_list(i)%CD
    1525    152943479 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1526    607892027 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1527              : 
    1528              :                CALL build_quartet_data(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1529    152943479 :                                        P, Q, W, m_max, F)
    1530    152943479 :                CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1531   2797974197 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1532    200172830 :                neris = neris + mysize
    1533              :             END DO
    1534   1134116021 :             DO i = 1, mysize
    1535   1134116021 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1536              :             END DO
    1537     47229351 :             tmp_max = tmp_max*max_contraction
    1538     47229351 :             IF (tmp_max < eps_schwarz) THEN
    1539     44770173 :                RETURN
    1540              :             END IF
    1541              : 
    1542    108035362 :             DO i = 1, ncoa
    1543     71096050 :                p1 = (i - 1)*ncob
    1544    198688082 :                DO j = 1, ncob
    1545     90652720 :                   p2 = (p1 + j - 1)*ncoc
    1546    525995141 :                   DO k = 1, ncoc
    1547    364246371 :                      p3 = (p2 + k - 1)*ncod
    1548   1383048762 :                      DO l = 1, ncod
    1549    928149671 :                         p4 = p3 + l
    1550   1292396042 :                         work2(i, j, k, l) = work(p4)
    1551              :                      END DO
    1552              :                   END DO
    1553              :                END DO
    1554              :             END DO
    1555              :          CASE (2)
    1556     35159307 :             DO i = 1, nproducts
    1557    100549340 :                A = pgf_product_list(i)%ra
    1558    100549340 :                B = pgf_product_list(i)%rb
    1559    100549340 :                C = pgf_product_list(i)%rc
    1560    100549340 :                D = pgf_product_list(i)%rd
    1561     25137335 :                Rho = pgf_product_list(i)%Rho
    1562     25137335 :                RhoInv = pgf_product_list(i)%RhoInv
    1563    100549340 :                P = pgf_product_list(i)%P
    1564    100549340 :                Q = pgf_product_list(i)%Q
    1565    100549340 :                W = pgf_product_list(i)%W
    1566    100549340 :                AB = pgf_product_list(i)%AB
    1567    100549340 :                CD = pgf_product_list(i)%CD
    1568     25137335 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1569    119744986 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1570              : 
    1571              :                CALL build_quartet_data(libint, B, A, C, D, &
    1572              :                                        ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1573     25137335 :                                        P, Q, W, m_max, F)
    1574              : 
    1575     25137335 :                CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
    1576    783143829 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1577     35159307 :                neris = neris + mysize
    1578              :             END DO
    1579    421690711 :             DO i = 1, mysize
    1580    421690711 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1581              :             END DO
    1582     10021972 :             tmp_max = tmp_max*max_contraction
    1583     10021972 :             IF (tmp_max < eps_schwarz) THEN
    1584              :                RETURN
    1585              :             END IF
    1586              : 
    1587     32154085 :             DO j = 1, ncob
    1588     25067570 :                p1 = (j - 1)*ncoa
    1589     60209581 :                DO i = 1, ncoa
    1590     28055496 :                   p2 = (p1 + i - 1)*ncoc
    1591    173270724 :                   DO k = 1, ncoc
    1592    120147658 :                      p3 = (p2 + k - 1)*ncod
    1593    473268832 :                      DO l = 1, ncod
    1594    325065678 :                         p4 = p3 + l
    1595    445213336 :                         work2(i, j, k, l) = work(p4)
    1596              :                      END DO
    1597              :                   END DO
    1598              :                END DO
    1599              :             END DO
    1600              :          CASE (3)
    1601     89745804 :             DO i = 1, nproducts
    1602    256072944 :                A = pgf_product_list(i)%ra
    1603    256072944 :                B = pgf_product_list(i)%rb
    1604    256072944 :                C = pgf_product_list(i)%rc
    1605    256072944 :                D = pgf_product_list(i)%rd
    1606     64018236 :                Rho = pgf_product_list(i)%Rho
    1607     64018236 :                RhoInv = pgf_product_list(i)%RhoInv
    1608    256072944 :                P = pgf_product_list(i)%P
    1609    256072944 :                Q = pgf_product_list(i)%Q
    1610    256072944 :                W = pgf_product_list(i)%W
    1611    256072944 :                AB = pgf_product_list(i)%AB
    1612    256072944 :                CD = pgf_product_list(i)%CD
    1613     64018236 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1614    242375185 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1615              : 
    1616              :                CALL build_quartet_data(libint, A, B, D, C, &
    1617              :                                        ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1618     64018236 :                                        P, Q, W, m_max, F)
    1619              : 
    1620     64018236 :                CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
    1621    940288743 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1622     89745804 :                neris = neris + mysize
    1623              :             END DO
    1624    461145875 :             DO i = 1, mysize
    1625    461145875 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1626              :             END DO
    1627     25727568 :             tmp_max = tmp_max*max_contraction
    1628     25727568 :             IF (tmp_max < eps_schwarz) THEN
    1629              :                RETURN
    1630              :             END IF
    1631              : 
    1632     54285102 :             DO i = 1, ncoa
    1633     36146485 :                p1 = (i - 1)*ncob
    1634     98089651 :                DO j = 1, ncob
    1635     43804549 :                   p2 = (p1 + j - 1)*ncod
    1636    282444246 :                   DO l = 1, ncod
    1637    202493212 :                      p3 = (p2 + l - 1)*ncoc
    1638    591259447 :                      DO k = 1, ncoc
    1639    344961686 :                         p4 = p3 + k
    1640    547454898 :                         work2(i, j, k, l) = work(p4)
    1641              :                      END DO
    1642              :                   END DO
    1643              :                END DO
    1644              :             END DO
    1645              :          CASE (4)
    1646     21755577 :             DO i = 1, nproducts
    1647     58286276 :                A = pgf_product_list(i)%ra
    1648     58286276 :                B = pgf_product_list(i)%rb
    1649     58286276 :                C = pgf_product_list(i)%rc
    1650     58286276 :                D = pgf_product_list(i)%rd
    1651     14571569 :                Rho = pgf_product_list(i)%Rho
    1652     14571569 :                RhoInv = pgf_product_list(i)%RhoInv
    1653     58286276 :                P = pgf_product_list(i)%P
    1654     58286276 :                Q = pgf_product_list(i)%Q
    1655     58286276 :                W = pgf_product_list(i)%W
    1656     58286276 :                AB = pgf_product_list(i)%AB
    1657     58286276 :                CD = pgf_product_list(i)%CD
    1658     14571569 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1659     66374819 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1660              : 
    1661              :                CALL build_quartet_data(libint, B, A, D, C, &
    1662              :                                        ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1663     14571569 :                                        P, Q, W, m_max, F)
    1664              : 
    1665     14571569 :                CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1666    389605651 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1667     21755577 :                neris = neris + mysize
    1668              :             END DO
    1669    246883327 :             DO i = 1, mysize
    1670    246883327 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1671              :             END DO
    1672      7184008 :             tmp_max = tmp_max*max_contraction
    1673      7184008 :             IF (tmp_max < eps_schwarz) THEN
    1674              :                RETURN
    1675              :             END IF
    1676              : 
    1677     23792400 :             DO j = 1, ncob
    1678     18514744 :                p1 = (j - 1)*ncoa
    1679     45156782 :                DO i = 1, ncoa
    1680     21364382 :                   p2 = (p1 + i - 1)*ncod
    1681    141134671 :                   DO l = 1, ncod
    1682    101255545 :                      p3 = (p2 + l - 1)*ncoc
    1683    317301752 :                      DO k = 1, ncoc
    1684    194681825 :                         p4 = p3 + k
    1685    295937370 :                         work2(i, j, k, l) = work(p4)
    1686              :                      END DO
    1687              :                   END DO
    1688              :                END DO
    1689              :             END DO
    1690              :          CASE (5)
    1691     95177433 :             DO i = 1, nproducts
    1692    272297700 :                A = pgf_product_list(i)%ra
    1693    272297700 :                B = pgf_product_list(i)%rb
    1694    272297700 :                C = pgf_product_list(i)%rc
    1695    272297700 :                D = pgf_product_list(i)%rd
    1696     68074425 :                Rho = pgf_product_list(i)%Rho
    1697     68074425 :                RhoInv = pgf_product_list(i)%RhoInv
    1698    272297700 :                P = pgf_product_list(i)%P
    1699    272297700 :                Q = pgf_product_list(i)%Q
    1700    272297700 :                W = pgf_product_list(i)%W
    1701    272297700 :                AB = pgf_product_list(i)%AB
    1702    272297700 :                CD = pgf_product_list(i)%CD
    1703     68074425 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1704    259733813 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1705              : 
    1706              :                CALL build_quartet_data(libint, C, D, A, B, &
    1707              :                                        EtaInv, ZetaInv, ZetapEtaInv, Rho, &
    1708     68074425 :                                        Q, P, W, m_max, F)
    1709              : 
    1710     68074425 :                CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
    1711   1120893819 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1712     95177433 :                neris = neris + mysize
    1713              :             END DO
    1714    580683249 :             DO i = 1, mysize
    1715    580683249 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1716              :             END DO
    1717     27103008 :             tmp_max = tmp_max*max_contraction
    1718     27103008 :             IF (tmp_max < eps_schwarz) THEN
    1719              :                RETURN
    1720              :             END IF
    1721              : 
    1722     48136658 :             DO k = 1, ncoc
    1723     28633752 :                p1 = (k - 1)*ncod
    1724     81089772 :                DO l = 1, ncod
    1725     32953114 :                   p2 = (p1 + l - 1)*ncoa
    1726    211791164 :                   DO i = 1, ncoa
    1727    150204298 :                      p3 = (p2 + i - 1)*ncob
    1728    623810494 :                      DO j = 1, ncob
    1729    440653082 :                         p4 = p3 + j
    1730    590857380 :                         work2(i, j, k, l) = work(p4)
    1731              :                      END DO
    1732              :                   END DO
    1733              :                END DO
    1734              :             END DO
    1735              :          CASE (6)
    1736     41005851 :             DO i = 1, nproducts
    1737    106144748 :                A = pgf_product_list(i)%ra
    1738    106144748 :                B = pgf_product_list(i)%rb
    1739    106144748 :                C = pgf_product_list(i)%rc
    1740    106144748 :                D = pgf_product_list(i)%rd
    1741     26536187 :                Rho = pgf_product_list(i)%Rho
    1742     26536187 :                RhoInv = pgf_product_list(i)%RhoInv
    1743    106144748 :                P = pgf_product_list(i)%P
    1744    106144748 :                Q = pgf_product_list(i)%Q
    1745    106144748 :                W = pgf_product_list(i)%W
    1746    106144748 :                AB = pgf_product_list(i)%AB
    1747    106144748 :                CD = pgf_product_list(i)%CD
    1748     26536187 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1749     95767000 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1750              : 
    1751              :                CALL build_quartet_data(libint, C, D, B, A, &
    1752              :                                        EtaInv, ZetaInv, ZetapEtaInv, Rho, &
    1753     26536187 :                                        Q, P, W, m_max, F)
    1754              : 
    1755     26536187 :                CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
    1756    380269096 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1757     41005851 :                neris = neris + mysize
    1758              :             END DO
    1759    229105358 :             DO i = 1, mysize
    1760    229105358 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1761              :             END DO
    1762     14469664 :             tmp_max = tmp_max*max_contraction
    1763     14469664 :             IF (tmp_max < eps_schwarz) THEN
    1764              :                RETURN
    1765              :             END IF
    1766              : 
    1767     21864135 :             DO k = 1, ncoc
    1768     12747523 :                p1 = (k - 1)*ncod
    1769     36742816 :                DO l = 1, ncod
    1770     14878681 :                   p2 = (p1 + l - 1)*ncob
    1771    102114770 :                   DO j = 1, ncob
    1772     74488566 :                      p3 = (p2 + j - 1)*ncoa
    1773    245163885 :                      DO i = 1, ncoa
    1774    155796638 :                         p4 = p3 + i
    1775    230285204 :                         work2(i, j, k, l) = work(p4)
    1776              :                      END DO
    1777              :                   END DO
    1778              :                END DO
    1779              :             END DO
    1780              :          CASE (7)
    1781     17296587 :             DO i = 1, nproducts
    1782     47427728 :                A = pgf_product_list(i)%ra
    1783     47427728 :                B = pgf_product_list(i)%rb
    1784     47427728 :                C = pgf_product_list(i)%rc
    1785     47427728 :                D = pgf_product_list(i)%rd
    1786     11856932 :                Rho = pgf_product_list(i)%Rho
    1787     11856932 :                RhoInv = pgf_product_list(i)%RhoInv
    1788     47427728 :                P = pgf_product_list(i)%P
    1789     47427728 :                Q = pgf_product_list(i)%Q
    1790     47427728 :                W = pgf_product_list(i)%W
    1791     47427728 :                AB = pgf_product_list(i)%AB
    1792     47427728 :                CD = pgf_product_list(i)%CD
    1793     11856932 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1794     63542942 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1795              : 
    1796              :                CALL build_quartet_data(libint, D, C, A, B, &
    1797              :                                        EtaInv, ZetaInv, ZetapEtaInv, Rho, &
    1798     11856932 :                                        Q, P, W, m_max, F)
    1799              : 
    1800     11856932 :                CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
    1801    566157832 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1802     17296587 :                neris = neris + mysize
    1803              :             END DO
    1804    354239258 :             DO i = 1, mysize
    1805    354239258 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1806              :             END DO
    1807      5439655 :             tmp_max = tmp_max*max_contraction
    1808      5439655 :             IF (tmp_max < eps_schwarz) THEN
    1809              :                RETURN
    1810              :             END IF
    1811              : 
    1812     16676220 :             DO l = 1, ncod
    1813     12930388 :                p1 = (l - 1)*ncoc
    1814     31285376 :                DO k = 1, ncoc
    1815     14609156 :                   p2 = (p1 + k - 1)*ncoa
    1816    103468368 :                   DO i = 1, ncoa
    1817     75928824 :                      p3 = (p2 + i - 1)*ncob
    1818    368420358 :                      DO j = 1, ncob
    1819    277882378 :                         p4 = p3 + j
    1820    353811202 :                         work2(i, j, k, l) = work(p4)
    1821              :                      END DO
    1822              :                   END DO
    1823              :                END DO
    1824              :             END DO
    1825              :          CASE (8)
    1826      5385265 :             DO i = 1, nproducts
    1827     12700576 :                A = pgf_product_list(i)%ra
    1828     12700576 :                B = pgf_product_list(i)%rb
    1829     12700576 :                C = pgf_product_list(i)%rc
    1830     12700576 :                D = pgf_product_list(i)%rd
    1831      3175144 :                Rho = pgf_product_list(i)%Rho
    1832      3175144 :                RhoInv = pgf_product_list(i)%RhoInv
    1833     12700576 :                P = pgf_product_list(i)%P
    1834     12700576 :                Q = pgf_product_list(i)%Q
    1835     12700576 :                W = pgf_product_list(i)%W
    1836     12700576 :                AB = pgf_product_list(i)%AB
    1837     12700576 :                CD = pgf_product_list(i)%CD
    1838      3175144 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1839     18037058 :                F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1840              : 
    1841              :                CALL build_quartet_data(libint, D, C, B, A, &
    1842              :                                        EtaInv, ZetaInv, ZetapEtaInv, Rho, &
    1843      3175144 :                                        Q, P, W, m_max, F)
    1844              : 
    1845      3175144 :                CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
    1846    153635212 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1847      5385265 :                neris = neris + mysize
    1848              :             END DO
    1849    113328321 :             DO i = 1, mysize
    1850    113328321 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1851              :             END DO
    1852      2210121 :             tmp_max = tmp_max*max_contraction
    1853      2210121 :             IF (tmp_max < eps_schwarz) THEN
    1854              :                RETURN
    1855              :             END IF
    1856              : 
    1857    146892164 :             DO l = 1, ncod
    1858      5858374 :                p1 = (l - 1)*ncoc
    1859     13425511 :                DO k = 1, ncoc
    1860      5918694 :                   p2 = (p1 + k - 1)*ncob
    1861     48359116 :                   DO j = 1, ncob
    1862     36582048 :                      p3 = (p2 + j - 1)*ncoa
    1863    128266872 :                      DO i = 1, ncoa
    1864     85766130 :                         p4 = p3 + i
    1865    122348178 :                         work2(i, j, k, l) = work(p4)
    1866              :                      END DO
    1867              :                   END DO
    1868              :                END DO
    1869              :             END DO
    1870              :          END SELECT
    1871              :       ELSE
    1872    145207422 :          DO i = 1, nproducts
    1873    433674868 :             A = pgf_product_list(i)%ra
    1874    433674868 :             B = pgf_product_list(i)%rb
    1875    433674868 :             C = pgf_product_list(i)%rc
    1876    433674868 :             D = pgf_product_list(i)%rd
    1877    108418717 :             Rho = pgf_product_list(i)%Rho
    1878    108418717 :             RhoInv = pgf_product_list(i)%RhoInv
    1879    433674868 :             P = pgf_product_list(i)%P
    1880    433674868 :             Q = pgf_product_list(i)%Q
    1881    433674868 :             W = pgf_product_list(i)%W
    1882    108418717 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1883    216837434 :             F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
    1884              : 
    1885              :             CALL build_quartet_data(libint, A, B, C, D, & !todo: check
    1886              :                                     ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1887    108418717 :                                     P, Q, W, m_max, F)
    1888    108418717 :             work(1) = work(1) + F(1)
    1889    145207422 :             neris = neris + mysize
    1890              :          END DO
    1891     36788705 :          work2(1, 1, 1, 1) = work(1)
    1892     36788705 :          tmp_max = max_contraction*ABS(work(1))
    1893     36788705 :          IF (tmp_max < eps_schwarz) RETURN
    1894              :       END IF
    1895              : 
    1896    131403879 :       IF (tmp_max < eps_schwarz) RETURN
    1897   8321252464 :       primitives_tmp = 0.0_dp
    1898              : 
    1899              :       CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
    1900              :                     n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2, &
    1901              :                     sphi_a, &
    1902              :                     sphi_b, &
    1903              :                     sphi_c, &
    1904              :                     sphi_d, &
    1905              :                     primitives_tmp, &
    1906    131403879 :                     buffer1, buffer2)
    1907              : 
    1908              :       primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1909              :                  s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1910              :                  s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1911              :                  s_offset_d + 1:s_offset_d + nsod*nl_d) = &
    1912              :          primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1913              :                     s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1914              :                     s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1915   8321252464 :                     s_offset_d + 1:s_offset_d + nsod*nl_d) + primitives_tmp(:, :, :, :)
    1916              : 
    1917              :    END SUBROUTINE evaluate_eri
    1918              : 
    1919              : ! **************************************************************************************************
    1920              : !> \brief Fill data structure used in libint
    1921              : !> \param libint ...
    1922              : !> \param A ...
    1923              : !> \param B ...
    1924              : !> \param C ...
    1925              : !> \param D ...
    1926              : !> \param ZetaInv ...
    1927              : !> \param EtaInv ...
    1928              : !> \param ZetapEtaInv ...
    1929              : !> \param Rho ...
    1930              : !> \param P ...
    1931              : !> \param Q ...
    1932              : !> \param W ...
    1933              : !> \param m_max ...
    1934              : !> \param F ...
    1935              : !> \par History
    1936              : !>      03.2007 created [Manuel Guidon]
    1937              : !> \author Manuel Guidon
    1938              : ! **************************************************************************************************
    1939    474732024 :    SUBROUTINE build_quartet_data(libint, A, B, C, D, &
    1940              :                                  ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1941    474732024 :                                  P, Q, W, m_max, F)
    1942              : 
    1943              :       TYPE(cp_libint_t)                                  :: libint
    1944              :       REAL(KIND=dp), INTENT(IN)                          :: A(3), B(3), C(3), D(3)
    1945              :       REAL(dp), INTENT(IN)                               :: ZetaInv, EtaInv, ZetapEtaInv, Rho, P(3), &
    1946              :                                                             Q(3), W(3)
    1947              :       INTEGER, INTENT(IN)                                :: m_max
    1948              :       REAL(KIND=dp), DIMENSION(:)                        :: F
    1949              : 
    1950    474732024 :       CALL cp_libint_set_params_eri(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, P, Q, W, m_max, F)
    1951    474732024 :    END SUBROUTINE build_quartet_data
    1952              : 
    1953              : END MODULE hfx_libint_interface
        

Generated by: LCOV version 2.0-1