LCOV - code coverage report
Current view: top level - src - hfx_libint_interface.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b8e0b09) Lines: 951 983 96.7 %
Date: 2024-08-31 06:31:37 Functions: 6 6 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief 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    81977862 :    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    81977862 :       Zeta = Zeta_A + Zeta_B
      98    81977862 :       ZetaInv = 1.0_dp/Zeta
      99    81977862 :       Eta = Zeta_C + Zeta_D
     100    81977862 :       EtaInv = 1.0_dp/Eta
     101    81977862 :       ZetapEtaInv = Zeta + Eta
     102    81977862 :       ZetapEtaInv = 1.0_dp/ZetapEtaInv
     103    81977862 :       Rho = Zeta*Eta*ZetapEtaInv
     104    81977862 :       RhoInv = 1.0_dp/Rho
     105             : 
     106   327911448 :       DO i = 1, 3
     107   245933586 :          P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv
     108   245933586 :          Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv
     109   245933586 :          AB(i) = A(i) - B(i)
     110   245933586 :          CD(i) = C(i) - D(i)
     111   245933586 :          PQ(i) = P(i) - Q(i)
     112   327911448 :          W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv
     113             :       END DO
     114             : 
     115   327911448 :       AB2 = DOT_PRODUCT(AB, AB)
     116   327911448 :       CD2 = DOT_PRODUCT(CD, CD)
     117   327911448 :       PQ2 = DOT_PRODUCT(PQ, PQ)
     118             : 
     119    81977862 :       S1234 = EXP((-Zeta_A*Zeta_B*ZetaInv*AB2) + (-Zeta_C*Zeta_D*EtaInv*CD2))
     120    81977862 :       T = Rho*PQ2
     121             : 
     122    94626900 :       SELECT CASE (potential_parameter%potential_type)
     123             :       CASE (do_potential_truncated)
     124    12649038 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     125    12649038 :          R1 = R11
     126    12649038 :          R2 = R22
     127    12649038 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     128           0 :             RETURN
     129             :          END IF
     130    12649038 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     131    12649038 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     132    12649038 :          factor = 2.0_dp*Pi*RhoInv
     133             :       CASE (do_potential_coulomb)
     134    62255592 :          CALL fgamma(m_max, T, F(1))
     135    62255592 :          factor = 2.0_dp*Pi*RhoInv
     136             :       CASE (do_potential_short)
     137     3082116 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     138     3082116 :          R1 = R11
     139     3082116 :          R2 = R22
     140     3082116 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     141             :             RETURN
     142             :          END IF
     143     3082116 :          CALL fgamma(m_max, T, F(1))
     144     3082116 :          omega2 = potential_parameter%omega**2
     145     3082116 :          omega_corr2 = omega2/(omega2 + Rho)
     146     3082116 :          omega_corr = SQRT(omega_corr2)
     147     3082116 :          T = T*omega_corr2
     148     3082116 :          CALL fgamma(m_max, T, Fm)
     149     3082116 :          tmp = -omega_corr
     150    11606112 :          DO i = 1, m_max + 1
     151     8523996 :             F(i) = F(i) + Fm(i)*tmp
     152    11606112 :             tmp = tmp*omega_corr2
     153             :          END DO
     154     3082116 :          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      939906 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     184      939906 :          R1 = R11
     185      939906 :          R2 = R22
     186      939906 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     187             :             RETURN
     188             :          END IF
     189      939906 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     190      939906 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     191             : 
     192             :          ! Coulomb
     193      939906 :          CALL fgamma(m_max, T, Fm)
     194             : 
     195     3380268 :          DO i = 1, m_max + 1
     196             :             F(i) = F(i)*(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
     197     3380268 :                    Fm(i)*potential_parameter%scale_longrange
     198             :          END DO
     199             : 
     200             :          ! longrange
     201      939906 :          omega2 = potential_parameter%omega**2
     202      939906 :          omega_corr2 = omega2/(omega2 + Rho)
     203      939906 :          omega_corr = SQRT(omega_corr2)
     204      939906 :          T = T*omega_corr2
     205      939906 :          CALL fgamma(m_max, T, Fm)
     206      939906 :          tmp = omega_corr
     207     3380268 :          DO i = 1, m_max + 1
     208     2440362 :             F(i) = F(i) + Fm(i)*tmp*potential_parameter%scale_longrange
     209     3380268 :             tmp = tmp*omega_corr2
     210             :          END DO
     211      939906 :          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    95809206 :          IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
     256             :       END SELECT
     257             : 
     258    81977862 :       tmp = (Pi*ZetapEtaInv)**3
     259    81977862 :       factor = factor*S1234*SQRT(tmp)
     260             : 
     261   316199892 :       DO i = 1, m_max + 1
     262   316199892 :          F(i) = F(i)*factor
     263             :       END DO
     264             : 
     265    81977862 :       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    93417313 :    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    93417313 :                                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    93417313 :                                           ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F)
     313             : 
     314    93417313 :    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    31781452 :    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    31781452 :                                  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    31781452 :                                  sphi_a, sphi_b, sphi_c, sphi_d, &
     388    31781452 :                                  work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
     389    31781452 :                                  use_virial, work_virial, work2_virial, primitives_tmp_virial, &
     390    31781452 :                                  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    31781452 :       m_max = n_a + n_b + n_c + n_d
     433    31781452 :       m_max = m_max + 1
     434    31781452 :       mysize = ncoa*ncob*ncoc*ncod
     435    63562904 :       a_mysize = mysize
     436             : 
     437  4152737008 :       work = 0.0_dp
     438  8001517108 :       work2 = 0.0_dp
     439             : 
     440    31781452 :       IF (use_virial) THEN
     441   850508488 :          work_virial = 0.0_dp
     442  1465410088 :          work2_virial = 0.0_dp
     443             :       END IF
     444             : 
     445    31781452 :       perm_case = 1
     446    31781452 :       IF (n_a < n_b) THEN
     447     2917159 :          perm_case = perm_case + 1
     448             :       END IF
     449    31781452 :       IF (n_c < n_d) THEN
     450     5032807 :          perm_case = perm_case + 2
     451             :       END IF
     452    31781452 :       IF (n_a + n_b > n_c + n_d) THEN
     453     5723891 :          perm_case = perm_case + 4
     454             :       END IF
     455             : 
     456             :       SELECT CASE (perm_case)
     457             :       CASE (1)
     458    77913540 :          DO i = 1, nproducts
     459   229652608 :             A = pgf_product_list(i)%ra
     460   229652608 :             B = pgf_product_list(i)%rb
     461   229652608 :             C = pgf_product_list(i)%rc
     462   229652608 :             D = pgf_product_list(i)%rd
     463    57413152 :             Rho = pgf_product_list(i)%Rho
     464    57413152 :             RhoInv = pgf_product_list(i)%RhoInv
     465   229652608 :             P = pgf_product_list(i)%P
     466   229652608 :             Q = pgf_product_list(i)%Q
     467   229652608 :             W = pgf_product_list(i)%W
     468   229652608 :             AB = pgf_product_list(i)%AB
     469   229652608 :             CD = pgf_product_list(i)%CD
     470    57413152 :             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    57413152 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     476    57413152 :             CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize)
     477   229652608 :             DO k = 4, 6
     478  1745316982 :                DO j = 1, mysize
     479             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     480             :                                                work_forces(j, k + 3) + &
     481  1687903830 :                                                work_forces(j, k + 6))
     482             :                END DO
     483             :             END DO
     484   746370976 :             DO k = 1, 12
     485  6809028472 :                DO j = 1, mysize
     486  6751615320 :                   work(j, k) = work(j, k) + work_forces(j, k)
     487             :                END DO
     488             :             END DO
     489    57413152 :             neris = neris + 12*mysize
     490    77913540 :             IF (use_virial) THEN
     491    10723659 :                CALL real_to_scaled(scoord(1:3), A, cell)
     492    10723659 :                CALL real_to_scaled(scoord(4:6), B, cell)
     493    10723659 :                CALL real_to_scaled(scoord(7:9), C, cell)
     494    10723659 :                CALL real_to_scaled(scoord(10:12), D, cell)
     495   139407567 :                DO k = 1, 12
     496  1164325875 :                   DO j = 1, mysize
     497  4228357140 :                      DO m = 1, 3
     498  4099673232 :                         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   266505044 :          DO n = 1, 12
     506             :             tmp_max = 0.0_dp
     507  2079639288 :             DO i = 1, mysize
     508  2079639288 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     509             :             END DO
     510   246004656 :             tmp_max = tmp_max*max_contraction
     511   246004656 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     512             : 
     513   605466056 :             DO i = 1, ncoa
     514   338961012 :                p1 = (i - 1)*ncob
     515   957892032 :                DO j = 1, ncob
     516   372926364 :                   p2 = (p1 + j - 1)*ncoc
     517  1710667248 :                   DO k = 1, ncoc
     518   998779872 :                      p3 = (p2 + k - 1)*ncod
     519  3205340868 :                      DO l = 1, ncod
     520  1833634632 :                         p4 = p3 + l
     521  2832414504 :                         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    20500388 :          IF (use_virial) THEN
     528     6234475 :             DO n = 1, 12
     529             :                tmp_max_virial = 0.0_dp
     530   112742940 :                DO i = 1, mysize
     531             :                   tmp_max_virial = MAX(tmp_max_virial, &
     532   112742940 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     533             :                END DO
     534     5754900 :                tmp_max_virial = tmp_max_virial*max_contraction
     535     5754900 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     536             : 
     537    16529839 :                DO i = 1, ncoa
     538    10295364 :                   p1 = (i - 1)*ncob
     539    29968524 :                   DO j = 1, ncob
     540    13918260 :                      p2 = (p1 + j - 1)*ncoc
     541    69611136 :                      DO k = 1, ncoc
     542    45397512 :                         p3 = (p2 + k - 1)*ncod
     543   166303812 :                         DO l = 1, ncod
     544   106988040 :                            p4 = p3 + l
     545   473349672 :                            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     5366299 :          DO i = 1, nproducts
     554    16818556 :             A = pgf_product_list(i)%ra
     555    16818556 :             B = pgf_product_list(i)%rb
     556    16818556 :             C = pgf_product_list(i)%rc
     557    16818556 :             D = pgf_product_list(i)%rd
     558     4204639 :             Rho = pgf_product_list(i)%Rho
     559     4204639 :             RhoInv = pgf_product_list(i)%RhoInv
     560    16818556 :             P = pgf_product_list(i)%P
     561    16818556 :             Q = pgf_product_list(i)%Q
     562    16818556 :             W = pgf_product_list(i)%W
     563    16818556 :             AB = pgf_product_list(i)%AB
     564    16818556 :             CD = pgf_product_list(i)%CD
     565     4204639 :             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     4204639 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     571     4204639 :             CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize)
     572    16818556 :             DO k = 4, 6
     573   299056819 :                DO j = 1, mysize
     574             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     575             :                                                work_forces(j, k + 3) + &
     576   294852180 :                                                work_forces(j, k + 6))
     577             :                END DO
     578             :             END DO
     579    54660307 :             DO k = 1, 12
     580  1183613359 :                DO j = 1, mysize
     581  1179408720 :                   work(j, k) = work(j, k) + work_forces(j, k)
     582             :                END DO
     583             :             END DO
     584     4204639 :             neris = neris + 12*mysize
     585     5366299 :             IF (use_virial) THEN
     586      674562 :                CALL real_to_scaled(scoord(1:3), B, cell)
     587      674562 :                CALL real_to_scaled(scoord(4:6), A, cell)
     588      674562 :                CALL real_to_scaled(scoord(7:9), C, cell)
     589      674562 :                CALL real_to_scaled(scoord(10:12), D, cell)
     590     8769306 :                DO k = 1, 12
     591   180424398 :                   DO j = 1, mysize
     592   694715112 :                      DO m = 1, 3
     593   686620368 :                         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    15101580 :          DO n = 1, 12
     601             :             tmp_max = 0.0_dp
     602   361438236 :             DO i = 1, mysize
     603   361438236 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     604             :             END DO
     605    13939920 :             tmp_max = tmp_max*max_contraction
     606    13939920 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     607             : 
     608    59319768 :             DO j = 1, ncob
     609    44218188 :                p1 = (j - 1)*ncoa
     610   103801032 :                DO i = 1, ncoa
     611    45642924 :                   p2 = (p1 + i - 1)*ncoc
     612   253822932 :                   DO k = 1, ncoc
     613   163961820 :                      p3 = (p2 + k - 1)*ncod
     614   557103060 :                      DO l = 1, ncod
     615   347498316 :                         p4 = p3 + l
     616   511460136 :                         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     1161660 :          IF (use_virial) THEN
     623     1274208 :             DO n = 1, 12
     624             :                tmp_max_virial = 0.0_dp
     625    34689024 :                DO i = 1, mysize
     626             :                   tmp_max_virial = MAX(tmp_max_virial, &
     627    34689024 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     628             :                END DO
     629     1176192 :                tmp_max_virial = tmp_max_virial*max_contraction
     630     1176192 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     631             : 
     632     5079264 :                DO j = 1, ncob
     633     3805056 :                   p1 = (j - 1)*ncoa
     634     8970624 :                   DO i = 1, ncoa
     635     3989376 :                      p2 = (p1 + i - 1)*ncoc
     636    22275504 :                      DO k = 1, ncoc
     637    14481072 :                         p3 = (p2 + k - 1)*ncod
     638    51983280 :                         DO l = 1, ncod
     639    33512832 :                            p4 = p3 + l
     640   148532400 :                            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    15779050 :          DO i = 1, nproducts
     649    48090580 :             A = pgf_product_list(i)%ra
     650    48090580 :             B = pgf_product_list(i)%rb
     651    48090580 :             C = pgf_product_list(i)%rc
     652    48090580 :             D = pgf_product_list(i)%rd
     653    12022645 :             Rho = pgf_product_list(i)%Rho
     654    12022645 :             RhoInv = pgf_product_list(i)%RhoInv
     655    48090580 :             P = pgf_product_list(i)%P
     656    48090580 :             Q = pgf_product_list(i)%Q
     657    48090580 :             W = pgf_product_list(i)%W
     658    48090580 :             AB = pgf_product_list(i)%AB
     659    48090580 :             CD = pgf_product_list(i)%CD
     660    12022645 :             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    12022645 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     666    12022645 :             CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize)
     667    48090580 :             DO k = 4, 6
     668   420409339 :                DO j = 1, mysize
     669             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     670             :                                                work_forces(j, k + 3) + &
     671   408386694 :                                                work_forces(j, k + 6))
     672             :                END DO
     673             :             END DO
     674   156294385 :             DO k = 1, 12
     675  1645569421 :                DO j = 1, mysize
     676  1633546776 :                   work(j, k) = work(j, k) + work_forces(j, k)
     677             :                END DO
     678             :             END DO
     679    12022645 :             neris = neris + 12*mysize
     680    15779050 :             IF (use_virial) THEN
     681     1825592 :                CALL real_to_scaled(scoord(1:3), A, cell)
     682     1825592 :                CALL real_to_scaled(scoord(4:6), B, cell)
     683     1825592 :                CALL real_to_scaled(scoord(7:9), D, cell)
     684     1825592 :                CALL real_to_scaled(scoord(10:12), C, cell)
     685    23732696 :                DO k = 1, 12
     686   205912496 :                   DO j = 1, mysize
     687   750626304 :                      DO m = 1, 3
     688   728719200 :                         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    48833265 :          DO n = 1, 12
     696             :             tmp_max = 0.0_dp
     697   519166620 :             DO i = 1, mysize
     698   519166620 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     699             :             END DO
     700    45076860 :             tmp_max = tmp_max*max_contraction
     701    45076860 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     702             : 
     703   126963021 :             DO i = 1, ncoa
     704    78129756 :                p1 = (i - 1)*ncob
     705   208514556 :                DO j = 1, ncob
     706    85307940 :                   p2 = (p1 + j - 1)*ncod
     707   492029568 :                   DO l = 1, ncod
     708   328591872 :                      p3 = (p2 + l - 1)*ncoc
     709   887989572 :                      DO k = 1, ncoc
     710   474089760 :                         p4 = p3 + k
     711   802681632 :                         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     3756405 :          IF (use_virial) THEN
     718     1829672 :             DO n = 1, 12
     719             :                tmp_max_virial = 0.0_dp
     720    32067096 :                DO i = 1, mysize
     721             :                   tmp_max_virial = MAX(tmp_max_virial, &
     722    32067096 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     723             :                END DO
     724     1688928 :                tmp_max_virial = tmp_max_virial*max_contraction
     725     1688928 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     726             : 
     727     5285648 :                DO i = 1, ncoa
     728     3455976 :                   p1 = (i - 1)*ncob
     729     9372720 :                   DO j = 1, ncob
     730     4227816 :                      p2 = (p1 + j - 1)*ncod
     731    25807560 :                      DO l = 1, ncod
     732    18123768 :                         p3 = (p2 + l - 1)*ncoc
     733    52729752 :                         DO k = 1, ncoc
     734    30378168 :                            p4 = p3 + k
     735   139636440 :                            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     2724516 :          DO i = 1, nproducts
     744     8341632 :             A = pgf_product_list(i)%ra
     745     8341632 :             B = pgf_product_list(i)%rb
     746     8341632 :             C = pgf_product_list(i)%rc
     747     8341632 :             D = pgf_product_list(i)%rd
     748     2085408 :             Rho = pgf_product_list(i)%Rho
     749     2085408 :             RhoInv = pgf_product_list(i)%RhoInv
     750     8341632 :             P = pgf_product_list(i)%P
     751     8341632 :             Q = pgf_product_list(i)%Q
     752     8341632 :             W = pgf_product_list(i)%W
     753     8341632 :             AB = pgf_product_list(i)%AB
     754     8341632 :             CD = pgf_product_list(i)%CD
     755     2085408 :             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     2085408 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     760     2085408 :             CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize)
     761     8341632 :             DO k = 4, 6
     762   112833549 :                DO j = 1, mysize
     763             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     764             :                                                work_forces(j, k + 3) + &
     765   110748141 :                                                work_forces(j, k + 6))
     766             :                END DO
     767             :             END DO
     768    27110304 :             DO k = 1, 12
     769   445077972 :                DO j = 1, mysize
     770   442992564 :                   work(j, k) = work(j, k) + work_forces(j, k)
     771             :                END DO
     772             :             END DO
     773     2085408 :             neris = neris + 12*mysize
     774     2724516 :             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     8308404 :          DO n = 1, 12
     790             :             tmp_max = 0.0_dp
     791   153053712 :             DO i = 1, mysize
     792   153053712 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     793             :             END DO
     794     7669296 :             tmp_max = tmp_max*max_contraction
     795     7669296 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     796             : 
     797    32305140 :             DO j = 1, ncob
     798    23996736 :                p1 = (j - 1)*ncoa
     799    56802960 :                DO i = 1, ncoa
     800    25136928 :                   p2 = (p1 + i - 1)*ncod
     801   145733616 :                   DO l = 1, ncod
     802    96599952 :                      p3 = (p2 + l - 1)*ncoc
     803   267121296 :                      DO k = 1, ncoc
     804   145384416 :                         p4 = p3 + k
     805   241984368 :                         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      639108 :          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    16486302 :          DO i = 1, nproducts
     838    49716668 :             A = pgf_product_list(i)%ra
     839    49716668 :             B = pgf_product_list(i)%rb
     840    49716668 :             C = pgf_product_list(i)%rc
     841    49716668 :             D = pgf_product_list(i)%rd
     842    12429167 :             Rho = pgf_product_list(i)%Rho
     843    12429167 :             RhoInv = pgf_product_list(i)%RhoInv
     844    49716668 :             P = pgf_product_list(i)%P
     845    49716668 :             Q = pgf_product_list(i)%Q
     846    49716668 :             W = pgf_product_list(i)%W
     847    49716668 :             AB = pgf_product_list(i)%AB
     848    49716668 :             CD = pgf_product_list(i)%CD
     849    12429167 :             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    12429167 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     854    12429167 :             CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize)
     855    49716668 :             DO k = 4, 6
     856   470763956 :                DO j = 1, mysize
     857             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     858             :                                                work_forces(j, k + 3) + &
     859   458334789 :                                                work_forces(j, k + 6))
     860             :                END DO
     861             :             END DO
     862   161579171 :             DO k = 1, 12
     863  1845768323 :                DO j = 1, mysize
     864  1833339156 :                   work(j, k) = work(j, k) + work_forces(j, k)
     865             :                END DO
     866             :             END DO
     867    12429167 :             neris = neris + 12*mysize
     868    16486302 :             IF (use_virial) THEN
     869     1995744 :                CALL real_to_scaled(scoord(1:3), C, cell)
     870     1995744 :                CALL real_to_scaled(scoord(4:6), D, cell)
     871     1995744 :                CALL real_to_scaled(scoord(7:9), A, cell)
     872     1995744 :                CALL real_to_scaled(scoord(10:12), B, cell)
     873    25944672 :                DO k = 1, 12
     874   267769560 :                   DO j = 1, mysize
     875   991248480 :                      DO m = 1, 3
     876   967299552 :                         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    52742755 :          DO n = 1, 12
     884             :             tmp_max = 0.0_dp
     885   576919308 :             DO i = 1, mysize
     886   576919308 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     887             :             END DO
     888    48685620 :             tmp_max = tmp_max*max_contraction
     889    48685620 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     890             : 
     891   118134931 :             DO k = 1, ncoc
     892    65392176 :                p1 = (k - 1)*ncod
     893   183078684 :                DO l = 1, ncod
     894    69000888 :                   p2 = (p1 + l - 1)*ncoa
     895   393086112 :                   DO i = 1, ncoa
     896   258693048 :                      p3 = (p2 + i - 1)*ncob
     897   855927624 :                      DO j = 1, ncob
     898   528233688 :                         p4 = p3 + j
     899   786926736 :                         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     4057135 :          IF (use_virial) THEN
     906     2252510 :             DO n = 1, 12
     907             :                tmp_max_virial = 0.0_dp
     908    44908116 :                DO i = 1, mysize
     909             :                   tmp_max_virial = MAX(tmp_max_virial, &
     910    44908116 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     911             :                END DO
     912     2079240 :                tmp_max_virial = tmp_max_virial*max_contraction
     913     2079240 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     914             : 
     915     5748830 :                DO k = 1, ncoc
     916     3496320 :                   p1 = (k - 1)*ncod
     917     9492504 :                   DO l = 1, ncod
     918     3916944 :                      p2 = (p1 + l - 1)*ncoa
     919    22983084 :                      DO i = 1, ncoa
     920    15569820 :                         p3 = (p2 + i - 1)*ncob
     921    62315640 :                         DO j = 1, ncob
     922    42828876 :                            p4 = p3 + j
     923   186885324 :                            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     4166113 :          DO i = 1, nproducts
     932    12546604 :             A = pgf_product_list(i)%ra
     933    12546604 :             B = pgf_product_list(i)%rb
     934    12546604 :             C = pgf_product_list(i)%rc
     935    12546604 :             D = pgf_product_list(i)%rd
     936     3136651 :             Rho = pgf_product_list(i)%Rho
     937     3136651 :             RhoInv = pgf_product_list(i)%RhoInv
     938    12546604 :             P = pgf_product_list(i)%P
     939    12546604 :             Q = pgf_product_list(i)%Q
     940    12546604 :             W = pgf_product_list(i)%W
     941    12546604 :             AB = pgf_product_list(i)%AB
     942    12546604 :             CD = pgf_product_list(i)%CD
     943     3136651 :             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     3136651 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     949     3136651 :             CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize)
     950    12546604 :             DO k = 4, 6
     951   118457794 :                DO j = 1, mysize
     952             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     953             :                                                work_forces(j, k + 3) + &
     954   115321143 :                                                work_forces(j, k + 6))
     955             :                END DO
     956             :             END DO
     957    40776463 :             DO k = 1, 12
     958   464421223 :                DO j = 1, mysize
     959   461284572 :                   work(j, k) = work(j, k) + work_forces(j, k)
     960             :                END DO
     961             :             END DO
     962     3136651 :             neris = neris + 12*mysize
     963     4166113 :             IF (use_virial) THEN
     964      381932 :                CALL real_to_scaled(scoord(1:3), C, cell)
     965      381932 :                CALL real_to_scaled(scoord(4:6), D, cell)
     966      381932 :                CALL real_to_scaled(scoord(7:9), B, cell)
     967      381932 :                CALL real_to_scaled(scoord(10:12), A, cell)
     968     4965116 :                DO k = 1, 12
     969    51073376 :                   DO j = 1, mysize
     970   189016224 :                      DO m = 1, 3
     971   184433040 :                         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    13383006 :          DO n = 1, 12
     979             :             tmp_max = 0.0_dp
     980   160046064 :             DO i = 1, mysize
     981   160046064 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     982             :             END DO
     983    12353544 :             tmp_max = tmp_max*max_contraction
     984    12353544 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     985             : 
     986    28897302 :             DO k = 1, ncoc
     987    15514296 :                p1 = (k - 1)*ncod
     988    45660072 :                DO l = 1, ncod
     989    17792232 :                   p2 = (p1 + l - 1)*ncob
     990   111098568 :                   DO j = 1, ncob
     991    77792040 :                      p3 = (p2 + j - 1)*ncoa
     992   243276792 :                      DO i = 1, ncoa
     993   147692520 :                         p4 = p3 + i
     994   225484560 :                         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     1029462 :          IF (use_virial) THEN
    1001      836940 :             DO n = 1, 12
    1002             :                tmp_max_virial = 0.0_dp
    1003    16324416 :                DO i = 1, mysize
    1004             :                   tmp_max_virial = MAX(tmp_max_virial, &
    1005    16324416 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1006             :                END DO
    1007      772560 :                tmp_max_virial = tmp_max_virial*max_contraction
    1008      772560 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1009             : 
    1010     1965852 :                DO k = 1, ncoc
    1011     1128912 :                   p1 = (k - 1)*ncod
    1012     3325296 :                   DO l = 1, ncod
    1013     1423824 :                      p2 = (p1 + l - 1)*ncob
    1014     9552144 :                      DO j = 1, ncob
    1015     6999408 :                         p3 = (p2 + j - 1)*ncoa
    1016    23975088 :                         DO i = 1, ncoa
    1017    15551856 :                            p4 = p3 + i
    1018    69206832 :                            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     2450212 :          DO i = 1, nproducts
    1027     7599388 :             A = pgf_product_list(i)%ra
    1028     7599388 :             B = pgf_product_list(i)%rb
    1029     7599388 :             C = pgf_product_list(i)%rc
    1030     7599388 :             D = pgf_product_list(i)%rd
    1031     1899847 :             Rho = pgf_product_list(i)%Rho
    1032     1899847 :             RhoInv = pgf_product_list(i)%RhoInv
    1033     7599388 :             P = pgf_product_list(i)%P
    1034     7599388 :             Q = pgf_product_list(i)%Q
    1035     7599388 :             W = pgf_product_list(i)%W
    1036     7599388 :             AB = pgf_product_list(i)%AB
    1037     7599388 :             CD = pgf_product_list(i)%CD
    1038     1899847 :             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     1899847 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1044     1899847 :             CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize)
    1045     7599388 :             DO k = 4, 6
    1046   182856172 :                DO j = 1, mysize
    1047             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1048             :                                                work_forces(j, k + 3) + &
    1049   180956325 :                                                work_forces(j, k + 6))
    1050             :                END DO
    1051             :             END DO
    1052    24698011 :             DO k = 1, 12
    1053   725725147 :                DO j = 1, mysize
    1054   723825300 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1055             :                END DO
    1056             :             END DO
    1057     1899847 :             neris = neris + 12*mysize
    1058     2450212 :             IF (use_virial) THEN
    1059      325994 :                CALL real_to_scaled(scoord(1:3), D, cell)
    1060      325994 :                CALL real_to_scaled(scoord(4:6), C, cell)
    1061      325994 :                CALL real_to_scaled(scoord(7:9), A, cell)
    1062      325994 :                CALL real_to_scaled(scoord(10:12), B, cell)
    1063     4237922 :                DO k = 1, 12
    1064   119774270 :                   DO j = 1, mysize
    1065   466057320 :                      DO m = 1, 3
    1066   462145392 :                         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     7154745 :          DO n = 1, 12
    1074             :             tmp_max = 0.0_dp
    1075   226299924 :             DO i = 1, mysize
    1076   226299924 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1077             :             END DO
    1078     6604380 :             tmp_max = tmp_max*max_contraction
    1079     6604380 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1080             : 
    1081    27487617 :             DO l = 1, ncod
    1082    20332872 :                p1 = (l - 1)*ncoc
    1083    47554668 :                DO k = 1, ncoc
    1084    20617416 :                   p2 = (p1 + k - 1)*ncoa
    1085   123494904 :                   DO i = 1, ncoa
    1086    82544616 :                      p3 = (p2 + i - 1)*ncob
    1087   322857576 :                      DO j = 1, ncob
    1088   219695544 :                         p4 = p3 + j
    1089   302240160 :                         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      550365 :          IF (use_virial) THEN
    1096      640822 :             DO n = 1, 12
    1097             :                tmp_max_virial = 0.0_dp
    1098    21936864 :                DO i = 1, mysize
    1099             :                   tmp_max_virial = MAX(tmp_max_virial, &
    1100    21936864 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1101             :                END DO
    1102      591528 :                tmp_max_virial = tmp_max_virial*max_contraction
    1103      591528 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1104             : 
    1105     2471998 :                DO l = 1, ncod
    1106     1831176 :                   p1 = (l - 1)*ncoc
    1107     4290744 :                   DO k = 1, ncoc
    1108     1868040 :                      p2 = (p1 + k - 1)*ncoa
    1109    10871928 :                      DO i = 1, ncoa
    1110     7172712 :                         p3 = (p2 + i - 1)*ncob
    1111    30386088 :                         DO j = 1, ncob
    1112    21345336 :                            p4 = p3 + j
    1113    92554056 :                            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      312733 :          DO i = 1, nproducts
    1122      903216 :             A = pgf_product_list(i)%ra
    1123      903216 :             B = pgf_product_list(i)%rb
    1124      903216 :             C = pgf_product_list(i)%rc
    1125      903216 :             D = pgf_product_list(i)%rd
    1126      225804 :             Rho = pgf_product_list(i)%Rho
    1127      225804 :             RhoInv = pgf_product_list(i)%RhoInv
    1128      903216 :             P = pgf_product_list(i)%P
    1129      903216 :             Q = pgf_product_list(i)%Q
    1130      903216 :             W = pgf_product_list(i)%W
    1131      903216 :             AB = pgf_product_list(i)%AB
    1132      903216 :             CD = pgf_product_list(i)%CD
    1133      225804 :             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      225804 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1139      225804 :             CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize)
    1140      903216 :             DO k = 4, 6
    1141    28733196 :                DO j = 1, mysize
    1142             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1143             :                                                work_forces(j, k + 3) + &
    1144    28507392 :                                                work_forces(j, k + 6))
    1145             :                END DO
    1146             :             END DO
    1147     2935452 :             DO k = 1, 12
    1148   114255372 :                DO j = 1, mysize
    1149   114029568 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1150             :                END DO
    1151             :             END DO
    1152      225804 :             neris = neris + 12*mysize
    1153      312733 :             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     1130077 :          DO n = 1, 12
    1169             :             tmp_max = 0.0_dp
    1170    44392404 :             DO i = 1, mysize
    1171    44392404 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1172             :             END DO
    1173     1043148 :             tmp_max = tmp_max*max_contraction
    1174     1043148 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1175             : 
    1176     4574665 :             DO l = 1, ncod
    1177     3444588 :                p1 = (l - 1)*ncoc
    1178     7932324 :                DO k = 1, ncoc
    1179     3444588 :                   p2 = (p1 + k - 1)*ncob
    1180    27556704 :                   DO j = 1, ncob
    1181    20667528 :                      p3 = (p2 + j - 1)*ncoa
    1182    67461372 :                      DO i = 1, ncoa
    1183    43349256 :                         p4 = p3 + i
    1184    64016784 :                         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    31868381 :          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    31781452 :       IF (.NOT. use_virial) THEN
    1218    30715998 :          IF (tmp_max_all < eps_schwarz) RETURN
    1219             :       END IF
    1220             : 
    1221    26437303 :       IF (tmp_max_all >= eps_schwarz) THEN
    1222   341522974 :          DO i = 1, 12
    1223 16305232248 :             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   315251976 :                           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 16331503246 :                           s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
    1241             :          END DO
    1242             :       END IF
    1243             : 
    1244    26437303 :       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    81977862 :    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    81977862 :       m_max = n_a + n_b + n_c + n_d
    1313    81977862 :       mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
    1314   163955724 :       a_mysize = mysize
    1315             : 
    1316    81977862 :       IF (m_max /= 0) THEN
    1317    51235482 :          perm_case = 1
    1318    51235482 :          IF (n_a < n_b) THEN
    1319    20829432 :             perm_case = perm_case + 1
    1320             :          END IF
    1321    51235482 :          IF (n_c < n_d) THEN
    1322    20829432 :             perm_case = perm_case + 2
    1323             :          END IF
    1324    51235482 :          IF (n_a + n_b > n_c + n_d) THEN
    1325           0 :             perm_case = perm_case + 4
    1326             :          END IF
    1327             : 
    1328    30406050 :          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    30406050 :                                            potential_parameter, libint, R1, R2)
    1332             : 
    1333    30406050 :             CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1334  2625005958 :             DO i = 1, mysize
    1335  2625005958 :                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    20829432 :                                            potential_parameter, libint, R1, R2)
    1356             : 
    1357    20829432 :             CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1358   944781876 :             DO i = 1, mysize
    1359   944781876 :                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    51235482 :             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    30742380 :                                         potential_parameter, libint, R1, R2)
    1397    30742380 :          max_val = ABS(get_ssss_f_val(libint))
    1398             :       END IF
    1399             : 
    1400    81977862 :    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   151366737 :    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   151366737 :                            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   151366737 :                            sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
    1463   151366737 :                            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   151366737 :       m_max = n_a + n_b + n_c + n_d
    1496   151366737 :       mysize = ncoa*ncob*ncoc*ncod
    1497   302733474 :       a_mysize = mysize
    1498             : 
    1499  3285595663 :       work = 0.0_dp
    1500   151366737 :       IF (m_max /= 0) THEN
    1501   119750279 :          perm_case = 1
    1502   119750279 :          IF (n_a < n_b) THEN
    1503    28710554 :             perm_case = perm_case + 1
    1504             :          END IF
    1505   119750279 :          IF (n_c < n_d) THEN
    1506    34926048 :             perm_case = perm_case + 2
    1507             :          END IF
    1508   119750279 :          IF (n_a + n_b > n_c + n_d) THEN
    1509    41910009 :             perm_case = perm_case + 4
    1510             :          END IF
    1511             :          SELECT CASE (perm_case)
    1512             :          CASE (1)
    1513   139029682 :             DO i = 1, nproducts
    1514   391129252 :                A = pgf_product_list(i)%ra
    1515   391129252 :                B = pgf_product_list(i)%rb
    1516   391129252 :                C = pgf_product_list(i)%rc
    1517   391129252 :                D = pgf_product_list(i)%rd
    1518    97782313 :                Rho = pgf_product_list(i)%Rho
    1519    97782313 :                RhoInv = pgf_product_list(i)%RhoInv
    1520   391129252 :                P = pgf_product_list(i)%P
    1521   391129252 :                Q = pgf_product_list(i)%Q
    1522   391129252 :                W = pgf_product_list(i)%W
    1523   391129252 :                AB = pgf_product_list(i)%AB
    1524   391129252 :                CD = pgf_product_list(i)%CD
    1525    97782313 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1526   392257310 :                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    97782313 :                                        P, Q, W, m_max, F)
    1530    97782313 :                CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1531  1947477832 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1532   139029682 :                neris = neris + mysize
    1533             :             END DO
    1534  1029316301 :             DO i = 1, mysize
    1535  1029316301 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1536             :             END DO
    1537    41247369 :             tmp_max = tmp_max*max_contraction
    1538    41247369 :             IF (tmp_max < eps_schwarz) THEN
    1539    35159896 :                RETURN
    1540             :             END IF
    1541             : 
    1542    97186102 :             DO i = 1, ncoa
    1543    64105235 :                p1 = (i - 1)*ncob
    1544   179275665 :                DO j = 1, ncob
    1545    82089563 :                   p2 = (p1 + j - 1)*ncoc
    1546   477303694 :                   DO k = 1, ncoc
    1547   331108896 :                      p3 = (p2 + k - 1)*ncod
    1548  1271808385 :                      DO l = 1, ncod
    1549   858609926 :                         p4 = p3 + l
    1550  1189718822 :                         work2(i, j, k, l) = work(p4)
    1551             :                      END DO
    1552             :                   END DO
    1553             :                END DO
    1554             :             END DO
    1555             :          CASE (2)
    1556    26092912 :             DO i = 1, nproducts
    1557    70312088 :                A = pgf_product_list(i)%ra
    1558    70312088 :                B = pgf_product_list(i)%rb
    1559    70312088 :                C = pgf_product_list(i)%rc
    1560    70312088 :                D = pgf_product_list(i)%rd
    1561    17578022 :                Rho = pgf_product_list(i)%Rho
    1562    17578022 :                RhoInv = pgf_product_list(i)%RhoInv
    1563    70312088 :                P = pgf_product_list(i)%P
    1564    70312088 :                Q = pgf_product_list(i)%Q
    1565    70312088 :                W = pgf_product_list(i)%W
    1566    70312088 :                AB = pgf_product_list(i)%AB
    1567    70312088 :                CD = pgf_product_list(i)%CD
    1568    17578022 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1569    84513253 :                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    17578022 :                                        P, Q, W, m_max, F)
    1574             : 
    1575    17578022 :                CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
    1576   596981568 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1577    26092912 :                neris = neris + mysize
    1578             :             END DO
    1579   383104619 :             DO i = 1, mysize
    1580   383104619 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1581             :             END DO
    1582     8514890 :             tmp_max = tmp_max*max_contraction
    1583     8514890 :             IF (tmp_max < eps_schwarz) THEN
    1584             :                RETURN
    1585             :             END IF
    1586             : 
    1587    28413696 :             DO j = 1, ncob
    1588    22205729 :                p1 = (j - 1)*ncoa
    1589    53507883 :                DO i = 1, ncoa
    1590    25094187 :                   p2 = (p1 + i - 1)*ncoc
    1591   156130848 :                   DO k = 1, ncoc
    1592   108830932 :                      p3 = (p2 + k - 1)*ncod
    1593   436319437 :                      DO l = 1, ncod
    1594   302394318 :                         p4 = p3 + l
    1595   411225250 :                         work2(i, j, k, l) = work(p4)
    1596             :                      END DO
    1597             :                   END DO
    1598             :                END DO
    1599             :             END DO
    1600             :          CASE (3)
    1601    65705238 :             DO i = 1, nproducts
    1602   175232976 :                A = pgf_product_list(i)%ra
    1603   175232976 :                B = pgf_product_list(i)%rb
    1604   175232976 :                C = pgf_product_list(i)%rc
    1605   175232976 :                D = pgf_product_list(i)%rd
    1606    43808244 :                Rho = pgf_product_list(i)%Rho
    1607    43808244 :                RhoInv = pgf_product_list(i)%RhoInv
    1608   175232976 :                P = pgf_product_list(i)%P
    1609   175232976 :                Q = pgf_product_list(i)%Q
    1610   175232976 :                W = pgf_product_list(i)%W
    1611   175232976 :                AB = pgf_product_list(i)%AB
    1612   175232976 :                CD = pgf_product_list(i)%CD
    1613    43808244 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1614   166835645 :                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    43808244 :                                        P, Q, W, m_max, F)
    1619             : 
    1620    43808244 :                CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
    1621   672469242 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1622    65705238 :                neris = neris + mysize
    1623             :             END DO
    1624   414137021 :             DO i = 1, mysize
    1625   414137021 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1626             :             END DO
    1627    21896994 :             tmp_max = tmp_max*max_contraction
    1628    21896994 :             IF (tmp_max < eps_schwarz) THEN
    1629             :                RETURN
    1630             :             END IF
    1631             : 
    1632    48274941 :             DO i = 1, ncoa
    1633    32340860 :                p1 = (i - 1)*ncob
    1634    87765065 :                DO j = 1, ncob
    1635    39490124 :                   p2 = (p1 + j - 1)*ncod
    1636   256242431 :                   DO l = 1, ncod
    1637   184411447 :                      p3 = (p2 + l - 1)*ncoc
    1638   542127500 :                      DO k = 1, ncoc
    1639   318225929 :                         p4 = p3 + k
    1640   502637376 :                         work2(i, j, k, l) = work(p4)
    1641             :                      END DO
    1642             :                   END DO
    1643             :                END DO
    1644             :             END DO
    1645             :          CASE (4)
    1646    17046671 :             DO i = 1, nproducts
    1647    43462616 :                A = pgf_product_list(i)%ra
    1648    43462616 :                B = pgf_product_list(i)%rb
    1649    43462616 :                C = pgf_product_list(i)%rc
    1650    43462616 :                D = pgf_product_list(i)%rd
    1651    10865654 :                Rho = pgf_product_list(i)%Rho
    1652    10865654 :                RhoInv = pgf_product_list(i)%RhoInv
    1653    43462616 :                P = pgf_product_list(i)%P
    1654    43462616 :                Q = pgf_product_list(i)%Q
    1655    43462616 :                W = pgf_product_list(i)%W
    1656    43462616 :                AB = pgf_product_list(i)%AB
    1657    43462616 :                CD = pgf_product_list(i)%CD
    1658    10865654 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1659    49940156 :                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    10865654 :                                        P, Q, W, m_max, F)
    1664             : 
    1665    10865654 :                CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1666   311361772 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1667    17046671 :                neris = neris + mysize
    1668             :             END DO
    1669   227209890 :             DO i = 1, mysize
    1670   227209890 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1671             :             END DO
    1672     6181017 :             tmp_max = tmp_max*max_contraction
    1673     6181017 :             IF (tmp_max < eps_schwarz) THEN
    1674             :                RETURN
    1675             :             END IF
    1676             : 
    1677    21488366 :             DO j = 1, ncob
    1678    16764631 :                p1 = (j - 1)*ncoa
    1679    41016451 :                DO i = 1, ncoa
    1680    19528085 :                   p2 = (p1 + i - 1)*ncod
    1681   130214764 :                   DO l = 1, ncod
    1682    93922048 :                      p3 = (p2 + l - 1)*ncoc
    1683   296984765 :                      DO k = 1, ncoc
    1684   183534632 :                         p4 = p3 + k
    1685   277456680 :                         work2(i, j, k, l) = work(p4)
    1686             :                      END DO
    1687             :                   END DO
    1688             :                END DO
    1689             :             END DO
    1690             :          CASE (5)
    1691    69675830 :             DO i = 1, nproducts
    1692   186466120 :                A = pgf_product_list(i)%ra
    1693   186466120 :                B = pgf_product_list(i)%rb
    1694   186466120 :                C = pgf_product_list(i)%rc
    1695   186466120 :                D = pgf_product_list(i)%rd
    1696    46616530 :                Rho = pgf_product_list(i)%Rho
    1697    46616530 :                RhoInv = pgf_product_list(i)%RhoInv
    1698   186466120 :                P = pgf_product_list(i)%P
    1699   186466120 :                Q = pgf_product_list(i)%Q
    1700   186466120 :                W = pgf_product_list(i)%W
    1701   186466120 :                AB = pgf_product_list(i)%AB
    1702   186466120 :                CD = pgf_product_list(i)%CD
    1703    46616530 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1704   180314205 :                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    46616530 :                                        Q, P, W, m_max, F)
    1709             : 
    1710    46616530 :                CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
    1711   840993142 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1712    69675830 :                neris = neris + mysize
    1713             :             END DO
    1714   525787087 :             DO i = 1, mysize
    1715   525787087 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1716             :             END DO
    1717    23059300 :             tmp_max = tmp_max*max_contraction
    1718    23059300 :             IF (tmp_max < eps_schwarz) THEN
    1719             :                RETURN
    1720             :             END IF
    1721             : 
    1722    42270484 :             DO k = 1, ncoc
    1723    25257179 :                p1 = (k - 1)*ncod
    1724    71636653 :                DO l = 1, ncod
    1725    29366169 :                   p2 = (p1 + l - 1)*ncoa
    1726   190147419 :                   DO i = 1, ncoa
    1727   135524071 :                      p3 = (p2 + i - 1)*ncob
    1728   572552853 :                      DO j = 1, ncob
    1729   407662613 :                         p4 = p3 + j
    1730   543186684 :                         work2(i, j, k, l) = work(p4)
    1731             :                      END DO
    1732             :                   END DO
    1733             :                END DO
    1734             :             END DO
    1735             :          CASE (6)
    1736    31867442 :             DO i = 1, nproducts
    1737    79459080 :                A = pgf_product_list(i)%ra
    1738    79459080 :                B = pgf_product_list(i)%rb
    1739    79459080 :                C = pgf_product_list(i)%rc
    1740    79459080 :                D = pgf_product_list(i)%rd
    1741    19864770 :                Rho = pgf_product_list(i)%Rho
    1742    19864770 :                RhoInv = pgf_product_list(i)%RhoInv
    1743    79459080 :                P = pgf_product_list(i)%P
    1744    79459080 :                Q = pgf_product_list(i)%Q
    1745    79459080 :                W = pgf_product_list(i)%W
    1746    79459080 :                AB = pgf_product_list(i)%AB
    1747    79459080 :                CD = pgf_product_list(i)%CD
    1748    19864770 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1749    72568375 :                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    19864770 :                                        Q, P, W, m_max, F)
    1754             : 
    1755    19864770 :                CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
    1756   295492901 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1757    31867442 :                neris = neris + mysize
    1758             :             END DO
    1759   204666543 :             DO i = 1, mysize
    1760   204666543 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1761             :             END DO
    1762    12002672 :             tmp_max = tmp_max*max_contraction
    1763    12002672 :             IF (tmp_max < eps_schwarz) THEN
    1764             :                RETURN
    1765             :             END IF
    1766             : 
    1767    18920230 :             DO k = 1, ncoc
    1768    11140197 :                p1 = (k - 1)*ncod
    1769    32079583 :                DO l = 1, ncod
    1770    13159353 :                   p2 = (p1 + l - 1)*ncob
    1771    91400325 :                   DO j = 1, ncob
    1772    67100775 :                      p3 = (p2 + j - 1)*ncoa
    1773   223748375 :                      DO i = 1, ncoa
    1774   143488247 :                         p4 = p3 + i
    1775   210589022 :                         work2(i, j, k, l) = work(p4)
    1776             :                      END DO
    1777             :                   END DO
    1778             :                END DO
    1779             :             END DO
    1780             :          CASE (7)
    1781    13674074 :             DO i = 1, nproducts
    1782    35352048 :                A = pgf_product_list(i)%ra
    1783    35352048 :                B = pgf_product_list(i)%rb
    1784    35352048 :                C = pgf_product_list(i)%rc
    1785    35352048 :                D = pgf_product_list(i)%rd
    1786     8838012 :                Rho = pgf_product_list(i)%Rho
    1787     8838012 :                RhoInv = pgf_product_list(i)%RhoInv
    1788    35352048 :                P = pgf_product_list(i)%P
    1789    35352048 :                Q = pgf_product_list(i)%Q
    1790    35352048 :                W = pgf_product_list(i)%W
    1791    35352048 :                AB = pgf_product_list(i)%AB
    1792    35352048 :                CD = pgf_product_list(i)%CD
    1793     8838012 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1794    47881408 :                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     8838012 :                                        Q, P, W, m_max, F)
    1799             : 
    1800     8838012 :                CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
    1801   465371597 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1802    13674074 :                neris = neris + mysize
    1803             :             END DO
    1804   331973799 :             DO i = 1, mysize
    1805   331973799 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1806             :             END DO
    1807     4836062 :             tmp_max = tmp_max*max_contraction
    1808     4836062 :             IF (tmp_max < eps_schwarz) THEN
    1809             :                RETURN
    1810             :             END IF
    1811             : 
    1812    15232201 :             DO l = 1, ncod
    1813    11837116 :                p1 = (l - 1)*ncoc
    1814    28730373 :                DO k = 1, ncoc
    1815    13498172 :                   p2 = (p1 + k - 1)*ncoa
    1816    96255558 :                   DO i = 1, ncoa
    1817    70920270 :                      p3 = (p2 + i - 1)*ncob
    1818   349371132 :                      DO j = 1, ncob
    1819   264952690 :                         p4 = p3 + j
    1820   335872960 :                         work2(i, j, k, l) = work(p4)
    1821             :                      END DO
    1822             :                   END DO
    1823             :                END DO
    1824             :             END DO
    1825             :          CASE (8)
    1826     4645528 :             DO i = 1, nproducts
    1827    10534212 :                A = pgf_product_list(i)%ra
    1828    10534212 :                B = pgf_product_list(i)%rb
    1829    10534212 :                C = pgf_product_list(i)%rc
    1830    10534212 :                D = pgf_product_list(i)%rd
    1831     2633553 :                Rho = pgf_product_list(i)%Rho
    1832     2633553 :                RhoInv = pgf_product_list(i)%RhoInv
    1833    10534212 :                P = pgf_product_list(i)%P
    1834    10534212 :                Q = pgf_product_list(i)%Q
    1835    10534212 :                W = pgf_product_list(i)%W
    1836    10534212 :                AB = pgf_product_list(i)%AB
    1837    10534212 :                CD = pgf_product_list(i)%CD
    1838     2633553 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1839    15008081 :                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     2633553 :                                        Q, P, W, m_max, F)
    1844             : 
    1845     2633553 :                CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
    1846   130897443 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1847     4645528 :                neris = neris + mysize
    1848             :             END DO
    1849   106167487 :             DO i = 1, mysize
    1850   106167487 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1851             :             END DO
    1852     2011975 :             tmp_max = tmp_max*max_contraction
    1853     2011975 :             IF (tmp_max < eps_schwarz) THEN
    1854             :                RETURN
    1855             :             END IF
    1856             : 
    1857   126837756 :             DO l = 1, ncod
    1858     5538850 :                p1 = (l - 1)*ncoc
    1859    12686647 :                DO k = 1, ncoc
    1860     5599170 :                   p2 = (p1 + k - 1)*ncob
    1861    45802924 :                   DO j = 1, ncob
    1862    34664904 :                      p3 = (p2 + j - 1)*ncoa
    1863   122210892 :                      DO i = 1, ncoa
    1864    81946818 :                         p4 = p3 + i
    1865   116611722 :                         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   100496747 :          DO i = 1, nproducts
    1873   275521156 :             A = pgf_product_list(i)%ra
    1874   275521156 :             B = pgf_product_list(i)%rb
    1875   275521156 :             C = pgf_product_list(i)%rc
    1876   275521156 :             D = pgf_product_list(i)%rd
    1877    68880289 :             Rho = pgf_product_list(i)%Rho
    1878    68880289 :             RhoInv = pgf_product_list(i)%RhoInv
    1879   275521156 :             P = pgf_product_list(i)%P
    1880   275521156 :             Q = pgf_product_list(i)%Q
    1881   275521156 :             W = pgf_product_list(i)%W
    1882    68880289 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1883   137760578 :             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    68880289 :                                     P, Q, W, m_max, F)
    1888    68880289 :             work(1) = work(1) + F(1)
    1889   100496747 :             neris = neris + mysize
    1890             :          END DO
    1891    31616458 :          work2(1, 1, 1, 1) = work(1)
    1892    31616458 :          tmp_max = max_contraction*ABS(work(1))
    1893    31616458 :          IF (tmp_max < eps_schwarz) RETURN
    1894             :       END IF
    1895             : 
    1896   116206841 :       IF (tmp_max < eps_schwarz) RETURN
    1897  7718217567 :       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   116206841 :                     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  7718217567 :                     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   316867387 :    SUBROUTINE build_quartet_data(libint, A, B, C, D, &
    1940             :                                  ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1941   316867387 :                                  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   316867387 :       CALL cp_libint_set_params_eri(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, P, Q, W, m_max, F)
    1951   316867387 :    END SUBROUTINE build_quartet_data
    1952             : 
    1953             : END MODULE hfx_libint_interface

Generated by: LCOV version 1.15