LCOV - code coverage report
Current view: top level - src - hfx_libint_interface.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:3158929) Lines: 951 983 96.7 %
Date: 2025-07-22 22:41:15 Functions: 6 6 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief 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    82517808 :    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    82517808 :       Zeta = Zeta_A + Zeta_B
      98    82517808 :       ZetaInv = 1.0_dp/Zeta
      99    82517808 :       Eta = Zeta_C + Zeta_D
     100    82517808 :       EtaInv = 1.0_dp/Eta
     101    82517808 :       ZetapEtaInv = Zeta + Eta
     102    82517808 :       ZetapEtaInv = 1.0_dp/ZetapEtaInv
     103    82517808 :       Rho = Zeta*Eta*ZetapEtaInv
     104    82517808 :       RhoInv = 1.0_dp/Rho
     105             : 
     106   330071232 :       DO i = 1, 3
     107   247553424 :          P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv
     108   247553424 :          Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv
     109   247553424 :          AB(i) = A(i) - B(i)
     110   247553424 :          CD(i) = C(i) - D(i)
     111   247553424 :          PQ(i) = P(i) - Q(i)
     112   330071232 :          W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv
     113             :       END DO
     114             : 
     115   330071232 :       AB2 = DOT_PRODUCT(AB, AB)
     116   330071232 :       CD2 = DOT_PRODUCT(CD, CD)
     117   330071232 :       PQ2 = DOT_PRODUCT(PQ, PQ)
     118             : 
     119    82517808 :       S1234 = EXP((-Zeta_A*Zeta_B*ZetaInv*AB2) + (-Zeta_C*Zeta_D*EtaInv*CD2))
     120    82517808 :       T = Rho*PQ2
     121             : 
     122    95404398 :       SELECT CASE (potential_parameter%potential_type)
     123             :       CASE (do_potential_truncated)
     124    12886590 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     125    12886590 :          R1 = R11
     126    12886590 :          R2 = R22
     127    12886590 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     128           0 :             RETURN
     129             :          END IF
     130    12886590 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     131    12886590 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     132    12886590 :          factor = 2.0_dp*Pi*RhoInv
     133             :       CASE (do_potential_coulomb)
     134    62557986 :          CALL fgamma(m_max, T, F(1))
     135    62557986 :          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    96586704 :          IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
     256             :       END SELECT
     257             : 
     258    82517808 :       tmp = (Pi*ZetapEtaInv)**3
     259    82517808 :       factor = factor*S1234*SQRT(tmp)
     260             : 
     261   318775392 :       DO i = 1, m_max + 1
     262   318775392 :          F(i) = F(i)*factor
     263             :       END DO
     264             : 
     265    82517808 :       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    98775197 :    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    98775197 :                                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    98775197 :                                           ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F)
     313             : 
     314    98775197 :    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    31938139 :    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    31938139 :                                  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    31938139 :                                  sphi_a, sphi_b, sphi_c, sphi_d, &
     388    31938139 :                                  work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
     389    31938139 :                                  use_virial, work_virial, work2_virial, primitives_tmp_virial, &
     390    31938139 :                                  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    31938139 :       m_max = n_a + n_b + n_c + n_d
     433    31938139 :       m_max = m_max + 1
     434    31938139 :       mysize = ncoa*ncob*ncoc*ncod
     435    63876278 :       a_mysize = mysize
     436             : 
     437  4183493011 :       work = 0.0_dp
     438  8057992015 :       work2 = 0.0_dp
     439             : 
     440    31938139 :       IF (use_virial) THEN
     441   850401248 :          work_virial = 0.0_dp
     442  1465200824 :          work2_virial = 0.0_dp
     443             :       END IF
     444             : 
     445    31938139 :       perm_case = 1
     446    31938139 :       IF (n_a < n_b) THEN
     447     2941150 :          perm_case = perm_case + 1
     448             :       END IF
     449    31938139 :       IF (n_c < n_d) THEN
     450     5065824 :          perm_case = perm_case + 2
     451             :       END IF
     452    31938139 :       IF (n_a + n_b > n_c + n_d) THEN
     453     5760416 :          perm_case = perm_case + 4
     454             :       END IF
     455             : 
     456             :       SELECT CASE (perm_case)
     457             :       CASE (1)
     458    80843629 :          DO i = 1, nproducts
     459   241048328 :             A = pgf_product_list(i)%ra
     460   241048328 :             B = pgf_product_list(i)%rb
     461   241048328 :             C = pgf_product_list(i)%rc
     462   241048328 :             D = pgf_product_list(i)%rd
     463    60262082 :             Rho = pgf_product_list(i)%Rho
     464    60262082 :             RhoInv = pgf_product_list(i)%RhoInv
     465   241048328 :             P = pgf_product_list(i)%P
     466   241048328 :             Q = pgf_product_list(i)%Q
     467   241048328 :             W = pgf_product_list(i)%W
     468   241048328 :             AB = pgf_product_list(i)%AB
     469   241048328 :             CD = pgf_product_list(i)%CD
     470    60262082 :             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    60262082 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     476    60262082 :             CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize)
     477   241048328 :             DO k = 4, 6
     478  1856400941 :                DO j = 1, mysize
     479             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     480             :                                                work_forces(j, k + 3) + &
     481  1796138859 :                                                work_forces(j, k + 6))
     482             :                END DO
     483             :             END DO
     484   783407066 :             DO k = 1, 12
     485  7244817518 :                DO j = 1, mysize
     486  7184555436 :                   work(j, k) = work(j, k) + work_forces(j, k)
     487             :                END DO
     488             :             END DO
     489    60262082 :             neris = neris + 12*mysize
     490    80843629 :             IF (use_virial) THEN
     491    10723386 :                CALL real_to_scaled(scoord(1:3), A, cell)
     492    10723386 :                CALL real_to_scaled(scoord(4:6), B, cell)
     493    10723386 :                CALL real_to_scaled(scoord(7:9), C, cell)
     494    10723386 :                CALL real_to_scaled(scoord(10:12), D, cell)
     495   139404018 :                DO k = 1, 12
     496  1164286866 :                   DO j = 1, mysize
     497  4228212024 :                      DO m = 1, 3
     498  4099531392 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     499             :                      END DO
     500             :                   END DO
     501             :                END DO
     502             :             END IF
     503             :          END DO
     504             : 
     505   267560111 :          DO n = 1, 12
     506             :             tmp_max = 0.0_dp
     507  2093421744 :             DO i = 1, mysize
     508  2093421744 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     509             :             END DO
     510   246978564 :             tmp_max = tmp_max*max_contraction
     511   246978564 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     512             : 
     513   608116199 :             DO i = 1, ncoa
     514   340556088 :                p1 = (i - 1)*ncob
     515   962323716 :                DO j = 1, ncob
     516   374789064 :                   p2 = (p1 + j - 1)*ncoc
     517  1720342260 :                   DO k = 1, ncoc
     518  1004997108 :                      p3 = (p2 + k - 1)*ncod
     519  3226229352 :                      DO l = 1, ncod
     520  1846443180 :                         p4 = p3 + l
     521  2851440288 :                         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    20581547 :          IF (use_virial) THEN
     528     6233656 :             DO n = 1, 12
     529             :                tmp_max_virial = 0.0_dp
     530   112731564 :                DO i = 1, mysize
     531             :                   tmp_max_virial = MAX(tmp_max_virial, &
     532   112731564 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     533             :                END DO
     534     5754144 :                tmp_max_virial = tmp_max_virial*max_contraction
     535     5754144 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     536             : 
     537    16527736 :                DO i = 1, ncoa
     538    10294080 :                   p1 = (i - 1)*ncob
     539    29965200 :                   DO j = 1, ncob
     540    13916976 :                      p2 = (p1 + j - 1)*ncoc
     541    69604716 :                      DO k = 1, ncoc
     542    45393660 :                         p3 = (p2 + k - 1)*ncod
     543   166288056 :                         DO l = 1, ncod
     544   106977420 :                            p4 = p3 + l
     545   473303340 :                            work2_virial(i, j, k, l, full_perm1(n), 1:3) = work_virial(p4, n, 1:3)
     546             :                         END DO
     547             :                      END DO
     548             :                   END DO
     549             :                END DO
     550             :             END DO
     551             :          END IF
     552             :       CASE (2)
     553     5720096 :          DO i = 1, nproducts
     554    18191652 :             A = pgf_product_list(i)%ra
     555    18191652 :             B = pgf_product_list(i)%rb
     556    18191652 :             C = pgf_product_list(i)%rc
     557    18191652 :             D = pgf_product_list(i)%rd
     558     4547913 :             Rho = pgf_product_list(i)%Rho
     559     4547913 :             RhoInv = pgf_product_list(i)%RhoInv
     560    18191652 :             P = pgf_product_list(i)%P
     561    18191652 :             Q = pgf_product_list(i)%Q
     562    18191652 :             W = pgf_product_list(i)%W
     563    18191652 :             AB = pgf_product_list(i)%AB
     564    18191652 :             CD = pgf_product_list(i)%CD
     565     4547913 :             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     4547913 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     571     4547913 :             CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize)
     572    18191652 :             DO k = 4, 6
     573   322971027 :                DO j = 1, mysize
     574             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     575             :                                                work_forces(j, k + 3) + &
     576   318423114 :                                                work_forces(j, k + 6))
     577             :                END DO
     578             :             END DO
     579    59122869 :             DO k = 1, 12
     580  1278240369 :                DO j = 1, mysize
     581  1273692456 :                   work(j, k) = work(j, k) + work_forces(j, k)
     582             :                END DO
     583             :             END DO
     584     4547913 :             neris = neris + 12*mysize
     585     5720096 :             IF (use_virial) THEN
     586      674510 :                CALL real_to_scaled(scoord(1:3), B, cell)
     587      674510 :                CALL real_to_scaled(scoord(4:6), A, cell)
     588      674510 :                CALL real_to_scaled(scoord(7:9), C, cell)
     589      674510 :                CALL real_to_scaled(scoord(10:12), D, cell)
     590     8768630 :                DO k = 1, 12
     591   180406874 :                   DO j = 1, mysize
     592   694647096 :                      DO m = 1, 3
     593   686552976 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     594             :                      END DO
     595             :                   END DO
     596             :                END DO
     597             :             END IF
     598             : 
     599             :          END DO
     600    15238379 :          DO n = 1, 12
     601             :             tmp_max = 0.0_dp
     602   364671456 :             DO i = 1, mysize
     603   364671456 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     604             :             END DO
     605    14066196 :             tmp_max = tmp_max*max_contraction
     606    14066196 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     607             : 
     608    59852675 :             DO j = 1, ncob
     609    44614296 :                p1 = (j - 1)*ncoa
     610   104731044 :                DO i = 1, ncoa
     611    46050552 :                   p2 = (p1 + i - 1)*ncoc
     612   256141260 :                   DO k = 1, ncoc
     613   165476412 :                      p3 = (p2 + k - 1)*ncod
     614   562132224 :                      DO l = 1, ncod
     615   350605260 :                         p4 = p3 + l
     616   516081672 :                         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     1172183 :          IF (use_virial) THEN
     623     1273922 :             DO n = 1, 12
     624             :                tmp_max_virial = 0.0_dp
     625    34681632 :                DO i = 1, mysize
     626             :                   tmp_max_virial = MAX(tmp_max_virial, &
     627    34681632 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     628             :                END DO
     629     1175928 :                tmp_max_virial = tmp_max_virial*max_contraction
     630     1175928 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     631             : 
     632     5078186 :                DO j = 1, ncob
     633     3804264 :                   p1 = (j - 1)*ncoa
     634     8968776 :                   DO i = 1, ncoa
     635     3988584 :                      p2 = (p1 + i - 1)*ncoc
     636    22271544 :                      DO k = 1, ncoc
     637    14478696 :                         p3 = (p2 + k - 1)*ncod
     638    51972984 :                         DO l = 1, ncod
     639    33505704 :                            p4 = p3 + l
     640   148501512 :                            work2_virial(i, j, k, l, full_perm2(n), 1:3) = work_virial(p4, n, 1:3)
     641             :                         END DO
     642             :                      END DO
     643             :                   END DO
     644             :                END DO
     645             :             END DO
     646             :          END IF
     647             :       CASE (3)
     648    16600248 :          DO i = 1, nproducts
     649    51281176 :             A = pgf_product_list(i)%ra
     650    51281176 :             B = pgf_product_list(i)%rb
     651    51281176 :             C = pgf_product_list(i)%rc
     652    51281176 :             D = pgf_product_list(i)%rd
     653    12820294 :             Rho = pgf_product_list(i)%Rho
     654    12820294 :             RhoInv = pgf_product_list(i)%RhoInv
     655    51281176 :             P = pgf_product_list(i)%P
     656    51281176 :             Q = pgf_product_list(i)%Q
     657    51281176 :             W = pgf_product_list(i)%W
     658    51281176 :             AB = pgf_product_list(i)%AB
     659    51281176 :             CD = pgf_product_list(i)%CD
     660    12820294 :             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    12820294 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     666    12820294 :             CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize)
     667    51281176 :             DO k = 4, 6
     668   453180784 :                DO j = 1, mysize
     669             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     670             :                                                work_forces(j, k + 3) + &
     671   440360490 :                                                work_forces(j, k + 6))
     672             :                END DO
     673             :             END DO
     674   166663822 :             DO k = 1, 12
     675  1774262254 :                DO j = 1, mysize
     676  1761441960 :                   work(j, k) = work(j, k) + work_forces(j, k)
     677             :                END DO
     678             :             END DO
     679    12820294 :             neris = neris + 12*mysize
     680    16600248 :             IF (use_virial) THEN
     681     1825532 :                CALL real_to_scaled(scoord(1:3), A, cell)
     682     1825532 :                CALL real_to_scaled(scoord(4:6), B, cell)
     683     1825532 :                CALL real_to_scaled(scoord(7:9), D, cell)
     684     1825532 :                CALL real_to_scaled(scoord(10:12), C, cell)
     685    23731916 :                DO k = 1, 12
     686   205909556 :                   DO j = 1, mysize
     687   750616944 :                      DO m = 1, 3
     688   728710560 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     689             :                      END DO
     690             :                   END DO
     691             :                END DO
     692             :             END IF
     693             : 
     694             :          END DO
     695    49139402 :          DO n = 1, 12
     696             :             tmp_max = 0.0_dp
     697   523588596 :             DO i = 1, mysize
     698   523588596 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     699             :             END DO
     700    45359448 :             tmp_max = tmp_max*max_contraction
     701    45359448 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     702             : 
     703   127791722 :             DO i = 1, ncoa
     704    78652320 :                p1 = (i - 1)*ncob
     705   209918304 :                DO j = 1, ncob
     706    85906536 :                   p2 = (p1 + j - 1)*ncod
     707   495686244 :                   DO l = 1, ncod
     708   331127388 :                      p3 = (p2 + l - 1)*ncoc
     709   895263072 :                      DO k = 1, ncoc
     710   478229148 :                         p4 = p3 + k
     711   809356536 :                         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     3779954 :          IF (use_virial) THEN
     718     1829516 :             DO n = 1, 12
     719             :                tmp_max_virial = 0.0_dp
     720    32066520 :                DO i = 1, mysize
     721             :                   tmp_max_virial = MAX(tmp_max_virial, &
     722    32066520 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     723             :                END DO
     724     1688784 :                tmp_max_virial = tmp_max_virial*max_contraction
     725     1688784 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     726             : 
     727     5285348 :                DO i = 1, ncoa
     728     3455832 :                   p1 = (i - 1)*ncob
     729     9372288 :                   DO j = 1, ncob
     730     4227672 :                      p2 = (p1 + j - 1)*ncod
     731    25806840 :                      DO l = 1, ncod
     732    18123336 :                         p3 = (p2 + l - 1)*ncoc
     733    52728744 :                         DO k = 1, ncoc
     734    30377736 :                            p4 = p3 + k
     735   139634280 :                            work2_virial(i, j, k, l, full_perm3(n), 1:3) = work_virial(p4, n, 1:3)
     736             :                         END DO
     737             :                      END DO
     738             :                   END DO
     739             :                END DO
     740             :             END DO
     741             :          END IF
     742             :       CASE (4)
     743     2884127 :          DO i = 1, nproducts
     744     8960352 :             A = pgf_product_list(i)%ra
     745     8960352 :             B = pgf_product_list(i)%rb
     746     8960352 :             C = pgf_product_list(i)%rc
     747     8960352 :             D = pgf_product_list(i)%rd
     748     2240088 :             Rho = pgf_product_list(i)%Rho
     749     2240088 :             RhoInv = pgf_product_list(i)%RhoInv
     750     8960352 :             P = pgf_product_list(i)%P
     751     8960352 :             Q = pgf_product_list(i)%Q
     752     8960352 :             W = pgf_product_list(i)%W
     753     8960352 :             AB = pgf_product_list(i)%AB
     754     8960352 :             CD = pgf_product_list(i)%CD
     755     2240088 :             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     2240088 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     760     2240088 :             CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize)
     761     8960352 :             DO k = 4, 6
     762   121804773 :                DO j = 1, mysize
     763             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     764             :                                                work_forces(j, k + 3) + &
     765   119564685 :                                                work_forces(j, k + 6))
     766             :                END DO
     767             :             END DO
     768    29121144 :             DO k = 1, 12
     769   480498828 :                DO j = 1, mysize
     770   478258740 :                   work(j, k) = work(j, k) + work_forces(j, k)
     771             :                END DO
     772             :             END DO
     773     2240088 :             neris = neris + 12*mysize
     774     2884127 :             IF (use_virial) THEN
     775      327543 :                CALL real_to_scaled(scoord(1:3), B, cell)
     776      327543 :                CALL real_to_scaled(scoord(4:6), A, cell)
     777      327543 :                CALL real_to_scaled(scoord(7:9), D, cell)
     778      327543 :                CALL real_to_scaled(scoord(10:12), C, cell)
     779     4258059 :                DO k = 1, 12
     780    59544879 :                   DO j = 1, mysize
     781   225077796 :                      DO m = 1, 3
     782   221147280 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     783             :                      END DO
     784             :                   END DO
     785             :                END DO
     786             :             END IF
     787             : 
     788             :          END DO
     789     8372507 :          DO n = 1, 12
     790             :             tmp_max = 0.0_dp
     791   154357368 :             DO i = 1, mysize
     792   154357368 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     793             :             END DO
     794     7728468 :             tmp_max = tmp_max*max_contraction
     795     7728468 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     796             : 
     797    32553671 :             DO j = 1, ncob
     798    24181164 :                p1 = (j - 1)*ncoa
     799    57240204 :                DO i = 1, ncoa
     800    25330572 :                   p2 = (p1 + i - 1)*ncod
     801   146886156 :                   DO l = 1, ncod
     802    97374420 :                      p3 = (p2 + l - 1)*ncoc
     803   269333892 :                      DO k = 1, ncoc
     804   146628900 :                         p4 = p3 + k
     805   244003320 :                         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      644039 :          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    17319892 :          DO i = 1, nproducts
     838    52954920 :             A = pgf_product_list(i)%ra
     839    52954920 :             B = pgf_product_list(i)%rb
     840    52954920 :             C = pgf_product_list(i)%rc
     841    52954920 :             D = pgf_product_list(i)%rd
     842    13238730 :             Rho = pgf_product_list(i)%Rho
     843    13238730 :             RhoInv = pgf_product_list(i)%RhoInv
     844    52954920 :             P = pgf_product_list(i)%P
     845    52954920 :             Q = pgf_product_list(i)%Q
     846    52954920 :             W = pgf_product_list(i)%W
     847    52954920 :             AB = pgf_product_list(i)%AB
     848    52954920 :             CD = pgf_product_list(i)%CD
     849    13238730 :             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    13238730 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     854    13238730 :             CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize)
     855    52954920 :             DO k = 4, 6
     856   505932669 :                DO j = 1, mysize
     857             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     858             :                                                work_forces(j, k + 3) + &
     859   492693939 :                                                work_forces(j, k + 6))
     860             :                END DO
     861             :             END DO
     862   172103490 :             DO k = 1, 12
     863  1984014486 :                DO j = 1, mysize
     864  1970775756 :                   work(j, k) = work(j, k) + work_forces(j, k)
     865             :                END DO
     866             :             END DO
     867    13238730 :             neris = neris + 12*mysize
     868    17319892 :             IF (use_virial) THEN
     869     1995603 :                CALL real_to_scaled(scoord(1:3), C, cell)
     870     1995603 :                CALL real_to_scaled(scoord(4:6), D, cell)
     871     1995603 :                CALL real_to_scaled(scoord(7:9), A, cell)
     872     1995603 :                CALL real_to_scaled(scoord(10:12), B, cell)
     873    25942839 :                DO k = 1, 12
     874   267747675 :                   DO j = 1, mysize
     875   991166580 :                      DO m = 1, 3
     876   967219344 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     877             :                      END DO
     878             :                   END DO
     879             :                END DO
     880             :             END IF
     881             : 
     882             :          END DO
     883    53055106 :          DO n = 1, 12
     884             :             tmp_max = 0.0_dp
     885   581580192 :             DO i = 1, mysize
     886   581580192 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     887             :             END DO
     888    48973944 :             tmp_max = tmp_max*max_contraction
     889    48973944 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     890             : 
     891   118885558 :             DO k = 1, ncoc
     892    65830452 :                p1 = (k - 1)*ncod
     893   184281864 :                DO l = 1, ncod
     894    69477468 :                   p2 = (p1 + l - 1)*ncoa
     895   395966568 :                   DO i = 1, ncoa
     896   260658648 :                      p3 = (p2 + i - 1)*ncob
     897   862742364 :                      DO j = 1, ncob
     898   532606248 :                         p4 = p3 + j
     899   793264896 :                         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     4081162 :          IF (use_virial) THEN
     906     2252003 :             DO n = 1, 12
     907             :                tmp_max_virial = 0.0_dp
     908    44899908 :                DO i = 1, mysize
     909             :                   tmp_max_virial = MAX(tmp_max_virial, &
     910    44899908 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     911             :                END DO
     912     2078772 :                tmp_max_virial = tmp_max_virial*max_contraction
     913     2078772 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     914             : 
     915     5747327 :                DO k = 1, ncoc
     916     3495324 :                   p1 = (k - 1)*ncod
     917     9490044 :                   DO l = 1, ncod
     918     3915948 :                      p2 = (p1 + l - 1)*ncoa
     919    22978104 :                      DO i = 1, ncoa
     920    15566832 :                         p3 = (p2 + i - 1)*ncob
     921    62303916 :                         DO j = 1, ncob
     922    42821136 :                            p4 = p3 + j
     923   186851376 :                            work2_virial(i, j, k, l, full_perm5(n), 1:3) = work_virial(p4, n, 1:3)
     924             :                         END DO
     925             :                      END DO
     926             :                   END DO
     927             :                END DO
     928             :             END DO
     929             :          END IF
     930             :       CASE (6)
     931     4432742 :          DO i = 1, nproducts
     932    13581276 :             A = pgf_product_list(i)%ra
     933    13581276 :             B = pgf_product_list(i)%rb
     934    13581276 :             C = pgf_product_list(i)%rc
     935    13581276 :             D = pgf_product_list(i)%rd
     936     3395319 :             Rho = pgf_product_list(i)%Rho
     937     3395319 :             RhoInv = pgf_product_list(i)%RhoInv
     938    13581276 :             P = pgf_product_list(i)%P
     939    13581276 :             Q = pgf_product_list(i)%Q
     940    13581276 :             W = pgf_product_list(i)%W
     941    13581276 :             AB = pgf_product_list(i)%AB
     942    13581276 :             CD = pgf_product_list(i)%CD
     943     3395319 :             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     3395319 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     949     3395319 :             CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize)
     950    13581276 :             DO k = 4, 6
     951   128089626 :                DO j = 1, mysize
     952             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     953             :                                                work_forces(j, k + 3) + &
     954   124694307 :                                                work_forces(j, k + 6))
     955             :                END DO
     956             :             END DO
     957    44139147 :             DO k = 1, 12
     958   502172547 :                DO j = 1, mysize
     959   498777228 :                   work(j, k) = work(j, k) + work_forces(j, k)
     960             :                END DO
     961             :             END DO
     962     3395319 :             neris = neris + 12*mysize
     963     4432742 :             IF (use_virial) THEN
     964      381872 :                CALL real_to_scaled(scoord(1:3), C, cell)
     965      381872 :                CALL real_to_scaled(scoord(4:6), D, cell)
     966      381872 :                CALL real_to_scaled(scoord(7:9), B, cell)
     967      381872 :                CALL real_to_scaled(scoord(10:12), A, cell)
     968     4964336 :                DO k = 1, 12
     969    51070436 :                   DO j = 1, mysize
     970   189006864 :                      DO m = 1, 3
     971   184424400 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
     972             :                      END DO
     973             :                   END DO
     974             :                END DO
     975             :             END IF
     976             : 
     977             :          END DO
     978    13486499 :          DO n = 1, 12
     979             :             tmp_max = 0.0_dp
     980   161255328 :             DO i = 1, mysize
     981   161255328 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     982             :             END DO
     983    12449076 :             tmp_max = tmp_max*max_contraction
     984    12449076 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     985             : 
     986    29118599 :             DO k = 1, ncoc
     987    15632100 :                p1 = (k - 1)*ncod
     988    46009644 :                DO l = 1, ncod
     989    17928468 :                   p2 = (p1 + l - 1)*ncob
     990   111931812 :                   DO j = 1, ncob
     991    78371244 :                      p3 = (p2 + j - 1)*ncoa
     992   245105964 :                      DO i = 1, ncoa
     993   148806252 :                         p4 = p3 + i
     994   227177496 :                         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     1037423 :          IF (use_virial) THEN
    1001      836784 :             DO n = 1, 12
    1002             :                tmp_max_virial = 0.0_dp
    1003    16323840 :                DO i = 1, mysize
    1004             :                   tmp_max_virial = MAX(tmp_max_virial, &
    1005    16323840 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1006             :                END DO
    1007      772416 :                tmp_max_virial = tmp_max_virial*max_contraction
    1008      772416 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1009             : 
    1010     1965552 :                DO k = 1, ncoc
    1011     1128768 :                   p1 = (k - 1)*ncod
    1012     3324864 :                   DO l = 1, ncod
    1013     1423680 :                      p2 = (p1 + l - 1)*ncob
    1014     9551424 :                      DO j = 1, ncob
    1015     6998976 :                         p3 = (p2 + j - 1)*ncoa
    1016    23974080 :                         DO i = 1, ncoa
    1017    15551424 :                            p4 = p3 + i
    1018    69204672 :                            work2_virial(i, j, k, l, full_perm6(n), 1:3) = work_virial(p4, n, 1:3)
    1019             :                         END DO
    1020             :                      END DO
    1021             :                   END DO
    1022             :                END DO
    1023             :             END DO
    1024             :          END IF
    1025             :       CASE (7)
    1026     2581909 :          DO i = 1, nproducts
    1027     8110332 :             A = pgf_product_list(i)%ra
    1028     8110332 :             B = pgf_product_list(i)%rb
    1029     8110332 :             C = pgf_product_list(i)%rc
    1030     8110332 :             D = pgf_product_list(i)%rd
    1031     2027583 :             Rho = pgf_product_list(i)%Rho
    1032     2027583 :             RhoInv = pgf_product_list(i)%RhoInv
    1033     8110332 :             P = pgf_product_list(i)%P
    1034     8110332 :             Q = pgf_product_list(i)%Q
    1035     8110332 :             W = pgf_product_list(i)%W
    1036     8110332 :             AB = pgf_product_list(i)%AB
    1037     8110332 :             CD = pgf_product_list(i)%CD
    1038     2027583 :             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     2027583 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1044     2027583 :             CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize)
    1045     8110332 :             DO k = 4, 6
    1046   195254919 :                DO j = 1, mysize
    1047             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1048             :                                                work_forces(j, k + 3) + &
    1049   193227336 :                                                work_forces(j, k + 6))
    1050             :                END DO
    1051             :             END DO
    1052    26358579 :             DO k = 1, 12
    1053   774936927 :                DO j = 1, mysize
    1054   772909344 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1055             :                END DO
    1056             :             END DO
    1057     2027583 :             neris = neris + 12*mysize
    1058     2581909 :             IF (use_virial) THEN
    1059      325942 :                CALL real_to_scaled(scoord(1:3), D, cell)
    1060      325942 :                CALL real_to_scaled(scoord(4:6), C, cell)
    1061      325942 :                CALL real_to_scaled(scoord(7:9), A, cell)
    1062      325942 :                CALL real_to_scaled(scoord(10:12), B, cell)
    1063     4237246 :                DO k = 1, 12
    1064   119756746 :                   DO j = 1, mysize
    1065   465989304 :                      DO m = 1, 3
    1066   462078000 :                         work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m)
    1067             :                      END DO
    1068             :                   END DO
    1069             :                END DO
    1070             :             END IF
    1071             : 
    1072             :          END DO
    1073     7206238 :          DO n = 1, 12
    1074             :             tmp_max = 0.0_dp
    1075   227976744 :             DO i = 1, mysize
    1076   227976744 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1077             :             END DO
    1078     6651912 :             tmp_max = tmp_max*max_contraction
    1079     6651912 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1080             : 
    1081    27686314 :             DO l = 1, ncod
    1082    20480076 :                p1 = (l - 1)*ncoc
    1083    47898912 :                DO k = 1, ncoc
    1084    20766924 :                   p2 = (p1 + k - 1)*ncoa
    1085   124440264 :                   DO i = 1, ncoa
    1086    83193264 :                      p3 = (p2 + i - 1)*ncob
    1087   325285020 :                      DO j = 1, ncob
    1088   221324832 :                         p4 = p3 + j
    1089   304518096 :                         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      554326 :          IF (use_virial) THEN
    1096      640536 :             DO n = 1, 12
    1097             :                tmp_max_virial = 0.0_dp
    1098    21929472 :                DO i = 1, mysize
    1099             :                   tmp_max_virial = MAX(tmp_max_virial, &
    1100    21929472 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1101             :                END DO
    1102      591264 :                tmp_max_virial = tmp_max_virial*max_contraction
    1103      591264 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1104             : 
    1105     2470920 :                DO l = 1, ncod
    1106     1830384 :                   p1 = (l - 1)*ncoc
    1107     4288896 :                   DO k = 1, ncoc
    1108     1867248 :                      p2 = (p1 + k - 1)*ncoa
    1109    10867968 :                      DO i = 1, ncoa
    1110     7170336 :                         p3 = (p2 + i - 1)*ncob
    1111    30375792 :                         DO j = 1, ncob
    1112    21338208 :                            p4 = p3 + j
    1113    92523168 :                            work2_virial(i, j, k, l, full_perm7(n), 1:3) = work_virial(p4, n, 1:3)
    1114             :                         END DO
    1115             :                      END DO
    1116             :                   END DO
    1117             :                END DO
    1118             :             END DO
    1119             :          END IF
    1120             :       CASE (8)
    1121      330693 :          DO i = 1, nproducts
    1122      972752 :             A = pgf_product_list(i)%ra
    1123      972752 :             B = pgf_product_list(i)%rb
    1124      972752 :             C = pgf_product_list(i)%rc
    1125      972752 :             D = pgf_product_list(i)%rd
    1126      243188 :             Rho = pgf_product_list(i)%Rho
    1127      243188 :             RhoInv = pgf_product_list(i)%RhoInv
    1128      972752 :             P = pgf_product_list(i)%P
    1129      972752 :             Q = pgf_product_list(i)%Q
    1130      972752 :             W = pgf_product_list(i)%W
    1131      972752 :             AB = pgf_product_list(i)%AB
    1132      972752 :             CD = pgf_product_list(i)%CD
    1133      243188 :             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      243188 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1139      243188 :             CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize)
    1140      972752 :             DO k = 4, 6
    1141    30973964 :                DO j = 1, mysize
    1142             :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1143             :                                                work_forces(j, k + 3) + &
    1144    30730776 :                                                work_forces(j, k + 6))
    1145             :                END DO
    1146             :             END DO
    1147     3161444 :             DO k = 1, 12
    1148   123166292 :                DO j = 1, mysize
    1149   122923104 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1150             :                END DO
    1151             :             END DO
    1152      243188 :             neris = neris + 12*mysize
    1153      330693 :             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     1137565 :          DO n = 1, 12
    1169             :             tmp_max = 0.0_dp
    1170    44703444 :             DO i = 1, mysize
    1171    44703444 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1172             :             END DO
    1173     1050060 :             tmp_max = tmp_max*max_contraction
    1174     1050060 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1175             : 
    1176     4605193 :             DO l = 1, ncod
    1177     3467628 :                p1 = (l - 1)*ncoc
    1178     7985316 :                DO k = 1, ncoc
    1179     3467628 :                   p2 = (p1 + k - 1)*ncob
    1180    27741024 :                   DO j = 1, ncob
    1181    20805768 :                      p3 = (p2 + j - 1)*ncoa
    1182    67926780 :                      DO i = 1, ncoa
    1183    43653384 :                         p4 = p3 + i
    1184    64459152 :                         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    32025644 :          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    31938139 :       IF (.NOT. use_virial) THEN
    1218    30872855 :          IF (tmp_max_all < eps_schwarz) RETURN
    1219             :       END IF
    1220             : 
    1221    26559621 :       IF (tmp_max_all >= eps_schwarz) THEN
    1222   343115110 :          DO i = 1, 12
    1223 16345484880 :             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   316721640 :                           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 16371878350 :                           s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
    1241             :          END DO
    1242             :       END IF
    1243             : 
    1244    26559621 :       IF (use_virial .AND. tmp_max_all_virial >= eps_schwarz) THEN
    1245    11813685 :          DO i = 1, 12
    1246    44528505 :             DO m = 1, 3
    1247  4937303412 :                primitives_tmp_virial(:, :, :, :) = 0.0_dp
    1248             :                CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
    1249             :                              n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2_virial(1, 1, 1, 1, i, m), &
    1250             :                              sphi_a, &
    1251             :                              sphi_b, &
    1252             :                              sphi_c, &
    1253             :                              sphi_d, &
    1254             :                              primitives_tmp_virial(1, 1, 1, 1), &
    1255    32714820 :                              buffer1, buffer2)
    1256             : 
    1257             :                primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1258             :                                  s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1259             :                                  s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1260             :                                  s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) = &
    1261             :                   primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
    1262             :                                     s_offset_b + 1:s_offset_b + nsob*nl_b, &
    1263             :                                     s_offset_c + 1:s_offset_c + nsoc*nl_c, &
    1264  4948208352 :                                     s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) + primitives_tmp_virial(:, :, :, :)
    1265             :             END DO
    1266             :          END DO
    1267             :       END IF
    1268             : 
    1269             :    END SUBROUTINE evaluate_deriv_eri
    1270             : 
    1271             : ! **************************************************************************************************
    1272             : !> \brief Evaluates max-abs values of  electron repulsion integrals for a primitive quartet
    1273             : !> \param libint ...
    1274             : !> \param A ...
    1275             : !> \param B ...
    1276             : !> \param C ...
    1277             : !> \param D ...
    1278             : !> \param Zeta_A ...
    1279             : !> \param Zeta_B ...
    1280             : !> \param Zeta_C ...
    1281             : !> \param Zeta_D ...
    1282             : !> \param n_a ...
    1283             : !> \param n_b ...
    1284             : !> \param n_c ...
    1285             : !> \param n_d ...
    1286             : !> \param max_val ...
    1287             : !> \param potential_parameter ...
    1288             : !> \param R1 ...
    1289             : !> \param R2 ...
    1290             : !> \param p_work ...
    1291             : !> \par History
    1292             : !>      03.2007 created [Manuel Guidon]
    1293             : !>      08.2007 refactured permutation part [Manuel Guidon]
    1294             : !> \author Manuel Guidon
    1295             : ! **************************************************************************************************
    1296    82517808 :    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    82517808 :       m_max = n_a + n_b + n_c + n_d
    1313    82517808 :       mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
    1314   165035616 :       a_mysize = mysize
    1315             : 
    1316    82517808 :       IF (m_max /= 0) THEN
    1317    51785124 :          perm_case = 1
    1318    51785124 :          IF (n_a < n_b) THEN
    1319    21035472 :             perm_case = perm_case + 1
    1320             :          END IF
    1321    51785124 :          IF (n_c < n_d) THEN
    1322    21035472 :             perm_case = perm_case + 2
    1323             :          END IF
    1324    51785124 :          IF (n_a + n_b > n_c + n_d) THEN
    1325           0 :             perm_case = perm_case + 4
    1326             :          END IF
    1327             : 
    1328    30749652 :          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    30749652 :                                            potential_parameter, libint, R1, R2)
    1332             : 
    1333    30749652 :             CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1334  2644645812 :             DO i = 1, mysize
    1335  2644645812 :                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    21035472 :                                            potential_parameter, libint, R1, R2)
    1356             : 
    1357    21035472 :             CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1358   950932776 :             DO i = 1, mysize
    1359   950932776 :                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    51785124 :             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    30732684 :                                         potential_parameter, libint, R1, R2)
    1397    30732684 :          max_val = ABS(get_ssss_f_val(libint))
    1398             :       END IF
    1399             : 
    1400    82517808 :    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   172605988 :    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   172605988 :                            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   172605988 :                            sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
    1463   172605988 :                            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   172605988 :       m_max = n_a + n_b + n_c + n_d
    1496   172605988 :       mysize = ncoa*ncob*ncoc*ncod
    1497   345211976 :       a_mysize = mysize
    1498             : 
    1499  3527797996 :       work = 0.0_dp
    1500   172605988 :       IF (m_max /= 0) THEN
    1501   136306924 :          perm_case = 1
    1502   136306924 :          IF (n_a < n_b) THEN
    1503    33214955 :             perm_case = perm_case + 1
    1504             :          END IF
    1505   136306924 :          IF (n_c < n_d) THEN
    1506    39690806 :             perm_case = perm_case + 2
    1507             :          END IF
    1508   136306924 :          IF (n_a + n_b > n_c + n_d) THEN
    1509    48257853 :             perm_case = perm_case + 4
    1510             :          END IF
    1511             :          SELECT CASE (perm_case)
    1512             :          CASE (1)
    1513   150016063 :             DO i = 1, nproducts
    1514   415763656 :                A = pgf_product_list(i)%ra
    1515   415763656 :                B = pgf_product_list(i)%rb
    1516   415763656 :                C = pgf_product_list(i)%rc
    1517   415763656 :                D = pgf_product_list(i)%rd
    1518   103940914 :                Rho = pgf_product_list(i)%Rho
    1519   103940914 :                RhoInv = pgf_product_list(i)%RhoInv
    1520   415763656 :                P = pgf_product_list(i)%P
    1521   415763656 :                Q = pgf_product_list(i)%Q
    1522   415763656 :                W = pgf_product_list(i)%W
    1523   415763656 :                AB = pgf_product_list(i)%AB
    1524   415763656 :                CD = pgf_product_list(i)%CD
    1525   103940914 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1526   416788524 :                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   103940914 :                                        P, Q, W, m_max, F)
    1530   103940914 :                CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1531  2046640882 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1532   150016063 :                neris = neris + mysize
    1533             :             END DO
    1534  1101974978 :             DO i = 1, mysize
    1535  1101974978 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1536             :             END DO
    1537    46075149 :             tmp_max = tmp_max*max_contraction
    1538    46075149 :             IF (tmp_max < eps_schwarz) THEN
    1539    44235841 :                RETURN
    1540             :             END IF
    1541             : 
    1542   104859744 :             DO i = 1, ncoa
    1543    68944772 :                p1 = (i - 1)*ncob
    1544   192637436 :                DO j = 1, ncob
    1545    87777692 :                   p2 = (p1 + j - 1)*ncoc
    1546   509540809 :                   DO k = 1, ncoc
    1547   352818345 :                      p3 = (p2 + k - 1)*ncod
    1548  1341415852 :                      DO l = 1, ncod
    1549   900819815 :                         p4 = p3 + l
    1550  1253638160 :                         work2(i, j, k, l) = work(p4)
    1551             :                      END DO
    1552             :                   END DO
    1553             :                END DO
    1554             :             END DO
    1555             :          CASE (2)
    1556    28740527 :             DO i = 1, nproducts
    1557    75922928 :                A = pgf_product_list(i)%ra
    1558    75922928 :                B = pgf_product_list(i)%rb
    1559    75922928 :                C = pgf_product_list(i)%rc
    1560    75922928 :                D = pgf_product_list(i)%rd
    1561    18980732 :                Rho = pgf_product_list(i)%Rho
    1562    18980732 :                RhoInv = pgf_product_list(i)%RhoInv
    1563    75922928 :                P = pgf_product_list(i)%P
    1564    75922928 :                Q = pgf_product_list(i)%Q
    1565    75922928 :                W = pgf_product_list(i)%W
    1566    75922928 :                AB = pgf_product_list(i)%AB
    1567    75922928 :                CD = pgf_product_list(i)%CD
    1568    18980732 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1569    90994066 :                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    18980732 :                                        P, Q, W, m_max, F)
    1574             : 
    1575    18980732 :                CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
    1576   629223435 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1577    28740527 :                neris = neris + mysize
    1578             :             END DO
    1579   411638325 :             DO i = 1, mysize
    1580   411638325 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1581             :             END DO
    1582     9759795 :             tmp_max = tmp_max*max_contraction
    1583     9759795 :             IF (tmp_max < eps_schwarz) THEN
    1584             :                RETURN
    1585             :             END IF
    1586             : 
    1587    31190339 :             DO j = 1, ncob
    1588    24322118 :                p1 = (j - 1)*ncoa
    1589    58440227 :                DO i = 1, ncoa
    1590    27249888 :                   p2 = (p1 + i - 1)*ncoc
    1591   168460530 :                   DO k = 1, ncoc
    1592   116888524 :                      p3 = (p2 + k - 1)*ncod
    1593   461540896 :                      DO l = 1, ncod
    1594   317402484 :                         p4 = p3 + l
    1595   434291008 :                         work2(i, j, k, l) = work(p4)
    1596             :                      END DO
    1597             :                   END DO
    1598             :                END DO
    1599             :             END DO
    1600             :          CASE (3)
    1601    72483411 :             DO i = 1, nproducts
    1602   189223676 :                A = pgf_product_list(i)%ra
    1603   189223676 :                B = pgf_product_list(i)%rb
    1604   189223676 :                C = pgf_product_list(i)%rc
    1605   189223676 :                D = pgf_product_list(i)%rd
    1606    47305919 :                Rho = pgf_product_list(i)%Rho
    1607    47305919 :                RhoInv = pgf_product_list(i)%RhoInv
    1608   189223676 :                P = pgf_product_list(i)%P
    1609   189223676 :                Q = pgf_product_list(i)%Q
    1610   189223676 :                W = pgf_product_list(i)%W
    1611   189223676 :                AB = pgf_product_list(i)%AB
    1612   189223676 :                CD = pgf_product_list(i)%CD
    1613    47305919 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1614   179743752 :                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    47305919 :                                        P, Q, W, m_max, F)
    1619             : 
    1620    47305919 :                CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
    1621   712703048 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1622    72483411 :                neris = neris + mysize
    1623             :             END DO
    1624   449364015 :             DO i = 1, mysize
    1625   449364015 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1626             :             END DO
    1627    25177492 :             tmp_max = tmp_max*max_contraction
    1628    25177492 :             IF (tmp_max < eps_schwarz) THEN
    1629             :                RETURN
    1630             :             END IF
    1631             : 
    1632    52852929 :             DO i = 1, ncoa
    1633    35176039 :                p1 = (i - 1)*ncob
    1634    95452180 :                DO j = 1, ncob
    1635    42599251 :                   p2 = (p1 + j - 1)*ncod
    1636   274730052 :                   DO l = 1, ncod
    1637   196954762 :                      p3 = (p2 + l - 1)*ncoc
    1638   575132773 :                      DO k = 1, ncoc
    1639   335578760 :                         p4 = p3 + k
    1640   532533522 :                         work2(i, j, k, l) = work(p4)
    1641             :                      END DO
    1642             :                   END DO
    1643             :                END DO
    1644             :             END DO
    1645             :          CASE (4)
    1646    18757051 :             DO i = 1, nproducts
    1647    46881664 :                A = pgf_product_list(i)%ra
    1648    46881664 :                B = pgf_product_list(i)%rb
    1649    46881664 :                C = pgf_product_list(i)%rc
    1650    46881664 :                D = pgf_product_list(i)%rd
    1651    11720416 :                Rho = pgf_product_list(i)%Rho
    1652    11720416 :                RhoInv = pgf_product_list(i)%RhoInv
    1653    46881664 :                P = pgf_product_list(i)%P
    1654    46881664 :                Q = pgf_product_list(i)%Q
    1655    46881664 :                W = pgf_product_list(i)%W
    1656    46881664 :                AB = pgf_product_list(i)%AB
    1657    46881664 :                CD = pgf_product_list(i)%CD
    1658    11720416 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1659    53672976 :                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    11720416 :                                        P, Q, W, m_max, F)
    1664             : 
    1665    11720416 :                CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1666   326471247 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1667    18757051 :                neris = neris + mysize
    1668             :             END DO
    1669   242107713 :             DO i = 1, mysize
    1670   242107713 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1671             :             END DO
    1672     7036635 :             tmp_max = tmp_max*max_contraction
    1673     7036635 :             IF (tmp_max < eps_schwarz) THEN
    1674             :                RETURN
    1675             :             END IF
    1676             : 
    1677    23213462 :             DO j = 1, ncob
    1678    18069976 :                p1 = (j - 1)*ncoa
    1679    44075512 :                DO i = 1, ncoa
    1680    20862050 :                   p2 = (p1 + i - 1)*ncod
    1681   137965417 :                   DO l = 1, ncod
    1682    99033391 :                      p3 = (p2 + l - 1)*ncoc
    1683   310435952 :                      DO k = 1, ncoc
    1684   190540511 :                         p4 = p3 + k
    1685   289573902 :                         work2(i, j, k, l) = work(p4)
    1686             :                      END DO
    1687             :                   END DO
    1688             :                END DO
    1689             :             END DO
    1690             :          CASE (5)
    1691    76904773 :             DO i = 1, nproducts
    1692   201480340 :                A = pgf_product_list(i)%ra
    1693   201480340 :                B = pgf_product_list(i)%rb
    1694   201480340 :                C = pgf_product_list(i)%rc
    1695   201480340 :                D = pgf_product_list(i)%rd
    1696    50370085 :                Rho = pgf_product_list(i)%Rho
    1697    50370085 :                RhoInv = pgf_product_list(i)%RhoInv
    1698   201480340 :                P = pgf_product_list(i)%P
    1699   201480340 :                Q = pgf_product_list(i)%Q
    1700   201480340 :                W = pgf_product_list(i)%W
    1701   201480340 :                AB = pgf_product_list(i)%AB
    1702   201480340 :                CD = pgf_product_list(i)%CD
    1703    50370085 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1704   194442227 :                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    50370085 :                                        Q, P, W, m_max, F)
    1709             : 
    1710    50370085 :                CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
    1711   888433702 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1712    76904773 :                neris = neris + mysize
    1713             :             END DO
    1714   566816676 :             DO i = 1, mysize
    1715   566816676 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1716             :             END DO
    1717    26534688 :             tmp_max = tmp_max*max_contraction
    1718    26534688 :             IF (tmp_max < eps_schwarz) THEN
    1719             :                RETURN
    1720             :             END IF
    1721             : 
    1722    46969200 :             DO k = 1, ncoc
    1723    27919827 :                p1 = (k - 1)*ncod
    1724    79107199 :                DO l = 1, ncod
    1725    32137999 :                   p2 = (p1 + l - 1)*ncoa
    1726   206598983 :                   DO i = 1, ncoa
    1727   146541157 :                      p3 = (p2 + i - 1)*ncob
    1728   609450535 :                      DO j = 1, ncob
    1729   430771379 :                         p4 = p3 + j
    1730   577312536 :                         work2(i, j, k, l) = work(p4)
    1731             :                      END DO
    1732             :                   END DO
    1733             :                END DO
    1734             :             END DO
    1735             :          CASE (6)
    1736    36182019 :             DO i = 1, nproducts
    1737    87742132 :                A = pgf_product_list(i)%ra
    1738    87742132 :                B = pgf_product_list(i)%rb
    1739    87742132 :                C = pgf_product_list(i)%rc
    1740    87742132 :                D = pgf_product_list(i)%rd
    1741    21935533 :                Rho = pgf_product_list(i)%Rho
    1742    21935533 :                RhoInv = pgf_product_list(i)%RhoInv
    1743    87742132 :                P = pgf_product_list(i)%P
    1744    87742132 :                Q = pgf_product_list(i)%Q
    1745    87742132 :                W = pgf_product_list(i)%W
    1746    87742132 :                AB = pgf_product_list(i)%AB
    1747    87742132 :                CD = pgf_product_list(i)%CD
    1748    21935533 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1749    79668318 :                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    21935533 :                                        Q, P, W, m_max, F)
    1754             : 
    1755    21935533 :                CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
    1756   314722314 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1757    36182019 :                neris = neris + mysize
    1758             :             END DO
    1759   224105568 :             DO i = 1, mysize
    1760   224105568 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1761             :             END DO
    1762    14246486 :             tmp_max = tmp_max*max_contraction
    1763    14246486 :             IF (tmp_max < eps_schwarz) THEN
    1764             :                RETURN
    1765             :             END IF
    1766             : 
    1767    21498247 :             DO k = 1, ncoc
    1768    12526655 :                p1 = (k - 1)*ncod
    1769    36107298 :                DO l = 1, ncod
    1770    14609051 :                   p2 = (p1 + l - 1)*ncob
    1771   100245727 :                   DO j = 1, ncob
    1772    73110021 :                      p3 = (p2 + j - 1)*ncoa
    1773   240447013 :                      DO i = 1, ncoa
    1774   152727941 :                         p4 = p3 + i
    1775   225837962 :                         work2(i, j, k, l) = work(p4)
    1776             :                      END DO
    1777             :                   END DO
    1778             :                END DO
    1779             :             END DO
    1780             :          CASE (7)
    1781    14675261 :             DO i = 1, nproducts
    1782    37482484 :                A = pgf_product_list(i)%ra
    1783    37482484 :                B = pgf_product_list(i)%rb
    1784    37482484 :                C = pgf_product_list(i)%rc
    1785    37482484 :                D = pgf_product_list(i)%rd
    1786     9370621 :                Rho = pgf_product_list(i)%Rho
    1787     9370621 :                RhoInv = pgf_product_list(i)%RhoInv
    1788    37482484 :                P = pgf_product_list(i)%P
    1789    37482484 :                Q = pgf_product_list(i)%Q
    1790    37482484 :                W = pgf_product_list(i)%W
    1791    37482484 :                AB = pgf_product_list(i)%AB
    1792    37482484 :                CD = pgf_product_list(i)%CD
    1793     9370621 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1794    50651904 :                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     9370621 :                                        Q, P, W, m_max, F)
    1799             : 
    1800     9370621 :                CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
    1801   483093315 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1802    14675261 :                neris = neris + mysize
    1803             :             END DO
    1804   347617502 :             DO i = 1, mysize
    1805   347617502 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1806             :             END DO
    1807     5304640 :             tmp_max = tmp_max*max_contraction
    1808     5304640 :             IF (tmp_max < eps_schwarz) THEN
    1809             :                RETURN
    1810             :             END IF
    1811             : 
    1812    16285965 :             DO l = 1, ncod
    1813    12632518 :                p1 = (l - 1)*ncoc
    1814    30585671 :                DO k = 1, ncoc
    1815    14299706 :                   p2 = (p1 + k - 1)*ncoa
    1816   101393391 :                   DO i = 1, ncoa
    1817    74461167 :                      p3 = (p2 + i - 1)*ncob
    1818   362285676 :                      DO j = 1, ncob
    1819   273524803 :                         p4 = p3 + j
    1820   347985970 :                         work2(i, j, k, l) = work(p4)
    1821             :                      END DO
    1822             :                   END DO
    1823             :                END DO
    1824             :             END DO
    1825             :          CASE (8)
    1826     4950060 :             DO i = 1, nproducts
    1827    11112084 :                A = pgf_product_list(i)%ra
    1828    11112084 :                B = pgf_product_list(i)%rb
    1829    11112084 :                C = pgf_product_list(i)%rc
    1830    11112084 :                D = pgf_product_list(i)%rd
    1831     2778021 :                Rho = pgf_product_list(i)%Rho
    1832     2778021 :                RhoInv = pgf_product_list(i)%RhoInv
    1833    11112084 :                P = pgf_product_list(i)%P
    1834    11112084 :                Q = pgf_product_list(i)%Q
    1835    11112084 :                W = pgf_product_list(i)%W
    1836    11112084 :                AB = pgf_product_list(i)%AB
    1837    11112084 :                CD = pgf_product_list(i)%CD
    1838     2778021 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1839    15792571 :                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     2778021 :                                        Q, P, W, m_max, F)
    1844             : 
    1845     2778021 :                CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
    1846   136032195 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1847     4950060 :                neris = neris + mysize
    1848             :             END DO
    1849   111575091 :             DO i = 1, mysize
    1850   111575091 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1851             :             END DO
    1852     2172039 :             tmp_max = tmp_max*max_contraction
    1853     2172039 :             IF (tmp_max < eps_schwarz) THEN
    1854             :                RETURN
    1855             :             END IF
    1856             : 
    1857   143698674 :             DO l = 1, ncod
    1858     5769988 :                p1 = (l - 1)*ncoc
    1859    13222058 :                DO k = 1, ncoc
    1860     5830308 :                   p2 = (p1 + k - 1)*ncob
    1861    47652028 :                   DO j = 1, ncob
    1862    36051732 :                      p3 = (p2 + j - 1)*ncoa
    1863   126375894 :                      DO i = 1, ncoa
    1864    84493854 :                         p4 = p3 + i
    1865   120545586 :                         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   109797931 :          DO i = 1, nproducts
    1873   293995468 :             A = pgf_product_list(i)%ra
    1874   293995468 :             B = pgf_product_list(i)%rb
    1875   293995468 :             C = pgf_product_list(i)%rc
    1876   293995468 :             D = pgf_product_list(i)%rd
    1877    73498867 :             Rho = pgf_product_list(i)%Rho
    1878    73498867 :             RhoInv = pgf_product_list(i)%RhoInv
    1879   293995468 :             P = pgf_product_list(i)%P
    1880   293995468 :             Q = pgf_product_list(i)%Q
    1881   293995468 :             W = pgf_product_list(i)%W
    1882    73498867 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1883   146997734 :             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    73498867 :                                     P, Q, W, m_max, F)
    1888    73498867 :             work(1) = work(1) + F(1)
    1889   109797931 :             neris = neris + mysize
    1890             :          END DO
    1891    36299064 :          work2(1, 1, 1, 1) = work(1)
    1892    36299064 :          tmp_max = max_contraction*ABS(work(1))
    1893    36299064 :          IF (tmp_max < eps_schwarz) RETURN
    1894             :       END IF
    1895             : 
    1896   128370147 :       IF (tmp_max < eps_schwarz) RETURN
    1897  7965812003 :       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   128370147 :                     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  7965812003 :                     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   339901108 :    SUBROUTINE build_quartet_data(libint, A, B, C, D, &
    1940             :                                  ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1941   339901108 :                                  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   339901108 :       CALL cp_libint_set_params_eri(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, P, Q, W, m_max, F)
    1951   339901108 :    END SUBROUTINE build_quartet_data
    1952             : 
    1953             : END MODULE hfx_libint_interface

Generated by: LCOV version 1.15