LCOV - code coverage report
Current view: top level - src - hfx_libint_interface.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:ccc2433) Lines: 951 983 96.7 %
Date: 2024-04-25 07:09:54 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    82215414 :    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    82215414 :       Zeta = Zeta_A + Zeta_B
      98    82215414 :       ZetaInv = 1.0_dp/Zeta
      99    82215414 :       Eta = Zeta_C + Zeta_D
     100    82215414 :       EtaInv = 1.0_dp/Eta
     101    82215414 :       ZetapEtaInv = Zeta + Eta
     102    82215414 :       ZetapEtaInv = 1.0_dp/ZetapEtaInv
     103    82215414 :       Rho = Zeta*Eta*ZetapEtaInv
     104    82215414 :       RhoInv = 1.0_dp/Rho
     105             : 
     106   328861656 :       DO i = 1, 3
     107   246646242 :          P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv
     108   246646242 :          Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv
     109   246646242 :          AB(i) = A(i) - B(i)
     110   246646242 :          CD(i) = C(i) - D(i)
     111   246646242 :          PQ(i) = P(i) - Q(i)
     112   328861656 :          W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv
     113             :       END DO
     114             : 
     115   328861656 :       AB2 = DOT_PRODUCT(AB, AB)
     116   328861656 :       CD2 = DOT_PRODUCT(CD, CD)
     117   328861656 :       PQ2 = DOT_PRODUCT(PQ, PQ)
     118             : 
     119    82215414 :       S1234 = EXP((-Zeta_A*Zeta_B*ZetaInv*AB2) + (-Zeta_C*Zeta_D*EtaInv*CD2))
     120    82215414 :       T = Rho*PQ2
     121             : 
     122    94745676 :       SELECT CASE (potential_parameter%potential_type)
     123             :       CASE (do_potential_truncated)
     124    12530262 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     125    12530262 :          R1 = R11
     126    12530262 :          R2 = R22
     127    12530262 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     128           0 :             RETURN
     129             :          END IF
     130    12530262 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     131    12530262 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     132    12530262 :          factor = 2.0_dp*Pi*RhoInv
     133             :       CASE (do_potential_coulomb)
     134    62611920 :          CALL fgamma(m_max, T, F(1))
     135    62611920 :          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    95927982 :          IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
     256             :       END SELECT
     257             : 
     258    82215414 :       tmp = (Pi*ZetapEtaInv)**3
     259    82215414 :       factor = factor*S1234*SQRT(tmp)
     260             : 
     261   317150100 :       DO i = 1, m_max + 1
     262   317150100 :          F(i) = F(i)*factor
     263             :       END DO
     264             : 
     265    82215414 :       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    93292540 :    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    93292540 :                                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    93292540 :                                           ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F)
     313             : 
     314    93292540 :    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    31773339 :    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    31773339 :                                  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    31773339 :                                  sphi_a, sphi_b, sphi_c, sphi_d, &
     388    31773339 :                                  work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
     389    31773339 :                                  use_virial, work_virial, work2_virial, primitives_tmp_virial, &
     390    31773339 :                                  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    31773339 :       m_max = n_a + n_b + n_c + n_d
     433    31773339 :       m_max = m_max + 1
     434    31773339 :       mysize = ncoa*ncob*ncoc*ncod
     435    63546678 :       a_mysize = mysize
     436             : 
     437  4151355207 :       work = 0.0_dp
     438  7998928155 :       work2 = 0.0_dp
     439             : 
     440    31773339 :       IF (use_virial) THEN
     441   850515572 :          work_virial = 0.0_dp
     442  1465421096 :          work2_virial = 0.0_dp
     443             :       END IF
     444             : 
     445    31773339 :       perm_case = 1
     446    31773339 :       IF (n_a < n_b) THEN
     447     2916283 :          perm_case = perm_case + 1
     448             :       END IF
     449    31773339 :       IF (n_c < n_d) THEN
     450     5030786 :          perm_case = perm_case + 2
     451             :       END IF
     452    31773339 :       IF (n_a + n_b > n_c + n_d) THEN
     453     5722142 :          perm_case = perm_case + 4
     454             :       END IF
     455             : 
     456             :       SELECT CASE (perm_case)
     457             :       CASE (1)
     458    77826291 :          DO i = 1, nproducts
     459   229320344 :             A = pgf_product_list(i)%ra
     460   229320344 :             B = pgf_product_list(i)%rb
     461   229320344 :             C = pgf_product_list(i)%rc
     462   229320344 :             D = pgf_product_list(i)%rd
     463    57330086 :             Rho = pgf_product_list(i)%Rho
     464    57330086 :             RhoInv = pgf_product_list(i)%RhoInv
     465   229320344 :             P = pgf_product_list(i)%P
     466   229320344 :             Q = pgf_product_list(i)%Q
     467   229320344 :             W = pgf_product_list(i)%W
     468   229320344 :             AB = pgf_product_list(i)%AB
     469   229320344 :             CD = pgf_product_list(i)%CD
     470    57330086 :             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    57330086 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     476    57330086 :             CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize)
     477   229320344 :             DO k = 4, 6
     478  1742843006 :                DO j = 1, mysize
     479             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     480             :                                                work_forces(j, k + 3) + &
     481  1685512920 :                                                work_forces(j, k + 6))
     482             :                END DO
     483             :             END DO
     484   745291118 :             DO k = 1, 12
     485  6799381766 :                DO j = 1, mysize
     486  6742051680 :                   work(j, k) = work(j, k) + work_forces(j, k)
     487             :                END DO
     488             :             END DO
     489    57330086 :             neris = neris + 12*mysize
     490    77826291 :             IF (use_virial) THEN
     491    10719799 :                CALL real_to_scaled(scoord(1:3), A, cell)
     492    10719799 :                CALL real_to_scaled(scoord(4:6), B, cell)
     493    10719799 :                CALL real_to_scaled(scoord(7:9), C, cell)
     494    10719799 :                CALL real_to_scaled(scoord(10:12), D, cell)
     495   139357387 :                DO k = 1, 12
     496  1164832423 :                   DO j = 1, mysize
     497  4230537732 :                      DO m = 1, 3
     498  4101900144 :                         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   266450665 :          DO n = 1, 12
     506             :             tmp_max = 0.0_dp
     507  2078889192 :             DO i = 1, mysize
     508  2078889192 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     509             :             END DO
     510   245954460 :             tmp_max = tmp_max*max_contraction
     511   245954460 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     512             : 
     513   605316769 :             DO i = 1, ncoa
     514   338866104 :                p1 = (i - 1)*ncob
     515   957639708 :                DO j = 1, ncob
     516   372819144 :                   p2 = (p1 + j - 1)*ncoc
     517  1710151428 :                   DO k = 1, ncoc
     518   998466180 :                      p3 = (p2 + k - 1)*ncod
     519  3204220056 :                      DO l = 1, ncod
     520  1832934732 :                         p4 = p3 + l
     521  2831400912 :                         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    20496205 :          IF (use_virial) THEN
     528     6234475 :             DO n = 1, 12
     529             :                tmp_max_virial = 0.0_dp
     530   112743228 :                DO i = 1, mysize
     531             :                   tmp_max_virial = MAX(tmp_max_virial, &
     532   112743228 :                                        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    16529887 :                DO i = 1, ncoa
     538    10295412 :                   p1 = (i - 1)*ncob
     539    29968620 :                   DO j = 1, ncob
     540    13918308 :                      p2 = (p1 + j - 1)*ncoc
     541    69611376 :                      DO k = 1, ncoc
     542    45397656 :                         p3 = (p2 + k - 1)*ncod
     543   166304292 :                         DO l = 1, ncod
     544   106988328 :                            p4 = p3 + l
     545   473350968 :                            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     5364238 :          DO i = 1, nproducts
     554    16811844 :             A = pgf_product_list(i)%ra
     555    16811844 :             B = pgf_product_list(i)%rb
     556    16811844 :             C = pgf_product_list(i)%rc
     557    16811844 :             D = pgf_product_list(i)%rd
     558     4202961 :             Rho = pgf_product_list(i)%Rho
     559     4202961 :             RhoInv = pgf_product_list(i)%RhoInv
     560    16811844 :             P = pgf_product_list(i)%P
     561    16811844 :             Q = pgf_product_list(i)%Q
     562    16811844 :             W = pgf_product_list(i)%W
     563    16811844 :             AB = pgf_product_list(i)%AB
     564    16811844 :             CD = pgf_product_list(i)%CD
     565     4202961 :             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     4202961 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     571     4202961 :             CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize)
     572    16811844 :             DO k = 4, 6
     573   298949343 :                DO j = 1, mysize
     574             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     575             :                                                work_forces(j, k + 3) + &
     576   294746382 :                                                work_forces(j, k + 6))
     577             :                END DO
     578             :             END DO
     579    54638493 :             DO k = 1, 12
     580  1183188489 :                DO j = 1, mysize
     581  1178985528 :                   work(j, k) = work(j, k) + work_forces(j, k)
     582             :                END DO
     583             :             END DO
     584     4202961 :             neris = neris + 12*mysize
     585     5364238 :             IF (use_virial) THEN
     586      673397 :                CALL real_to_scaled(scoord(1:3), B, cell)
     587      673397 :                CALL real_to_scaled(scoord(4:6), A, cell)
     588      673397 :                CALL real_to_scaled(scoord(7:9), C, cell)
     589      673397 :                CALL real_to_scaled(scoord(10:12), D, cell)
     590     8754161 :                DO k = 1, 12
     591   180257837 :                   DO j = 1, mysize
     592   694095468 :                      DO m = 1, 3
     593   686014704 :                         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    15096601 :          DO n = 1, 12
     601             :             tmp_max = 0.0_dp
     602   361339788 :             DO i = 1, mysize
     603   361339788 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     604             :             END DO
     605    13935324 :             tmp_max = tmp_max*max_contraction
     606    13935324 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     607             : 
     608    59300461 :             DO j = 1, ncob
     609    44203860 :                p1 = (j - 1)*ncoa
     610   103767780 :                DO i = 1, ncoa
     611    45628596 :                   p2 = (p1 + i - 1)*ncoc
     612   253749456 :                   DO k = 1, ncoc
     613   163917000 :                      p3 = (p2 + k - 1)*ncod
     614   556950060 :                      DO l = 1, ncod
     615   347404464 :                         p4 = p3 + l
     616   511321464 :                         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     1161277 :          IF (use_virial) THEN
     623     1274234 :             DO n = 1, 12
     624             :                tmp_max_virial = 0.0_dp
     625    34689696 :                DO i = 1, mysize
     626             :                   tmp_max_virial = MAX(tmp_max_virial, &
     627    34689696 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     628             :                END DO
     629     1176216 :                tmp_max_virial = tmp_max_virial*max_contraction
     630     1176216 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     631             : 
     632     5079362 :                DO j = 1, ncob
     633     3805128 :                   p1 = (j - 1)*ncoa
     634     8970792 :                   DO i = 1, ncoa
     635     3989448 :                      p2 = (p1 + i - 1)*ncoc
     636    22275864 :                      DO k = 1, ncoc
     637    14481288 :                         p3 = (p2 + k - 1)*ncod
     638    51984216 :                         DO l = 1, ncod
     639    33513480 :                            p4 = p3 + l
     640   148535208 :                            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    15767506 :          DO i = 1, nproducts
     649    48050772 :             A = pgf_product_list(i)%ra
     650    48050772 :             B = pgf_product_list(i)%rb
     651    48050772 :             C = pgf_product_list(i)%rc
     652    48050772 :             D = pgf_product_list(i)%rd
     653    12012693 :             Rho = pgf_product_list(i)%Rho
     654    12012693 :             RhoInv = pgf_product_list(i)%RhoInv
     655    48050772 :             P = pgf_product_list(i)%P
     656    48050772 :             Q = pgf_product_list(i)%Q
     657    48050772 :             W = pgf_product_list(i)%W
     658    48050772 :             AB = pgf_product_list(i)%AB
     659    48050772 :             CD = pgf_product_list(i)%CD
     660    12012693 :             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    12012693 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     666    12012693 :             CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize)
     667    48050772 :             DO k = 4, 6
     668   420534321 :                DO j = 1, mysize
     669             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     670             :                                                work_forces(j, k + 3) + &
     671   408521628 :                                                work_forces(j, k + 6))
     672             :                END DO
     673             :             END DO
     674   156165009 :             DO k = 1, 12
     675  1646099205 :                DO j = 1, mysize
     676  1634086512 :                   work(j, k) = work(j, k) + work_forces(j, k)
     677             :                END DO
     678             :             END DO
     679    12012693 :             neris = neris + 12*mysize
     680    15767506 :             IF (use_virial) THEN
     681     1825312 :                CALL real_to_scaled(scoord(1:3), A, cell)
     682     1825312 :                CALL real_to_scaled(scoord(4:6), B, cell)
     683     1825312 :                CALL real_to_scaled(scoord(7:9), D, cell)
     684     1825312 :                CALL real_to_scaled(scoord(10:12), C, cell)
     685    23729056 :                DO k = 1, 12
     686   205930600 :                   DO j = 1, mysize
     687   750709920 :                      DO m = 1, 3
     688   728806176 :                         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    48812569 :          DO n = 1, 12
     696             :             tmp_max = 0.0_dp
     697   519006144 :             DO i = 1, mysize
     698   519006144 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     699             :             END DO
     700    45057756 :             tmp_max = tmp_max*max_contraction
     701    45057756 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     702             : 
     703   126903697 :             DO i = 1, ncoa
     704    78091128 :                p1 = (i - 1)*ncob
     705   208416900 :                DO j = 1, ncob
     706    85268016 :                   p2 = (p1 + j - 1)*ncod
     707   491823900 :                   DO l = 1, ncod
     708   328464756 :                      p3 = (p2 + l - 1)*ncoc
     709   887681160 :                      DO k = 1, ncoc
     710   473948388 :                         p4 = p3 + k
     711   802413144 :                         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     3754813 :          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     2722223 :          DO i = 1, nproducts
     744     8333284 :             A = pgf_product_list(i)%ra
     745     8333284 :             B = pgf_product_list(i)%rb
     746     8333284 :             C = pgf_product_list(i)%rc
     747     8333284 :             D = pgf_product_list(i)%rd
     748     2083321 :             Rho = pgf_product_list(i)%Rho
     749     2083321 :             RhoInv = pgf_product_list(i)%RhoInv
     750     8333284 :             P = pgf_product_list(i)%P
     751     8333284 :             Q = pgf_product_list(i)%Q
     752     8333284 :             W = pgf_product_list(i)%W
     753     8333284 :             AB = pgf_product_list(i)%AB
     754     8333284 :             CD = pgf_product_list(i)%CD
     755     2083321 :             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     2083321 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     760     2083321 :             CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize)
     761     8333284 :             DO k = 4, 6
     762   112827172 :                DO j = 1, mysize
     763             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     764             :                                                work_forces(j, k + 3) + &
     765   110743851 :                                                work_forces(j, k + 6))
     766             :                END DO
     767             :             END DO
     768    27083173 :             DO k = 1, 12
     769   445058725 :                DO j = 1, mysize
     770   442975404 :                   work(j, k) = work(j, k) + work_forces(j, k)
     771             :                END DO
     772             :             END DO
     773     2083321 :             neris = neris + 12*mysize
     774     2722223 :             IF (use_virial) THEN
     775      326198 :                CALL real_to_scaled(scoord(1:3), B, cell)
     776      326198 :                CALL real_to_scaled(scoord(4:6), A, cell)
     777      326198 :                CALL real_to_scaled(scoord(7:9), D, cell)
     778      326198 :                CALL real_to_scaled(scoord(10:12), C, cell)
     779     4240574 :                DO k = 1, 12
     780    59380838 :                   DO j = 1, mysize
     781   224475432 :                      DO m = 1, 3
     782   220561056 :                         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     8305726 :          DO n = 1, 12
     790             :             tmp_max = 0.0_dp
     791   153020676 :             DO i = 1, mysize
     792   153020676 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     793             :             END DO
     794     7666824 :             tmp_max = tmp_max*max_contraction
     795     7666824 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     796             : 
     797    32294938 :             DO j = 1, ncob
     798    23989212 :                p1 = (j - 1)*ncoa
     799    56785440 :                DO i = 1, ncoa
     800    25129404 :                   p2 = (p1 + i - 1)*ncod
     801   145693188 :                   DO l = 1, ncod
     802    96574572 :                      p3 = (p2 + l - 1)*ncoc
     803   267057828 :                      DO k = 1, ncoc
     804   145353852 :                         p4 = p3 + k
     805   241928424 :                         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      638902 :          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    16461057 :          DO i = 1, nproducts
     838    49620744 :             A = pgf_product_list(i)%ra
     839    49620744 :             B = pgf_product_list(i)%rb
     840    49620744 :             C = pgf_product_list(i)%rc
     841    49620744 :             D = pgf_product_list(i)%rd
     842    12405186 :             Rho = pgf_product_list(i)%Rho
     843    12405186 :             RhoInv = pgf_product_list(i)%RhoInv
     844    49620744 :             P = pgf_product_list(i)%P
     845    49620744 :             Q = pgf_product_list(i)%Q
     846    49620744 :             W = pgf_product_list(i)%W
     847    49620744 :             AB = pgf_product_list(i)%AB
     848    49620744 :             CD = pgf_product_list(i)%CD
     849    12405186 :             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    12405186 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     854    12405186 :             CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize)
     855    49620744 :             DO k = 4, 6
     856   469790964 :                DO j = 1, mysize
     857             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     858             :                                                work_forces(j, k + 3) + &
     859   457385778 :                                                work_forces(j, k + 6))
     860             :                END DO
     861             :             END DO
     862   161267418 :             DO k = 1, 12
     863  1841948298 :                DO j = 1, mysize
     864  1829543112 :                   work(j, k) = work(j, k) + work_forces(j, k)
     865             :                END DO
     866             :             END DO
     867    12405186 :             neris = neris + 12*mysize
     868    16461057 :             IF (use_virial) THEN
     869     1995233 :                CALL real_to_scaled(scoord(1:3), C, cell)
     870     1995233 :                CALL real_to_scaled(scoord(4:6), D, cell)
     871     1995233 :                CALL real_to_scaled(scoord(7:9), A, cell)
     872     1995233 :                CALL real_to_scaled(scoord(10:12), B, cell)
     873    25938029 :                DO k = 1, 12
     874   267827393 :                   DO j = 1, mysize
     875   991500252 :                      DO m = 1, 3
     876   967557456 :                         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    52726323 :          DO n = 1, 12
     884             :             tmp_max = 0.0_dp
     885   576724500 :             DO i = 1, mysize
     886   576724500 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     887             :             END DO
     888    48670452 :             tmp_max = tmp_max*max_contraction
     889    48670452 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     890             : 
     891   118097943 :             DO k = 1, ncoc
     892    65371620 :                p1 = (k - 1)*ncod
     893   183020964 :                DO l = 1, ncod
     894    68978892 :                   p2 = (p1 + l - 1)*ncoa
     895   392965296 :                   DO i = 1, ncoa
     896   258614784 :                      p3 = (p2 + i - 1)*ncob
     897   855647724 :                      DO j = 1, ncob
     898   528054048 :                         p4 = p3 + j
     899   786668832 :                         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     4055871 :          IF (use_virial) THEN
     906     2252549 :             DO n = 1, 12
     907             :                tmp_max_virial = 0.0_dp
     908    44908836 :                DO i = 1, mysize
     909             :                   tmp_max_virial = MAX(tmp_max_virial, &
     910    44908836 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     911             :                END DO
     912     2079276 :                tmp_max_virial = tmp_max_virial*max_contraction
     913     2079276 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     914             : 
     915     5748953 :                DO k = 1, ncoc
     916     3496404 :                   p1 = (k - 1)*ncod
     917     9492708 :                   DO l = 1, ncod
     918     3917028 :                      p2 = (p1 + l - 1)*ncoa
     919    22983504 :                      DO i = 1, ncoa
     920    15570072 :                         p3 = (p2 + i - 1)*ncob
     921    62316660 :                         DO j = 1, ncob
     922    42829560 :                            p4 = p3 + j
     923   186888312 :                            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     4166028 :          DO i = 1, nproducts
     932    12547312 :             A = pgf_product_list(i)%ra
     933    12547312 :             B = pgf_product_list(i)%rb
     934    12547312 :             C = pgf_product_list(i)%rc
     935    12547312 :             D = pgf_product_list(i)%rd
     936     3136828 :             Rho = pgf_product_list(i)%Rho
     937     3136828 :             RhoInv = pgf_product_list(i)%RhoInv
     938    12547312 :             P = pgf_product_list(i)%P
     939    12547312 :             Q = pgf_product_list(i)%Q
     940    12547312 :             W = pgf_product_list(i)%W
     941    12547312 :             AB = pgf_product_list(i)%AB
     942    12547312 :             CD = pgf_product_list(i)%CD
     943     3136828 :             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     3136828 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     949     3136828 :             CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize)
     950    12547312 :             DO k = 4, 6
     951   118595986 :                DO j = 1, mysize
     952             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     953             :                                                work_forces(j, k + 3) + &
     954   115459158 :                                                work_forces(j, k + 6))
     955             :                END DO
     956             :             END DO
     957    40778764 :             DO k = 1, 12
     958   464973460 :                DO j = 1, mysize
     959   461836632 :                   work(j, k) = work(j, k) + work_forces(j, k)
     960             :                END DO
     961             :             END DO
     962     3136828 :             neris = neris + 12*mysize
     963     4166028 :             IF (use_virial) THEN
     964      380875 :                CALL real_to_scaled(scoord(1:3), C, cell)
     965      380875 :                CALL real_to_scaled(scoord(4:6), D, cell)
     966      380875 :                CALL real_to_scaled(scoord(7:9), B, cell)
     967      380875 :                CALL real_to_scaled(scoord(10:12), A, cell)
     968     4951375 :                DO k = 1, 12
     969    51021583 :                   DO j = 1, mysize
     970   188851332 :                      DO m = 1, 3
     971   184280832 :                         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    13379600 :          DO n = 1, 12
     979             :             tmp_max = 0.0_dp
     980   159997236 :             DO i = 1, mysize
     981   159997236 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     982             :             END DO
     983    12350400 :             tmp_max = tmp_max*max_contraction
     984    12350400 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     985             : 
     986    28889828 :             DO k = 1, ncoc
     987    15510228 :                p1 = (k - 1)*ncod
     988    45647928 :                DO l = 1, ncod
     989    17787300 :                   p2 = (p1 + l - 1)*ncob
     990   111067212 :                   DO j = 1, ncob
     991    77769684 :                      p3 = (p2 + j - 1)*ncoa
     992   243203820 :                      DO i = 1, ncoa
     993   147646836 :                         p4 = p3 + i
     994   225416520 :                         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     1029200 :          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     2445452 :          DO i = 1, nproducts
    1027     7581140 :             A = pgf_product_list(i)%ra
    1028     7581140 :             B = pgf_product_list(i)%rb
    1029     7581140 :             C = pgf_product_list(i)%rc
    1030     7581140 :             D = pgf_product_list(i)%rd
    1031     1895285 :             Rho = pgf_product_list(i)%Rho
    1032     1895285 :             RhoInv = pgf_product_list(i)%RhoInv
    1033     7581140 :             P = pgf_product_list(i)%P
    1034     7581140 :             Q = pgf_product_list(i)%Q
    1035     7581140 :             W = pgf_product_list(i)%W
    1036     7581140 :             AB = pgf_product_list(i)%AB
    1037     7581140 :             CD = pgf_product_list(i)%CD
    1038     1895285 :             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     1895285 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1044     1895285 :             CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize)
    1045     7581140 :             DO k = 4, 6
    1046   182525588 :                DO j = 1, mysize
    1047             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1048             :                                                work_forces(j, k + 3) + &
    1049   180630303 :                                                work_forces(j, k + 6))
    1050             :                END DO
    1051             :             END DO
    1052    24638705 :             DO k = 1, 12
    1053   724416497 :                DO j = 1, mysize
    1054   722521212 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1055             :                END DO
    1056             :             END DO
    1057     1895285 :             neris = neris + 12*mysize
    1058     2445452 :             IF (use_virial) THEN
    1059      326134 :                CALL real_to_scaled(scoord(1:3), D, cell)
    1060      326134 :                CALL real_to_scaled(scoord(4:6), C, cell)
    1061      326134 :                CALL real_to_scaled(scoord(7:9), A, cell)
    1062      326134 :                CALL real_to_scaled(scoord(10:12), B, cell)
    1063     4239742 :                DO k = 1, 12
    1064   119835166 :                   DO j = 1, mysize
    1065   466295304 :                      DO m = 1, 3
    1066   462381696 :                         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     7152171 :          DO n = 1, 12
    1074             :             tmp_max = 0.0_dp
    1075   226223892 :             DO i = 1, mysize
    1076   226223892 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1077             :             END DO
    1078     6602004 :             tmp_max = tmp_max*max_contraction
    1079     6602004 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1080             : 
    1081    27477879 :             DO l = 1, ncod
    1082    20325708 :                p1 = (l - 1)*ncoc
    1083    47537964 :                DO k = 1, ncoc
    1084    20610252 :                   p2 = (p1 + k - 1)*ncoa
    1085   123454872 :                   DO i = 1, ncoa
    1086    82518912 :                      p3 = (p2 + i - 1)*ncob
    1087   322751052 :                      DO j = 1, ncob
    1088   219621888 :                         p4 = p3 + j
    1089   302140800 :                         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      550167 :          IF (use_virial) THEN
    1096      640848 :             DO n = 1, 12
    1097             :                tmp_max_virial = 0.0_dp
    1098    21937536 :                DO i = 1, mysize
    1099             :                   tmp_max_virial = MAX(tmp_max_virial, &
    1100    21937536 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1101             :                END DO
    1102      591552 :                tmp_max_virial = tmp_max_virial*max_contraction
    1103      591552 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1104             : 
    1105     2472096 :                DO l = 1, ncod
    1106     1831248 :                   p1 = (l - 1)*ncoc
    1107     4290912 :                   DO k = 1, ncoc
    1108     1868112 :                      p2 = (p1 + k - 1)*ncoa
    1109    10872288 :                      DO i = 1, ncoa
    1110     7172928 :                         p3 = (p2 + i - 1)*ncob
    1111    30387024 :                         DO j = 1, ncob
    1112    21345984 :                            p4 = p3 + j
    1113    92556864 :                            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      313084 :          DO i = 1, nproducts
    1122      904720 :             A = pgf_product_list(i)%ra
    1123      904720 :             B = pgf_product_list(i)%rb
    1124      904720 :             C = pgf_product_list(i)%rc
    1125      904720 :             D = pgf_product_list(i)%rd
    1126      226180 :             Rho = pgf_product_list(i)%Rho
    1127      226180 :             RhoInv = pgf_product_list(i)%RhoInv
    1128      904720 :             P = pgf_product_list(i)%P
    1129      904720 :             Q = pgf_product_list(i)%Q
    1130      904720 :             W = pgf_product_list(i)%W
    1131      904720 :             AB = pgf_product_list(i)%AB
    1132      904720 :             CD = pgf_product_list(i)%CD
    1133      226180 :             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      226180 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1139      226180 :             CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize)
    1140      904720 :             DO k = 4, 6
    1141    28774876 :                DO j = 1, mysize
    1142             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1143             :                                                work_forces(j, k + 3) + &
    1144    28548696 :                                                work_forces(j, k + 6))
    1145             :                END DO
    1146             :             END DO
    1147     2940340 :             DO k = 1, 12
    1148   114420964 :                DO j = 1, mysize
    1149   114194784 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1150             :                END DO
    1151             :             END DO
    1152      226180 :             neris = neris + 12*mysize
    1153      313084 :             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     1129752 :          DO n = 1, 12
    1169             :             tmp_max = 0.0_dp
    1170    44380440 :             DO i = 1, mysize
    1171    44380440 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1172             :             END DO
    1173     1042848 :             tmp_max = tmp_max*max_contraction
    1174     1042848 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1175             : 
    1176     4573404 :             DO l = 1, ncod
    1177     3443652 :                p1 = (l - 1)*ncoc
    1178     7930152 :                DO k = 1, ncoc
    1179     3443652 :                   p2 = (p1 + k - 1)*ncob
    1180    27549216 :                   DO j = 1, ncob
    1181    20661912 :                      p3 = (p2 + j - 1)*ncoa
    1182    67443156 :                      DO i = 1, ncoa
    1183    43337592 :                         p4 = p3 + i
    1184    63999504 :                         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    31860243 :          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    31773339 :       IF (.NOT. use_virial) THEN
    1218    30707878 :          IF (tmp_max_all < eps_schwarz) RETURN
    1219             :       END IF
    1220             : 
    1221    26428071 :       IF (tmp_max_all >= eps_schwarz) THEN
    1222   341403816 :          DO i = 1, 12
    1223 16297360152 :             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   315141984 :                           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 16323621984 :                           s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
    1241             :          END DO
    1242             :       END IF
    1243             : 
    1244    26428071 :       IF (use_virial .AND. tmp_max_all_virial >= eps_schwarz) THEN
    1245    11813711 :          DO i = 1, 12
    1246    44528603 :             DO m = 1, 3
    1247  4937375340 :                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    32714892 :                              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  4948280304 :                                     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    82215414 :    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    82215414 :       m_max = n_a + n_b + n_c + n_d
    1313    82215414 :       mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
    1314   164430828 :       a_mysize = mysize
    1315             : 
    1316    82215414 :       IF (m_max /= 0) THEN
    1317    51395466 :          perm_case = 1
    1318    51395466 :          IF (n_a < n_b) THEN
    1319    20893668 :             perm_case = perm_case + 1
    1320             :          END IF
    1321    51395466 :          IF (n_c < n_d) THEN
    1322    20893668 :             perm_case = perm_case + 2
    1323             :          END IF
    1324    51395466 :          IF (n_a + n_b > n_c + n_d) THEN
    1325           0 :             perm_case = perm_case + 4
    1326             :          END IF
    1327             : 
    1328    30501798 :          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    30501798 :                                            potential_parameter, libint, R1, R2)
    1332             : 
    1333    30501798 :             CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1334  2631875574 :             DO i = 1, mysize
    1335  2631875574 :                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    20893668 :                                            potential_parameter, libint, R1, R2)
    1356             : 
    1357    20893668 :             CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1358   947594928 :             DO i = 1, mysize
    1359   947594928 :                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    51395466 :             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    30819948 :                                         potential_parameter, libint, R1, R2)
    1397    30819948 :          max_val = ABS(get_ssss_f_val(libint))
    1398             :       END IF
    1399             : 
    1400    82215414 :    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   152267346 :    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   152267346 :                            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   152267346 :                            sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
    1463   152267346 :                            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   152267346 :       m_max = n_a + n_b + n_c + n_d
    1496   152267346 :       mysize = ncoa*ncob*ncoc*ncod
    1497   304534692 :       a_mysize = mysize
    1498             : 
    1499  3300466681 :       work = 0.0_dp
    1500   152267346 :       IF (m_max /= 0) THEN
    1501   120509174 :          perm_case = 1
    1502   120509174 :          IF (n_a < n_b) THEN
    1503    28886688 :             perm_case = perm_case + 1
    1504             :          END IF
    1505   120509174 :          IF (n_c < n_d) THEN
    1506    35115212 :             perm_case = perm_case + 2
    1507             :          END IF
    1508   120509174 :          IF (n_a + n_b > n_c + n_d) THEN
    1509    42117012 :             perm_case = perm_case + 4
    1510             :          END IF
    1511             :          SELECT CASE (perm_case)
    1512             :          CASE (1)
    1513   141431501 :             DO i = 1, nproducts
    1514   399466720 :                A = pgf_product_list(i)%ra
    1515   399466720 :                B = pgf_product_list(i)%rb
    1516   399466720 :                C = pgf_product_list(i)%rc
    1517   399466720 :                D = pgf_product_list(i)%rd
    1518    99866680 :                Rho = pgf_product_list(i)%Rho
    1519    99866680 :                RhoInv = pgf_product_list(i)%RhoInv
    1520   399466720 :                P = pgf_product_list(i)%P
    1521   399466720 :                Q = pgf_product_list(i)%Q
    1522   399466720 :                W = pgf_product_list(i)%W
    1523   399466720 :                AB = pgf_product_list(i)%AB
    1524   399466720 :                CD = pgf_product_list(i)%CD
    1525    99866680 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1526   400738693 :                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    99866680 :                                        P, Q, W, m_max, F)
    1530    99866680 :                CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1531  1984571302 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1532   141431501 :                neris = neris + mysize
    1533             :             END DO
    1534  1035486171 :             DO i = 1, mysize
    1535  1035486171 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1536             :             END DO
    1537    41564821 :             tmp_max = tmp_max*max_contraction
    1538    41564821 :             IF (tmp_max < eps_schwarz) THEN
    1539    35306896 :                RETURN
    1540             :             END IF
    1541             : 
    1542    98005966 :             DO i = 1, ncoa
    1543    64644506 :                p1 = (i - 1)*ncob
    1544   180752712 :                DO j = 1, ncob
    1545    82746746 :                   p2 = (p1 + j - 1)*ncoc
    1546   480999868 :                   DO k = 1, ncoc
    1547   333608616 :                      p3 = (p2 + k - 1)*ncod
    1548  1280164624 :                      DO l = 1, ncod
    1549   863809262 :                         p4 = p3 + l
    1550  1197417878 :                         work2(i, j, k, l) = work(p4)
    1551             :                      END DO
    1552             :                   END DO
    1553             :                END DO
    1554             :             END DO
    1555             :          CASE (2)
    1556    26534171 :             DO i = 1, nproducts
    1557    71795068 :                A = pgf_product_list(i)%ra
    1558    71795068 :                B = pgf_product_list(i)%rb
    1559    71795068 :                C = pgf_product_list(i)%rc
    1560    71795068 :                D = pgf_product_list(i)%rd
    1561    17948767 :                Rho = pgf_product_list(i)%Rho
    1562    17948767 :                RhoInv = pgf_product_list(i)%RhoInv
    1563    71795068 :                P = pgf_product_list(i)%P
    1564    71795068 :                Q = pgf_product_list(i)%Q
    1565    71795068 :                W = pgf_product_list(i)%W
    1566    71795068 :                AB = pgf_product_list(i)%AB
    1567    71795068 :                CD = pgf_product_list(i)%CD
    1568    17948767 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1569    86234732 :                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    17948767 :                                        P, Q, W, m_max, F)
    1574             : 
    1575    17948767 :                CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
    1576   605842628 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1577    26534171 :                neris = neris + mysize
    1578             :             END DO
    1579   384958999 :             DO i = 1, mysize
    1580   384958999 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1581             :             END DO
    1582     8585404 :             tmp_max = tmp_max*max_contraction
    1583     8585404 :             IF (tmp_max < eps_schwarz) THEN
    1584             :                RETURN
    1585             :             END IF
    1586             : 
    1587    28648635 :             DO j = 1, ncob
    1588    22384187 :                p1 = (j - 1)*ncoa
    1589    53926272 :                DO i = 1, ncoa
    1590    25277637 :                   p2 = (p1 + i - 1)*ncoc
    1591   157177386 :                   DO k = 1, ncoc
    1592   109515562 :                      p3 = (p2 + k - 1)*ncod
    1593   438493471 :                      DO l = 1, ncod
    1594   303700272 :                         p4 = p3 + l
    1595   413215834 :                         work2(i, j, k, l) = work(p4)
    1596             :                      END DO
    1597             :                   END DO
    1598             :                END DO
    1599             :             END DO
    1600             :          CASE (3)
    1601    66677540 :             DO i = 1, nproducts
    1602   178599768 :                A = pgf_product_list(i)%ra
    1603   178599768 :                B = pgf_product_list(i)%rb
    1604   178599768 :                C = pgf_product_list(i)%rc
    1605   178599768 :                D = pgf_product_list(i)%rd
    1606    44649942 :                Rho = pgf_product_list(i)%Rho
    1607    44649942 :                RhoInv = pgf_product_list(i)%RhoInv
    1608   178599768 :                P = pgf_product_list(i)%P
    1609   178599768 :                Q = pgf_product_list(i)%Q
    1610   178599768 :                W = pgf_product_list(i)%W
    1611   178599768 :                AB = pgf_product_list(i)%AB
    1612   178599768 :                CD = pgf_product_list(i)%CD
    1613    44649942 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1614   170048652 :                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    44649942 :                                        P, Q, W, m_max, F)
    1619             : 
    1620    44649942 :                CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
    1621   684621453 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1622    66677540 :                neris = neris + mysize
    1623             :             END DO
    1624   416183866 :             DO i = 1, mysize
    1625   416183866 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1626             :             END DO
    1627    22027598 :             tmp_max = tmp_max*max_contraction
    1628    22027598 :             IF (tmp_max < eps_schwarz) THEN
    1629             :                RETURN
    1630             :             END IF
    1631             : 
    1632    48593868 :             DO i = 1, ncoa
    1633    32548902 :                p1 = (i - 1)*ncob
    1634    88319514 :                DO j = 1, ncob
    1635    39725646 :                   p2 = (p1 + j - 1)*ncod
    1636   257676826 :                   DO l = 1, ncod
    1637   185402278 :                      p3 = (p2 + l - 1)*ncoc
    1638   544901808 :                      DO k = 1, ncoc
    1639   319773884 :                         p4 = p3 + k
    1640   505176162 :                         work2(i, j, k, l) = work(p4)
    1641             :                      END DO
    1642             :                   END DO
    1643             :                END DO
    1644             :             END DO
    1645             :          CASE (4)
    1646    17247580 :             DO i = 1, nproducts
    1647    44132964 :                A = pgf_product_list(i)%ra
    1648    44132964 :                B = pgf_product_list(i)%rb
    1649    44132964 :                C = pgf_product_list(i)%rc
    1650    44132964 :                D = pgf_product_list(i)%rd
    1651    11033241 :                Rho = pgf_product_list(i)%Rho
    1652    11033241 :                RhoInv = pgf_product_list(i)%RhoInv
    1653    44132964 :                P = pgf_product_list(i)%P
    1654    44132964 :                Q = pgf_product_list(i)%Q
    1655    44132964 :                W = pgf_product_list(i)%W
    1656    44132964 :                AB = pgf_product_list(i)%AB
    1657    44132964 :                CD = pgf_product_list(i)%CD
    1658    11033241 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1659    50683186 :                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    11033241 :                                        P, Q, W, m_max, F)
    1664             : 
    1665    11033241 :                CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1666   314809139 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1667    17247580 :                neris = neris + mysize
    1668             :             END DO
    1669   227952835 :             DO i = 1, mysize
    1670   227952835 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1671             :             END DO
    1672     6214339 :             tmp_max = tmp_max*max_contraction
    1673     6214339 :             IF (tmp_max < eps_schwarz) THEN
    1674             :                RETURN
    1675             :             END IF
    1676             : 
    1677    21602369 :             DO j = 1, ncob
    1678    16850830 :                p1 = (j - 1)*ncoa
    1679    41221321 :                DO i = 1, ncoa
    1680    19618952 :                   p2 = (p1 + i - 1)*ncod
    1681   130747996 :                   DO l = 1, ncod
    1682    94278214 :                      p3 = (p2 + l - 1)*ncoc
    1683   298004828 :                      DO k = 1, ncoc
    1684   184107662 :                         p4 = p3 + k
    1685   278385876 :                         work2(i, j, k, l) = work(p4)
    1686             :                      END DO
    1687             :                   END DO
    1688             :                END DO
    1689             :             END DO
    1690             :          CASE (5)
    1691    70613739 :             DO i = 1, nproducts
    1692   189757092 :                A = pgf_product_list(i)%ra
    1693   189757092 :                B = pgf_product_list(i)%rb
    1694   189757092 :                C = pgf_product_list(i)%rc
    1695   189757092 :                D = pgf_product_list(i)%rd
    1696    47439273 :                Rho = pgf_product_list(i)%Rho
    1697    47439273 :                RhoInv = pgf_product_list(i)%RhoInv
    1698   189757092 :                P = pgf_product_list(i)%P
    1699   189757092 :                Q = pgf_product_list(i)%Q
    1700   189757092 :                W = pgf_product_list(i)%W
    1701   189757092 :                AB = pgf_product_list(i)%AB
    1702   189757092 :                CD = pgf_product_list(i)%CD
    1703    47439273 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1704   183457452 :                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    47439273 :                                        Q, P, W, m_max, F)
    1709             : 
    1710    47439273 :                CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
    1711   853178142 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1712    70613739 :                neris = neris + mysize
    1713             :             END DO
    1714   527811096 :             DO i = 1, mysize
    1715   527811096 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1716             :             END DO
    1717    23174466 :             tmp_max = tmp_max*max_contraction
    1718    23174466 :             IF (tmp_max < eps_schwarz) THEN
    1719             :                RETURN
    1720             :             END IF
    1721             : 
    1722    42516134 :             DO k = 1, ncoc
    1723    25403701 :                p1 = (k - 1)*ncod
    1724    72036607 :                DO l = 1, ncod
    1725    29520473 :                   p2 = (p1 + l - 1)*ncoa
    1726   191074441 :                   DO i = 1, ncoa
    1727   136150267 :                      p3 = (p2 + i - 1)*ncob
    1728   574663559 :                      DO j = 1, ncob
    1729   408992819 :                         p4 = p3 + j
    1730   545143086 :                         work2(i, j, k, l) = work(p4)
    1731             :                      END DO
    1732             :                   END DO
    1733             :                END DO
    1734             :             END DO
    1735             :          CASE (6)
    1736    32223351 :             DO i = 1, nproducts
    1737    80616320 :                A = pgf_product_list(i)%ra
    1738    80616320 :                B = pgf_product_list(i)%rb
    1739    80616320 :                C = pgf_product_list(i)%rc
    1740    80616320 :                D = pgf_product_list(i)%rd
    1741    20154080 :                Rho = pgf_product_list(i)%Rho
    1742    20154080 :                RhoInv = pgf_product_list(i)%RhoInv
    1743    80616320 :                P = pgf_product_list(i)%P
    1744    80616320 :                Q = pgf_product_list(i)%Q
    1745    80616320 :                W = pgf_product_list(i)%W
    1746    80616320 :                AB = pgf_product_list(i)%AB
    1747    80616320 :                CD = pgf_product_list(i)%CD
    1748    20154080 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1749    73575793 :                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    20154080 :                                        Q, P, W, m_max, F)
    1754             : 
    1755    20154080 :                CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
    1756   299416795 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1757    32223351 :                neris = neris + mysize
    1758             :             END DO
    1759   205483427 :             DO i = 1, mysize
    1760   205483427 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1761             :             END DO
    1762    12069271 :             tmp_max = tmp_max*max_contraction
    1763    12069271 :             IF (tmp_max < eps_schwarz) THEN
    1764             :                RETURN
    1765             :             END IF
    1766             : 
    1767    19006913 :             DO k = 1, ncoc
    1768    11186272 :                p1 = (k - 1)*ncod
    1769    32213031 :                DO l = 1, ncod
    1770    13206118 :                   p2 = (p1 + l - 1)*ncob
    1771    91672370 :                   DO j = 1, ncob
    1772    67279980 :                      p3 = (p2 + j - 1)*ncoa
    1773   224242938 :                      DO i = 1, ncoa
    1774   143756840 :                         p4 = p3 + i
    1775   211036820 :                         work2(i, j, k, l) = work(p4)
    1776             :                      END DO
    1777             :                   END DO
    1778             :                END DO
    1779             :             END DO
    1780             :          CASE (7)
    1781    13823808 :             DO i = 1, nproducts
    1782    35872828 :                A = pgf_product_list(i)%ra
    1783    35872828 :                B = pgf_product_list(i)%rb
    1784    35872828 :                C = pgf_product_list(i)%rc
    1785    35872828 :                D = pgf_product_list(i)%rd
    1786     8968207 :                Rho = pgf_product_list(i)%Rho
    1787     8968207 :                RhoInv = pgf_product_list(i)%RhoInv
    1788    35872828 :                P = pgf_product_list(i)%P
    1789    35872828 :                Q = pgf_product_list(i)%Q
    1790    35872828 :                W = pgf_product_list(i)%W
    1791    35872828 :                AB = pgf_product_list(i)%AB
    1792    35872828 :                CD = pgf_product_list(i)%CD
    1793     8968207 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1794    48555662 :                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     8968207 :                                        Q, P, W, m_max, F)
    1799             : 
    1800     8968207 :                CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
    1801   469723695 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1802    13823808 :                neris = neris + mysize
    1803             :             END DO
    1804   332686953 :             DO i = 1, mysize
    1805   332686953 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1806             :             END DO
    1807     4855601 :             tmp_max = tmp_max*max_contraction
    1808     4855601 :             IF (tmp_max < eps_schwarz) THEN
    1809             :                RETURN
    1810             :             END IF
    1811             : 
    1812    15284209 :             DO l = 1, ncod
    1813    11876293 :                p1 = (l - 1)*ncoc
    1814    28821546 :                DO k = 1, ncoc
    1815    13537337 :                   p2 = (p1 + k - 1)*ncoa
    1816    96511146 :                   DO i = 1, ncoa
    1817    71097516 :                      p3 = (p2 + i - 1)*ncob
    1818   349990923 :                      DO j = 1, ncob
    1819   265356070 :                         p4 = p3 + j
    1820   336453586 :                         work2(i, j, k, l) = work(p4)
    1821             :                      END DO
    1822             :                   END DO
    1823             :                END DO
    1824             :             END DO
    1825             :          CASE (8)
    1826     4674682 :             DO i = 1, nproducts
    1827    10628032 :                A = pgf_product_list(i)%ra
    1828    10628032 :                B = pgf_product_list(i)%rb
    1829    10628032 :                C = pgf_product_list(i)%rc
    1830    10628032 :                D = pgf_product_list(i)%rd
    1831     2657008 :                Rho = pgf_product_list(i)%Rho
    1832     2657008 :                RhoInv = pgf_product_list(i)%RhoInv
    1833    10628032 :                P = pgf_product_list(i)%P
    1834    10628032 :                Q = pgf_product_list(i)%Q
    1835    10628032 :                W = pgf_product_list(i)%W
    1836    10628032 :                AB = pgf_product_list(i)%AB
    1837    10628032 :                CD = pgf_product_list(i)%CD
    1838     2657008 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1839    15139352 :                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     2657008 :                                        Q, P, W, m_max, F)
    1844             : 
    1845     2657008 :                CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
    1846   131880226 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1847     4674682 :                neris = neris + mysize
    1848             :             END DO
    1849   106386990 :             DO i = 1, mysize
    1850   106386990 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1851             :             END DO
    1852     2017674 :             tmp_max = tmp_max*max_contraction
    1853     2017674 :             IF (tmp_max < eps_schwarz) THEN
    1854             :                RETURN
    1855             :             END IF
    1856             : 
    1857   127605278 :             DO l = 1, ncod
    1858     5545348 :                p1 = (l - 1)*ncoc
    1859    12701772 :                DO k = 1, ncoc
    1860     5605668 :                   p2 = (p1 + k - 1)*ncob
    1861    45854908 :                   DO j = 1, ncob
    1862    34703892 :                      p3 = (p2 + j - 1)*ncoa
    1863   122356674 :                      DO i = 1, ncoa
    1864    82047114 :                         p4 = p3 + i
    1865   116751006 :                         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   101540891 :          DO i = 1, nproducts
    1873   279130876 :             A = pgf_product_list(i)%ra
    1874   279130876 :             B = pgf_product_list(i)%rb
    1875   279130876 :             C = pgf_product_list(i)%rc
    1876   279130876 :             D = pgf_product_list(i)%rd
    1877    69782719 :             Rho = pgf_product_list(i)%Rho
    1878    69782719 :             RhoInv = pgf_product_list(i)%RhoInv
    1879   279130876 :             P = pgf_product_list(i)%P
    1880   279130876 :             Q = pgf_product_list(i)%Q
    1881   279130876 :             W = pgf_product_list(i)%W
    1882    69782719 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1883   139565438 :             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    69782719 :                                     P, Q, W, m_max, F)
    1888    69782719 :             work(1) = work(1) + F(1)
    1889   101540891 :             neris = neris + mysize
    1890             :          END DO
    1891    31758172 :          work2(1, 1, 1, 1) = work(1)
    1892    31758172 :          tmp_max = max_contraction*ABS(work(1))
    1893    31758172 :          IF (tmp_max < eps_schwarz) RETURN
    1894             :       END IF
    1895             : 
    1896   116960450 :       IF (tmp_max < eps_schwarz) RETURN
    1897  7738251460 :       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   116960450 :                     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  7738251460 :                     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   322499917 :    SUBROUTINE build_quartet_data(libint, A, B, C, D, &
    1940             :                                  ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1941   322499917 :                                  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   322499917 :       CALL cp_libint_set_params_eri(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, P, Q, W, m_max, F)
    1951   322499917 :    END SUBROUTINE build_quartet_data
    1952             : 
    1953             : END MODULE hfx_libint_interface

Generated by: LCOV version 1.15