LCOV - code coverage report
Current view: top level - src - hfx_libint_interface.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:936074a) Lines: 96.7 % 983 951
Test Date: 2025-12-04 06:27:48 Functions: 100.0 % 6 6

            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     83698296 :    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     83698296 :       Zeta = Zeta_A + Zeta_B
      98     83698296 :       ZetaInv = 1.0_dp/Zeta
      99     83698296 :       Eta = Zeta_C + Zeta_D
     100     83698296 :       EtaInv = 1.0_dp/Eta
     101     83698296 :       ZetapEtaInv = Zeta + Eta
     102     83698296 :       ZetapEtaInv = 1.0_dp/ZetapEtaInv
     103     83698296 :       Rho = Zeta*Eta*ZetapEtaInv
     104     83698296 :       RhoInv = 1.0_dp/Rho
     105              : 
     106    334793184 :       DO i = 1, 3
     107    251094888 :          P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv
     108    251094888 :          Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv
     109    251094888 :          AB(i) = A(i) - B(i)
     110    251094888 :          CD(i) = C(i) - D(i)
     111    251094888 :          PQ(i) = P(i) - Q(i)
     112    334793184 :          W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv
     113              :       END DO
     114              : 
     115    334793184 :       AB2 = DOT_PRODUCT(AB, AB)
     116    334793184 :       CD2 = DOT_PRODUCT(CD, CD)
     117    334793184 :       PQ2 = DOT_PRODUCT(PQ, PQ)
     118              : 
     119     83698296 :       S1234 = EXP((-Zeta_A*Zeta_B*ZetaInv*AB2) + (-Zeta_C*Zeta_D*EtaInv*CD2))
     120     83698296 :       T = Rho*PQ2
     121              : 
     122     96904248 :       SELECT CASE (potential_parameter%potential_type)
     123              :       CASE (do_potential_truncated)
     124     13205952 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     125     13205952 :          R1 = R11
     126     13205952 :          R2 = R22
     127     13205952 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     128            0 :             RETURN
     129              :          END IF
     130     13205952 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     131     13205952 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     132     13205952 :          factor = 2.0_dp*Pi*RhoInv
     133              :       CASE (do_potential_coulomb)
     134     63390024 :          CALL fgamma(m_max, T, F(1))
     135     63390024 :          factor = 2.0_dp*Pi*RhoInv
     136              :       CASE (do_potential_short)
     137      3091812 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     138      3091812 :          R1 = R11
     139      3091812 :          R2 = R22
     140      3091812 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     141              :             RETURN
     142              :          END IF
     143      3091812 :          CALL fgamma(m_max, T, F(1))
     144      3091812 :          omega2 = potential_parameter%omega**2
     145      3091812 :          omega_corr2 = omega2/(omega2 + Rho)
     146      3091812 :          omega_corr = SQRT(omega_corr2)
     147      3091812 :          T = T*omega_corr2
     148      3091812 :          CALL fgamma(m_max, T, Fm)
     149      3091812 :          tmp = -omega_corr
     150     11625504 :          DO i = 1, m_max + 1
     151      8533692 :             F(i) = F(i) + Fm(i)*tmp
     152     11625504 :             tmp = tmp*omega_corr2
     153              :          END DO
     154      3091812 :          factor = 2.0_dp*Pi*RhoInv
     155              :       CASE (do_potential_long)
     156       978084 :          omega2 = potential_parameter%omega**2
     157       978084 :          omega_corr2 = omega2/(omega2 + Rho)
     158       978084 :          omega_corr = SQRT(omega_corr2)
     159       978084 :          T = T*omega_corr2
     160       978084 :          CALL fgamma(m_max, T, F(1))
     161       978084 :          tmp = omega_corr
     162      3565704 :          DO i = 1, m_max + 1
     163      2587620 :             F(i) = F(i)*tmp
     164      3565704 :             tmp = tmp*omega_corr2
     165              :          END DO
     166       978084 :          factor = 2.0_dp*Pi*RhoInv
     167              :       CASE (do_potential_mix_cl)
     168      1558026 :          CALL fgamma(m_max, T, F(1))
     169      1558026 :          omega2 = potential_parameter%omega**2
     170      1558026 :          omega_corr2 = omega2/(omega2 + Rho)
     171      1558026 :          omega_corr = SQRT(omega_corr2)
     172      1558026 :          T = T*omega_corr2
     173      1558026 :          CALL fgamma(m_max, T, Fm)
     174      1558026 :          tmp = omega_corr
     175      5416428 :          DO i = 1, m_max + 1
     176      3858402 :             F(i) = F(i)*potential_parameter%scale_coulomb + Fm(i)*tmp*potential_parameter%scale_longrange
     177      5416428 :             tmp = tmp*omega_corr2
     178              :          END DO
     179      1558026 :          factor = 2.0_dp*Pi*RhoInv
     180              :       CASE (do_potential_mix_cl_trunc)
     181              : 
     182              :          ! truncated
     183       959298 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     184       959298 :          R1 = R11
     185       959298 :          R2 = R22
     186       959298 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     187              :             RETURN
     188              :          END IF
     189       959298 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     190       959298 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     191              : 
     192              :          ! Coulomb
     193       959298 :          CALL fgamma(m_max, T, Fm)
     194              : 
     195      3419052 :          DO i = 1, m_max + 1
     196              :             F(i) = F(i)*(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
     197      3419052 :                    Fm(i)*potential_parameter%scale_longrange
     198              :          END DO
     199              : 
     200              :          ! longrange
     201       959298 :          omega2 = potential_parameter%omega**2
     202       959298 :          omega_corr2 = omega2/(omega2 + Rho)
     203       959298 :          omega_corr = SQRT(omega_corr2)
     204       959298 :          T = T*omega_corr2
     205       959298 :          CALL fgamma(m_max, T, Fm)
     206       959298 :          tmp = omega_corr
     207      3419052 :          DO i = 1, m_max + 1
     208      2459754 :             F(i) = F(i) + Fm(i)*tmp*potential_parameter%scale_longrange
     209      3419052 :             tmp = tmp*omega_corr2
     210              :          END DO
     211       959298 :          factor = 2.0_dp*Pi*RhoInv
     212              : 
     213              :       CASE (do_potential_gaussian)
     214            0 :          omega2 = potential_parameter%omega**2
     215            0 :          T = -omega2*T/(Rho + omega2)
     216            0 :          tmp = 1.0_dp
     217            0 :          DO i = 1, m_max + 1
     218            0 :             F(i) = EXP(T)*tmp
     219            0 :             tmp = tmp*omega2/(Rho + omega2)
     220              :          END DO
     221            0 :          factor = (Pi/(Rho + omega2))**(1.5_dp)
     222              :       CASE (do_potential_mix_lg)
     223       272700 :          omega2 = potential_parameter%omega**2
     224       272700 :          omega_corr2 = omega2/(omega2 + Rho)
     225       272700 :          omega_corr = SQRT(omega_corr2)
     226       272700 :          T = T*omega_corr2
     227       272700 :          CALL fgamma(m_max, T, Fm)
     228       272700 :          tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange
     229       981720 :          DO i = 1, m_max + 1
     230       709020 :             Fm(i) = Fm(i)*tmp
     231       981720 :             tmp = tmp*omega_corr2
     232              :          END DO
     233              :          T = Rho*PQ2
     234       272700 :          T = -omega2*T/(Rho + omega2)
     235       272700 :          tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
     236       981720 :          DO i = 1, m_max + 1
     237       709020 :             F(i) = EXP(T)*tmp + Fm(i)
     238       981720 :             tmp = tmp*omega2/(Rho + omega2)
     239              :          END DO
     240       242400 :          factor = 1.0_dp
     241              :       CASE (do_potential_id)
     242       242400 :          ssss = -Zeta_A*Zeta_B*ZetaInv*AB2
     243              : 
     244       242400 :          num = (Zeta_A + Zeta_B)*Zeta_C
     245       242400 :          den = Zeta_A + Zeta_B + Zeta_C
     246       969600 :          ssss = ssss - num/den*SUM((P - C)**2)
     247              : 
     248       969600 :          G(:) = (Zeta_A*A(:) + Zeta_B*B(:) + Zeta_C*C(:))/den
     249       242400 :          num = den*Zeta_D
     250       242400 :          den = den + Zeta_D
     251       969600 :          ssss = ssss - num/den*SUM((G - D)**2)
     252              : 
     253      5332800 :          F(:) = EXP(ssss)
     254       242400 :          factor = 1.0_dp
     255     98105946 :          IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
     256              :       END SELECT
     257              : 
     258     83698296 :       tmp = (Pi*ZetapEtaInv)**3
     259     83698296 :       factor = factor*S1234*SQRT(tmp)
     260              : 
     261    323189496 :       DO i = 1, m_max + 1
     262    323189496 :          F(i) = F(i)*factor
     263              :       END DO
     264              : 
     265     83698296 :       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     99205823 :    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     99205823 :                                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     99205823 :                                           ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F)
     313              : 
     314     99205823 :    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     32373636 :    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     32373636 :                                  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     32373636 :                                  sphi_a, sphi_b, sphi_c, sphi_d, &
     388     32373636 :                                  work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
     389     32373636 :                                  use_virial, work_virial, work2_virial, primitives_tmp_virial, &
     390     32373636 :                                  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     32373636 :       m_max = n_a + n_b + n_c + n_d
     433     32373636 :       m_max = m_max + 1
     434     32373636 :       mysize = ncoa*ncob*ncoc*ncod
     435     64747272 :       a_mysize = mysize
     436              : 
     437   4310902044 :       work = 0.0_dp
     438   8274148116 :       work2 = 0.0_dp
     439              : 
     440     32373636 :       IF (use_virial) THEN
     441    850401248 :          work_virial = 0.0_dp
     442   1465200824 :          work2_virial = 0.0_dp
     443              :       END IF
     444              : 
     445     32373636 :       perm_case = 1
     446     32373636 :       IF (n_a < n_b) THEN
     447      3012902 :          perm_case = perm_case + 1
     448              :       END IF
     449     32373636 :       IF (n_c < n_d) THEN
     450      5168834 :          perm_case = perm_case + 2
     451              :       END IF
     452     32373636 :       IF (n_a + n_b > n_c + n_d) THEN
     453      5872900 :          perm_case = perm_case + 4
     454              :       END IF
     455              : 
     456              :       SELECT CASE (perm_case)
     457              :       CASE (1)
     458     81269253 :          DO i = 1, nproducts
     459    241892468 :             A = pgf_product_list(i)%ra
     460    241892468 :             B = pgf_product_list(i)%rb
     461    241892468 :             C = pgf_product_list(i)%rc
     462    241892468 :             D = pgf_product_list(i)%rd
     463     60473117 :             Rho = pgf_product_list(i)%Rho
     464     60473117 :             RhoInv = pgf_product_list(i)%RhoInv
     465    241892468 :             P = pgf_product_list(i)%P
     466    241892468 :             Q = pgf_product_list(i)%Q
     467    241892468 :             W = pgf_product_list(i)%W
     468    241892468 :             AB = pgf_product_list(i)%AB
     469    241892468 :             CD = pgf_product_list(i)%CD
     470     60473117 :             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     60473117 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     476     60473117 :             CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize)
     477    241892468 :             DO k = 4, 6
     478   1868183360 :                DO j = 1, mysize
     479              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     480              :                                                work_forces(j, k + 3) + &
     481   1807710243 :                                                work_forces(j, k + 6))
     482              :                END DO
     483              :             END DO
     484    786150521 :             DO k = 1, 12
     485   7291314089 :                DO j = 1, mysize
     486   7230840972 :                   work(j, k) = work(j, k) + work_forces(j, k)
     487              :                END DO
     488              :             END DO
     489     60473117 :             neris = neris + 12*mysize
     490     81269253 :             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    270349768 :          DO n = 1, 12
     506              :             tmp_max = 0.0_dp
     507   2139902736 :             DO i = 1, mysize
     508   2139902736 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     509              :             END DO
     510    249553632 :             tmp_max = tmp_max*max_contraction
     511    249553632 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     512              : 
     513    615242152 :             DO i = 1, ncoa
     514    344892384 :                p1 = (i - 1)*ncob
     515    974615160 :                DO j = 1, ncob
     516    380169144 :                   p2 = (p1 + j - 1)*ncoc
     517   1749204984 :                   DO k = 1, ncoc
     518   1024143456 :                      p3 = (p2 + k - 1)*ncod
     519   3294661704 :                      DO l = 1, ncod
     520   1890349104 :                         p4 = p3 + l
     521   2914492560 :                         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     20796136 :          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      5773026 :          DO i = 1, nproducts
     554     18297236 :             A = pgf_product_list(i)%ra
     555     18297236 :             B = pgf_product_list(i)%rb
     556     18297236 :             C = pgf_product_list(i)%rc
     557     18297236 :             D = pgf_product_list(i)%rd
     558      4574309 :             Rho = pgf_product_list(i)%Rho
     559      4574309 :             RhoInv = pgf_product_list(i)%RhoInv
     560     18297236 :             P = pgf_product_list(i)%P
     561     18297236 :             Q = pgf_product_list(i)%Q
     562     18297236 :             W = pgf_product_list(i)%W
     563     18297236 :             AB = pgf_product_list(i)%AB
     564     18297236 :             CD = pgf_product_list(i)%CD
     565      4574309 :             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      4574309 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     571      4574309 :             CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize)
     572     18297236 :             DO k = 4, 6
     573    326559881 :                DO j = 1, mysize
     574              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     575              :                                                work_forces(j, k + 3) + &
     576    321985572 :                                                work_forces(j, k + 6))
     577              :                END DO
     578              :             END DO
     579     59466017 :             DO k = 1, 12
     580   1292516597 :                DO j = 1, mysize
     581   1287942288 :                   work(j, k) = work(j, k) + work_forces(j, k)
     582              :                END DO
     583              :             END DO
     584      4574309 :             neris = neris + 12*mysize
     585      5773026 :             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     15583321 :          DO n = 1, 12
     601              :             tmp_max = 0.0_dp
     602    378963552 :             DO i = 1, mysize
     603    378963552 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     604              :             END DO
     605     14384604 :             tmp_max = tmp_max*max_contraction
     606     14384604 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     607              : 
     608     61328881 :             DO j = 1, ncob
     609     45745560 :                p1 = (j - 1)*ncoa
     610    107441580 :                DO i = 1, ncoa
     611     47311416 :                   p2 = (p1 + i - 1)*ncoc
     612    263863188 :                   DO k = 1, ncoc
     613    170806212 :                      p3 = (p2 + k - 1)*ncod
     614    582696576 :                      DO l = 1, ncod
     615    364578948 :                         p4 = p3 + l
     616    535385160 :                         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      1198717 :          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     16732645 :          DO i = 1, nproducts
     649     51545316 :             A = pgf_product_list(i)%ra
     650     51545316 :             B = pgf_product_list(i)%rb
     651     51545316 :             C = pgf_product_list(i)%rc
     652     51545316 :             D = pgf_product_list(i)%rd
     653     12886329 :             Rho = pgf_product_list(i)%Rho
     654     12886329 :             RhoInv = pgf_product_list(i)%RhoInv
     655     51545316 :             P = pgf_product_list(i)%P
     656     51545316 :             Q = pgf_product_list(i)%Q
     657     51545316 :             W = pgf_product_list(i)%W
     658     51545316 :             AB = pgf_product_list(i)%AB
     659     51545316 :             CD = pgf_product_list(i)%CD
     660     12886329 :             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     12886329 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     666     12886329 :             CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize)
     667     51545316 :             DO k = 4, 6
     668    457769523 :                DO j = 1, mysize
     669              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     670              :                                                work_forces(j, k + 3) + &
     671    444883194 :                                                work_forces(j, k + 6))
     672              :                END DO
     673              :             END DO
     674    167522277 :             DO k = 1, 12
     675   1792419105 :                DO j = 1, mysize
     676   1779532776 :                   work(j, k) = work(j, k) + work_forces(j, k)
     677              :                END DO
     678              :             END DO
     679     12886329 :             neris = neris + 12*mysize
     680     16732645 :             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     50002108 :          DO n = 1, 12
     696              :             tmp_max = 0.0_dp
     697    541721964 :             DO i = 1, mysize
     698    541721964 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     699              :             END DO
     700     46155792 :             tmp_max = tmp_max*max_contraction
     701     46155792 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     702              : 
     703    130233748 :             DO i = 1, ncoa
     704     80231640 :                p1 = (i - 1)*ncob
     705    214294920 :                DO j = 1, ncob
     706     87907488 :                   p2 = (p1 + j - 1)*ncod
     707    508848564 :                   DO l = 1, ncod
     708    340709436 :                      p3 = (p2 + l - 1)*ncoc
     709    924183096 :                      DO k = 1, ncoc
     710    495566172 :                         p4 = p3 + k
     711    836275608 :                         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      3846316 :          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      2915097 :          DO i = 1, nproducts
     744      9022120 :             A = pgf_product_list(i)%ra
     745      9022120 :             B = pgf_product_list(i)%rb
     746      9022120 :             C = pgf_product_list(i)%rc
     747      9022120 :             D = pgf_product_list(i)%rd
     748      2255530 :             Rho = pgf_product_list(i)%Rho
     749      2255530 :             RhoInv = pgf_product_list(i)%RhoInv
     750      9022120 :             P = pgf_product_list(i)%P
     751      9022120 :             Q = pgf_product_list(i)%Q
     752      9022120 :             W = pgf_product_list(i)%W
     753      9022120 :             AB = pgf_product_list(i)%AB
     754      9022120 :             CD = pgf_product_list(i)%CD
     755      2255530 :             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      2255530 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     760      2255530 :             CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize)
     761      9022120 :             DO k = 4, 6
     762    123577855 :                DO j = 1, mysize
     763              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     764              :                                                work_forces(j, k + 3) + &
     765    121322325 :                                                work_forces(j, k + 6))
     766              :                END DO
     767              :             END DO
     768     29321890 :             DO k = 1, 12
     769    487544830 :                DO j = 1, mysize
     770    485289300 :                   work(j, k) = work(j, k) + work_forces(j, k)
     771              :                END DO
     772              :             END DO
     773      2255530 :             neris = neris + 12*mysize
     774      2915097 :             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      8574371 :          DO n = 1, 12
     790              :             tmp_max = 0.0_dp
     791    161407320 :             DO i = 1, mysize
     792    161407320 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     793              :             END DO
     794      7914804 :             tmp_max = tmp_max*max_contraction
     795      7914804 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     796              : 
     797     33392303 :             DO j = 1, ncob
     798     24817932 :                p1 = (j - 1)*ncoa
     799     58798572 :                DO i = 1, ncoa
     800     26065836 :                   p2 = (p1 + i - 1)*ncod
     801    151747020 :                   DO l = 1, ncod
     802    100863252 :                      p3 = (p2 + l - 1)*ncoc
     803    280421604 :                      DO k = 1, ncoc
     804    153492516 :                         p4 = p3 + k
     805    254355768 :                         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       659567 :          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     17453778 :          DO i = 1, nproducts
     838     53221648 :             A = pgf_product_list(i)%ra
     839     53221648 :             B = pgf_product_list(i)%rb
     840     53221648 :             C = pgf_product_list(i)%rc
     841     53221648 :             D = pgf_product_list(i)%rd
     842     13305412 :             Rho = pgf_product_list(i)%Rho
     843     13305412 :             RhoInv = pgf_product_list(i)%RhoInv
     844     53221648 :             P = pgf_product_list(i)%P
     845     53221648 :             Q = pgf_product_list(i)%Q
     846     53221648 :             W = pgf_product_list(i)%W
     847     53221648 :             AB = pgf_product_list(i)%AB
     848     53221648 :             CD = pgf_product_list(i)%CD
     849     13305412 :             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     13305412 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     854     13305412 :             CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize)
     855     53221648 :             DO k = 4, 6
     856    511055041 :                DO j = 1, mysize
     857              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     858              :                                                work_forces(j, k + 3) + &
     859    497749629 :                                                work_forces(j, k + 6))
     860              :                END DO
     861              :             END DO
     862    172970356 :             DO k = 1, 12
     863   2004303928 :                DO j = 1, mysize
     864   1990998516 :                   work(j, k) = work(j, k) + work_forces(j, k)
     865              :                END DO
     866              :             END DO
     867     13305412 :             neris = neris + 12*mysize
     868     17453778 :             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     53928758 :          DO n = 1, 12
     884              :             tmp_max = 0.0_dp
     885    601866816 :             DO i = 1, mysize
     886    601866816 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     887              :             END DO
     888     49780392 :             tmp_max = tmp_max*max_contraction
     889     49780392 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     890              : 
     891    121091426 :             DO k = 1, ncoc
     892     67162668 :                p1 = (k - 1)*ncod
     893    187991208 :                DO l = 1, ncod
     894     71048148 :                   p2 = (p1 + l - 1)*ncoa
     895    406166232 :                   DO i = 1, ncoa
     896    267955416 :                      p3 = (p2 + i - 1)*ncob
     897    891089988 :                      DO j = 1, ncob
     898    552086424 :                         p4 = p3 + j
     899    820041840 :                         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      4148366 :          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      4480916 :          DO i = 1, nproducts
     932     13677332 :             A = pgf_product_list(i)%ra
     933     13677332 :             B = pgf_product_list(i)%rb
     934     13677332 :             C = pgf_product_list(i)%rc
     935     13677332 :             D = pgf_product_list(i)%rd
     936      3419333 :             Rho = pgf_product_list(i)%Rho
     937      3419333 :             RhoInv = pgf_product_list(i)%RhoInv
     938     13677332 :             P = pgf_product_list(i)%P
     939     13677332 :             Q = pgf_product_list(i)%Q
     940     13677332 :             W = pgf_product_list(i)%W
     941     13677332 :             AB = pgf_product_list(i)%AB
     942     13677332 :             CD = pgf_product_list(i)%CD
     943      3419333 :             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      3419333 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     949      3419333 :             CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize)
     950     13677332 :             DO k = 4, 6
     951    130273430 :                DO j = 1, mysize
     952              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     953              :                                                work_forces(j, k + 3) + &
     954    126854097 :                                                work_forces(j, k + 6))
     955              :                END DO
     956              :             END DO
     957     44451329 :             DO k = 1, 12
     958    510835721 :                DO j = 1, mysize
     959    507416388 :                   work(j, k) = work(j, k) + work_forces(j, k)
     960              :                END DO
     961              :             END DO
     962      3419333 :             neris = neris + 12*mysize
     963      4480916 :             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     13800579 :          DO n = 1, 12
     979              :             tmp_max = 0.0_dp
     980    169918056 :             DO i = 1, mysize
     981    169918056 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     982              :             END DO
     983     12738996 :             tmp_max = tmp_max*max_contraction
     984     12738996 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     985              : 
     986     29933511 :             DO k = 1, ncoc
     987     16132932 :                p1 = (k - 1)*ncod
     988     47459340 :                DO l = 1, ncod
     989     18587412 :                   p2 = (p1 + l - 1)*ncob
     990    116655372 :                   DO j = 1, ncob
     991     81935028 :                      p3 = (p2 + j - 1)*ncoa
     992    257701500 :                      DO i = 1, ncoa
     993    157179060 :                         p4 = p3 + i
     994    239114088 :                         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      1061583 :          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      2613015 :          DO i = 1, nproducts
    1027      8172396 :             A = pgf_product_list(i)%ra
    1028      8172396 :             B = pgf_product_list(i)%rb
    1029      8172396 :             C = pgf_product_list(i)%rc
    1030      8172396 :             D = pgf_product_list(i)%rd
    1031      2043099 :             Rho = pgf_product_list(i)%Rho
    1032      2043099 :             RhoInv = pgf_product_list(i)%RhoInv
    1033      8172396 :             P = pgf_product_list(i)%P
    1034      8172396 :             Q = pgf_product_list(i)%Q
    1035      8172396 :             W = pgf_product_list(i)%W
    1036      8172396 :             AB = pgf_product_list(i)%AB
    1037      8172396 :             CD = pgf_product_list(i)%CD
    1038      2043099 :             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      2043099 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1044      2043099 :             CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize)
    1045      8172396 :             DO k = 4, 6
    1046    197525367 :                DO j = 1, mysize
    1047              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1048              :                                                work_forces(j, k + 3) + &
    1049    195482268 :                                                work_forces(j, k + 6))
    1050              :                END DO
    1051              :             END DO
    1052     26560287 :             DO k = 1, 12
    1053    783972171 :                DO j = 1, mysize
    1054    781929072 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1055              :                END DO
    1056              :             END DO
    1057      2043099 :             neris = neris + 12*mysize
    1058      2613015 :             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      7408908 :          DO n = 1, 12
    1074              :             tmp_max = 0.0_dp
    1075    237028896 :             DO i = 1, mysize
    1076    237028896 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1077              :             END DO
    1078      6838992 :             tmp_max = tmp_max*max_contraction
    1079      6838992 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1080              : 
    1081     28491480 :             DO l = 1, ncod
    1082     21082572 :                p1 = (l - 1)*ncoc
    1083     49322088 :                DO k = 1, ncoc
    1084     21400524 :                   p2 = (p1 + k - 1)*ncoa
    1085    128660184 :                   DO i = 1, ncoa
    1086     86177088 :                      p3 = (p2 + i - 1)*ncob
    1087    337767516 :                      DO j = 1, ncob
    1088    230189904 :                         p4 = p3 + j
    1089    316366992 :                         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       569916 :          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       341729 :          DO i = 1, nproducts
    1122       994776 :             A = pgf_product_list(i)%ra
    1123       994776 :             B = pgf_product_list(i)%rb
    1124       994776 :             C = pgf_product_list(i)%rc
    1125       994776 :             D = pgf_product_list(i)%rd
    1126       248694 :             Rho = pgf_product_list(i)%Rho
    1127       248694 :             RhoInv = pgf_product_list(i)%RhoInv
    1128       994776 :             P = pgf_product_list(i)%P
    1129       994776 :             Q = pgf_product_list(i)%Q
    1130       994776 :             W = pgf_product_list(i)%W
    1131       994776 :             AB = pgf_product_list(i)%AB
    1132       994776 :             CD = pgf_product_list(i)%CD
    1133       248694 :             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       248694 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1139       248694 :             CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize)
    1140       994776 :             DO k = 4, 6
    1141     31730280 :                DO j = 1, mysize
    1142              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1143              :                                                work_forces(j, k + 3) + &
    1144     31481586 :                                                work_forces(j, k + 6))
    1145              :                END DO
    1146              :             END DO
    1147      3233022 :             DO k = 1, 12
    1148    126175038 :                DO j = 1, mysize
    1149    125926344 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1150              :                END DO
    1151              :             END DO
    1152       248694 :             neris = neris + 12*mysize
    1153       341729 :             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      1209455 :          DO n = 1, 12
    1169              :             tmp_max = 0.0_dp
    1170     47719068 :             DO i = 1, mysize
    1171     47719068 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1172              :             END DO
    1173      1116420 :             tmp_max = tmp_max*max_contraction
    1174      1116420 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1175              : 
    1176      4901651 :             DO l = 1, ncod
    1177      3692196 :                p1 = (l - 1)*ncoc
    1178      8500812 :                DO k = 1, ncoc
    1179      3692196 :                   p2 = (p1 + k - 1)*ncob
    1180     29537568 :                   DO j = 1, ncob
    1181     22153176 :                      p3 = (p2 + j - 1)*ncoa
    1182     72448020 :                      DO i = 1, ncoa
    1183     46602648 :                         p4 = p3 + i
    1184     68755824 :                         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     32466671 :          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     32373636 :       IF (.NOT. use_virial) THEN
    1218     31308352 :          IF (tmp_max_all < eps_schwarz) RETURN
    1219              :       END IF
    1220              : 
    1221     26908466 :       IF (tmp_max_all >= eps_schwarz) THEN
    1222    347651226 :          DO i = 1, 12
    1223  16768732032 :             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    320908824 :                           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  16795474434 :                           s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
    1241              :          END DO
    1242              :       END IF
    1243              : 
    1244     26908466 :       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     83698296 :    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     83698296 :       m_max = n_a + n_b + n_c + n_d
    1313     83698296 :       mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
    1314    167396592 :       a_mysize = mysize
    1315              : 
    1316     83698296 :       IF (m_max /= 0) THEN
    1317     52482024 :          perm_case = 1
    1318     52482024 :          IF (n_a < n_b) THEN
    1319     21319686 :             perm_case = perm_case + 1
    1320              :          END IF
    1321     52482024 :          IF (n_c < n_d) THEN
    1322     21319686 :             perm_case = perm_case + 2
    1323              :          END IF
    1324     52482024 :          IF (n_a + n_b > n_c + n_d) THEN
    1325            0 :             perm_case = perm_case + 4
    1326              :          END IF
    1327              : 
    1328     31162338 :          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     31162338 :                                            potential_parameter, libint, R1, R2)
    1332              : 
    1333     31162338 :             CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1334   2676179022 :             DO i = 1, mysize
    1335   2676179022 :                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     21319686 :                                            potential_parameter, libint, R1, R2)
    1356              : 
    1357     21319686 :             CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1358    963095802 :             DO i = 1, mysize
    1359    963095802 :                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     52482024 :             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     31216272 :                                         potential_parameter, libint, R1, R2)
    1397     31216272 :          max_val = ABS(get_ssss_f_val(libint))
    1398              :       END IF
    1399              : 
    1400     83698296 :    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    173669481 :    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    173669481 :                            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    173669481 :                            sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
    1463    173669481 :                            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    173669481 :       m_max = n_a + n_b + n_c + n_d
    1496    173669481 :       mysize = ncoa*ncob*ncoc*ncod
    1497    347338962 :       a_mysize = mysize
    1498              : 
    1499   3556031927 :       work = 0.0_dp
    1500    173669481 :       IF (m_max /= 0) THEN
    1501    137205406 :          perm_case = 1
    1502    137205406 :          IF (n_a < n_b) THEN
    1503     33438725 :             perm_case = perm_case + 1
    1504              :          END IF
    1505    137205406 :          IF (n_c < n_d) THEN
    1506     39958092 :             perm_case = perm_case + 2
    1507              :          END IF
    1508    137205406 :          IF (n_a + n_b > n_c + n_d) THEN
    1509     48567356 :             perm_case = perm_case + 4
    1510              :          END IF
    1511              :          SELECT CASE (perm_case)
    1512              :          CASE (1)
    1513    198471447 :             DO i = 1, nproducts
    1514    608389112 :                A = pgf_product_list(i)%ra
    1515    608389112 :                B = pgf_product_list(i)%rb
    1516    608389112 :                C = pgf_product_list(i)%rc
    1517    608389112 :                D = pgf_product_list(i)%rd
    1518    152097278 :                Rho = pgf_product_list(i)%Rho
    1519    152097278 :                RhoInv = pgf_product_list(i)%RhoInv
    1520    608389112 :                P = pgf_product_list(i)%P
    1521    608389112 :                Q = pgf_product_list(i)%Q
    1522    608389112 :                W = pgf_product_list(i)%W
    1523    608389112 :                AB = pgf_product_list(i)%AB
    1524    608389112 :                CD = pgf_product_list(i)%CD
    1525    152097278 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1526    604243123 :                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    152097278 :                                        P, Q, W, m_max, F)
    1530    152097278 :                CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1531   2775152804 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1532    198471447 :                neris = neris + mysize
    1533              :             END DO
    1534   1111138371 :             DO i = 1, mysize
    1535   1111138371 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1536              :             END DO
    1537     46374169 :             tmp_max = tmp_max*max_contraction
    1538     46374169 :             IF (tmp_max < eps_schwarz) THEN
    1539     44445201 :                RETURN
    1540              :             END IF
    1541              : 
    1542    105648970 :             DO i = 1, ncoa
    1543     69484449 :                p1 = (i - 1)*ncob
    1544    194194495 :                DO j = 1, ncob
    1545     88545525 :                   p2 = (p1 + j - 1)*ncoc
    1546    513972300 :                   DO k = 1, ncoc
    1547    355942326 :                      p3 = (p2 + k - 1)*ncod
    1548   1352962025 :                      DO l = 1, ncod
    1549    908474174 :                         p4 = p3 + l
    1550   1264416500 :                         work2(i, j, k, l) = work(p4)
    1551              :                      END DO
    1552              :                   END DO
    1553              :                END DO
    1554              :             END DO
    1555              :          CASE (2)
    1556     34798507 :             DO i = 1, nproducts
    1557     99828696 :                A = pgf_product_list(i)%ra
    1558     99828696 :                B = pgf_product_list(i)%rb
    1559     99828696 :                C = pgf_product_list(i)%rc
    1560     99828696 :                D = pgf_product_list(i)%rd
    1561     24957174 :                Rho = pgf_product_list(i)%Rho
    1562     24957174 :                RhoInv = pgf_product_list(i)%RhoInv
    1563     99828696 :                P = pgf_product_list(i)%P
    1564     99828696 :                Q = pgf_product_list(i)%Q
    1565     99828696 :                W = pgf_product_list(i)%W
    1566     99828696 :                AB = pgf_product_list(i)%AB
    1567     99828696 :                CD = pgf_product_list(i)%CD
    1568     24957174 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1569    118855802 :                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     24957174 :                                        P, Q, W, m_max, F)
    1574              : 
    1575     24957174 :                CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
    1576    776620522 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1577     34798507 :                neris = neris + mysize
    1578              :             END DO
    1579    415153957 :             DO i = 1, mysize
    1580    415153957 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1581              :             END DO
    1582      9841333 :             tmp_max = tmp_max*max_contraction
    1583      9841333 :             IF (tmp_max < eps_schwarz) THEN
    1584              :                RETURN
    1585              :             END IF
    1586              : 
    1587     31493412 :             DO j = 1, ncob
    1588     24559103 :                p1 = (j - 1)*ncoa
    1589     59008665 :                DO i = 1, ncoa
    1590     27515253 :                   p2 = (p1 + i - 1)*ncoc
    1591    170073219 :                   DO k = 1, ncoc
    1592    117998863 :                      p3 = (p2 + k - 1)*ncod
    1593    465692047 :                      DO l = 1, ncod
    1594    320177931 :                         p4 = p3 + l
    1595    438176794 :                         work2(i, j, k, l) = work(p4)
    1596              :                      END DO
    1597              :                   END DO
    1598              :                END DO
    1599              :             END DO
    1600              :          CASE (3)
    1601     88969978 :             DO i = 1, nproducts
    1602    254527520 :                A = pgf_product_list(i)%ra
    1603    254527520 :                B = pgf_product_list(i)%rb
    1604    254527520 :                C = pgf_product_list(i)%rc
    1605    254527520 :                D = pgf_product_list(i)%rd
    1606     63631880 :                Rho = pgf_product_list(i)%Rho
    1607     63631880 :                RhoInv = pgf_product_list(i)%RhoInv
    1608    254527520 :                P = pgf_product_list(i)%P
    1609    254527520 :                Q = pgf_product_list(i)%Q
    1610    254527520 :                W = pgf_product_list(i)%W
    1611    254527520 :                AB = pgf_product_list(i)%AB
    1612    254527520 :                CD = pgf_product_list(i)%CD
    1613     63631880 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1614    240811819 :                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     63631880 :                                        P, Q, W, m_max, F)
    1619              : 
    1620     63631880 :                CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
    1621    932417975 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1622     88969978 :                neris = neris + mysize
    1623              :             END DO
    1624    453213646 :             DO i = 1, mysize
    1625    453213646 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1626              :             END DO
    1627     25338098 :             tmp_max = tmp_max*max_contraction
    1628     25338098 :             IF (tmp_max < eps_schwarz) THEN
    1629              :                RETURN
    1630              :             END IF
    1631              : 
    1632     53247410 :             DO i = 1, ncoa
    1633     35446397 :                p1 = (i - 1)*ncob
    1634     96202825 :                DO j = 1, ncob
    1635     42955415 :                   p2 = (p1 + j - 1)*ncod
    1636    277067246 :                   DO l = 1, ncod
    1637    198665434 :                      p3 = (p2 + l - 1)*ncoc
    1638    580201673 :                      DO k = 1, ncoc
    1639    338580824 :                         p4 = p3 + k
    1640    537246258 :                         work2(i, j, k, l) = work(p4)
    1641              :                      END DO
    1642              :                   END DO
    1643              :                END DO
    1644              :             END DO
    1645              :          CASE (4)
    1646     21556783 :             DO i = 1, nproducts
    1647     57889332 :                A = pgf_product_list(i)%ra
    1648     57889332 :                B = pgf_product_list(i)%rb
    1649     57889332 :                C = pgf_product_list(i)%rc
    1650     57889332 :                D = pgf_product_list(i)%rd
    1651     14472333 :                Rho = pgf_product_list(i)%Rho
    1652     14472333 :                RhoInv = pgf_product_list(i)%RhoInv
    1653     57889332 :                P = pgf_product_list(i)%P
    1654     57889332 :                Q = pgf_product_list(i)%Q
    1655     57889332 :                W = pgf_product_list(i)%W
    1656     57889332 :                AB = pgf_product_list(i)%AB
    1657     57889332 :                CD = pgf_product_list(i)%CD
    1658     14472333 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1659     65908302 :                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     14472333 :                                        P, Q, W, m_max, F)
    1664              : 
    1665     14472333 :                CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1666    386572136 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1667     21556783 :                neris = neris + mysize
    1668              :             END DO
    1669    243834712 :             DO i = 1, mysize
    1670    243834712 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1671              :             END DO
    1672      7084450 :             tmp_max = tmp_max*max_contraction
    1673      7084450 :             IF (tmp_max < eps_schwarz) THEN
    1674              :                RETURN
    1675              :             END IF
    1676              : 
    1677     23402669 :             DO j = 1, ncob
    1678     18215956 :                p1 = (j - 1)*ncoa
    1679     44434183 :                DO i = 1, ncoa
    1680     21031514 :                   p2 = (p1 + i - 1)*ncod
    1681    139053799 :                   DO l = 1, ncod
    1682     99806329 :                      p3 = (p2 + l - 1)*ncoc
    1683    312901928 :                      DO k = 1, ncoc
    1684    192064085 :                         p4 = p3 + k
    1685    291870414 :                         work2(i, j, k, l) = work(p4)
    1686              :                      END DO
    1687              :                   END DO
    1688              :                END DO
    1689              :             END DO
    1690              :          CASE (5)
    1691     94388219 :             DO i = 1, nproducts
    1692    270726204 :                A = pgf_product_list(i)%ra
    1693    270726204 :                B = pgf_product_list(i)%rb
    1694    270726204 :                C = pgf_product_list(i)%rc
    1695    270726204 :                D = pgf_product_list(i)%rd
    1696     67681551 :                Rho = pgf_product_list(i)%Rho
    1697     67681551 :                RhoInv = pgf_product_list(i)%RhoInv
    1698    270726204 :                P = pgf_product_list(i)%P
    1699    270726204 :                Q = pgf_product_list(i)%Q
    1700    270726204 :                W = pgf_product_list(i)%W
    1701    270726204 :                AB = pgf_product_list(i)%AB
    1702    270726204 :                CD = pgf_product_list(i)%CD
    1703     67681551 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1704    258122159 :                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     67681551 :                                        Q, P, W, m_max, F)
    1709              : 
    1710     67681551 :                CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
    1711   1111782984 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1712     94388219 :                neris = neris + mysize
    1713              :             END DO
    1714    571531949 :             DO i = 1, mysize
    1715    571531949 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1716              :             END DO
    1717     26706668 :             tmp_max = tmp_max*max_contraction
    1718     26706668 :             IF (tmp_max < eps_schwarz) THEN
    1719              :                RETURN
    1720              :             END IF
    1721              : 
    1722     47299137 :             DO k = 1, ncoc
    1723     28123812 :                p1 = (k - 1)*ncod
    1724     79680475 :                DO l = 1, ncod
    1725     32381338 :                   p2 = (p1 + l - 1)*ncoa
    1726    208188200 :                   DO i = 1, ncoa
    1727    147683050 :                      p3 = (p2 + i - 1)*ncob
    1728    614200312 :                      DO j = 1, ncob
    1729    434135924 :                         p4 = p3 + j
    1730    581818974 :                         work2(i, j, k, l) = work(p4)
    1731              :                      END DO
    1732              :                   END DO
    1733              :                END DO
    1734              :             END DO
    1735              :          CASE (6)
    1736     40717475 :             DO i = 1, nproducts
    1737    105569324 :                A = pgf_product_list(i)%ra
    1738    105569324 :                B = pgf_product_list(i)%rb
    1739    105569324 :                C = pgf_product_list(i)%rc
    1740    105569324 :                D = pgf_product_list(i)%rd
    1741     26392331 :                Rho = pgf_product_list(i)%Rho
    1742     26392331 :                RhoInv = pgf_product_list(i)%RhoInv
    1743    105569324 :                P = pgf_product_list(i)%P
    1744    105569324 :                Q = pgf_product_list(i)%Q
    1745    105569324 :                W = pgf_product_list(i)%W
    1746    105569324 :                AB = pgf_product_list(i)%AB
    1747    105569324 :                CD = pgf_product_list(i)%CD
    1748     26392331 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1749     95212264 :                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     26392331 :                                        Q, P, W, m_max, F)
    1754              : 
    1755     26392331 :                CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
    1756    377340334 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1757     40717475 :                neris = neris + mysize
    1758              :             END DO
    1759    226166140 :             DO i = 1, mysize
    1760    226166140 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1761              :             END DO
    1762     14325144 :             tmp_max = tmp_max*max_contraction
    1763     14325144 :             IF (tmp_max < eps_schwarz) THEN
    1764              :                RETURN
    1765              :             END IF
    1766              : 
    1767     21622078 :             DO k = 1, ncoc
    1768     12605408 :                p1 = (k - 1)*ncod
    1769     36331932 :                DO l = 1, ncod
    1770     14709854 :                   p2 = (p1 + l - 1)*ncob
    1771    100978870 :                   DO j = 1, ncob
    1772     73663608 :                      p3 = (p2 + j - 1)*ncoa
    1773    242414950 :                      DO i = 1, ncoa
    1774    154041488 :                         p4 = p3 + i
    1775    227705096 :                         work2(i, j, k, l) = work(p4)
    1776              :                      END DO
    1777              :                   END DO
    1778              :                END DO
    1779              :             END DO
    1780              :          CASE (7)
    1781     17113144 :             DO i = 1, nproducts
    1782     47061592 :                A = pgf_product_list(i)%ra
    1783     47061592 :                B = pgf_product_list(i)%rb
    1784     47061592 :                C = pgf_product_list(i)%rc
    1785     47061592 :                D = pgf_product_list(i)%rd
    1786     11765398 :                Rho = pgf_product_list(i)%Rho
    1787     11765398 :                RhoInv = pgf_product_list(i)%RhoInv
    1788     47061592 :                P = pgf_product_list(i)%P
    1789     47061592 :                Q = pgf_product_list(i)%Q
    1790     47061592 :                W = pgf_product_list(i)%W
    1791     47061592 :                AB = pgf_product_list(i)%AB
    1792     47061592 :                CD = pgf_product_list(i)%CD
    1793     11765398 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1794     63047949 :                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     11765398 :                                        Q, P, W, m_max, F)
    1799              : 
    1800     11765398 :                CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
    1801    561714897 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1802     17113144 :                neris = neris + mysize
    1803              :             END DO
    1804    349775500 :             DO i = 1, mysize
    1805    349775500 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1806              :             END DO
    1807      5347746 :             tmp_max = tmp_max*max_contraction
    1808      5347746 :             IF (tmp_max < eps_schwarz) THEN
    1809              :                RETURN
    1810              :             END IF
    1811              : 
    1812     16404561 :             DO l = 1, ncod
    1813     12723211 :                p1 = (l - 1)*ncoc
    1814     30799856 :                DO k = 1, ncoc
    1815     14395295 :                   p2 = (p1 + k - 1)*ncoa
    1816    102044406 :                   DO i = 1, ncoa
    1817     74925900 :                      p3 = (p2 + i - 1)*ncob
    1818    364305231 :                      DO j = 1, ncob
    1819    274984036 :                         p4 = p3 + j
    1820    349909936 :                         work2(i, j, k, l) = work(p4)
    1821              :                      END DO
    1822              :                   END DO
    1823              :                END DO
    1824              :             END DO
    1825              :          CASE (8)
    1826      5340698 :             DO i = 1, nproducts
    1827     12611600 :                A = pgf_product_list(i)%ra
    1828     12611600 :                B = pgf_product_list(i)%rb
    1829     12611600 :                C = pgf_product_list(i)%rc
    1830     12611600 :                D = pgf_product_list(i)%rd
    1831      3152900 :                Rho = pgf_product_list(i)%Rho
    1832      3152900 :                RhoInv = pgf_product_list(i)%RhoInv
    1833     12611600 :                P = pgf_product_list(i)%P
    1834     12611600 :                Q = pgf_product_list(i)%Q
    1835     12611600 :                W = pgf_product_list(i)%W
    1836     12611600 :                AB = pgf_product_list(i)%AB
    1837     12611600 :                CD = pgf_product_list(i)%CD
    1838      3152900 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1839     17910501 :                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      3152900 :                                        Q, P, W, m_max, F)
    1844              : 
    1845      3152900 :                CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
    1846    152600702 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1847      5340698 :                neris = neris + mysize
    1848              :             END DO
    1849    112289502 :             DO i = 1, mysize
    1850    112289502 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1851              :             END DO
    1852      2187798 :             tmp_max = tmp_max*max_contraction
    1853      2187798 :             IF (tmp_max < eps_schwarz) THEN
    1854              :                RETURN
    1855              :             END IF
    1856              : 
    1857    144646113 :             DO l = 1, ncod
    1858      5807629 :                p1 = (l - 1)*ncoc
    1859     13308656 :                DO k = 1, ncoc
    1860      5867949 :                   p2 = (p1 + k - 1)*ncob
    1861     47953156 :                   DO j = 1, ncob
    1862     36277578 :                      p3 = (p2 + j - 1)*ncoa
    1863    127184403 :                      DO i = 1, ncoa
    1864     85038876 :                         p4 = p3 + i
    1865    121316454 :                         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    144564117 :          DO i = 1, nproducts
    1873    432400168 :             A = pgf_product_list(i)%ra
    1874    432400168 :             B = pgf_product_list(i)%rb
    1875    432400168 :             C = pgf_product_list(i)%rc
    1876    432400168 :             D = pgf_product_list(i)%rd
    1877    108100042 :             Rho = pgf_product_list(i)%Rho
    1878    108100042 :             RhoInv = pgf_product_list(i)%RhoInv
    1879    432400168 :             P = pgf_product_list(i)%P
    1880    432400168 :             Q = pgf_product_list(i)%Q
    1881    432400168 :             W = pgf_product_list(i)%W
    1882    108100042 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1883    216200084 :             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    108100042 :                                     P, Q, W, m_max, F)
    1888    108100042 :             work(1) = work(1) + F(1)
    1889    144564117 :             neris = neris + mysize
    1890              :          END DO
    1891     36464075 :          work2(1, 1, 1, 1) = work(1)
    1892     36464075 :          tmp_max = max_contraction*ABS(work(1))
    1893     36464075 :          IF (tmp_max < eps_schwarz) RETURN
    1894              :       END IF
    1895              : 
    1896    129224280 :       IF (tmp_max < eps_schwarz) RETURN
    1897   8008727068 :       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    129224280 :                     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   8008727068 :                     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    472250887 :    SUBROUTINE build_quartet_data(libint, A, B, C, D, &
    1940              :                                  ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1941    472250887 :                                  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    472250887 :       CALL cp_libint_set_params_eri(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, P, Q, W, m_max, F)
    1951    472250887 :    END SUBROUTINE build_quartet_data
    1952              : 
    1953              : END MODULE hfx_libint_interface
        

Generated by: LCOV version 2.0-1