LCOV - code coverage report
Current view: top level - src - hfx_libint_interface.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:06f838d) Lines: 96.7 % 975 943
Test Date: 2026-06-05 07:04:50 Functions: 100.0 % 6 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Interface to the Libint-Library
      10              : !> \par History
      11              : !>      11.2006 created [Manuel Guidon]
      12              : !>      11.2019 Fixed potential_id initial values (A. Bussy)
      13              : !> \author Manuel Guidon
      14              : ! **************************************************************************************************
      15              : MODULE hfx_libint_interface
      16              : 
      17              :    USE cell_types,                      ONLY: cell_type,&
      18              :                                               real_to_scaled
      19              :    USE gamma,                           ONLY: fgamma => fgamma_0
      20              :    USE hfx_contraction_methods,         ONLY: contract
      21              :    USE hfx_types,                       ONLY: hfx_pgf_product_list,&
      22              :                                               hfx_potential_type
      23              :    USE input_constants,                 ONLY: &
      24              :         do_potential_coulomb, do_potential_gaussian, do_potential_id, do_potential_long, &
      25              :         do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, &
      26              :         do_potential_truncated
      27              :    USE kinds,                           ONLY: dp,&
      28              :                                               int_8
      29              :    USE libint_wrapper,                  ONLY: cp_libint_get_derivs,&
      30              :                                               cp_libint_get_eris,&
      31              :                                               cp_libint_set_params_eri,&
      32              :                                               cp_libint_set_params_eri_deriv,&
      33              :                                               cp_libint_set_params_eri_screen,&
      34              :                                               cp_libint_t,&
      35              :                                               get_ssss_f_val,&
      36              :                                               prim_data_f_size
      37              :    USE mathconstants,                   ONLY: pi
      38              :    USE orbital_pointers,                ONLY: nco
      39              :    USE t_c_g0,                          ONLY: t_c_g0_n
      40              : #include "./base/base_uses.f90"
      41              : 
      42              :    IMPLICIT NONE
      43              :    PRIVATE
      44              :    PUBLIC :: evaluate_eri, &
      45              :              evaluate_deriv_eri, &
      46              :              evaluate_eri_screen
      47              : 
      48              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm1 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
      49              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm2 = [4, 5, 6, 1, 2, 3, 7, 8, 9, 10, 11, 12]
      50              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm3 = [1, 2, 3, 4, 5, 6, 10, 11, 12, 7, 8, 9]
      51              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm4 = [4, 5, 6, 1, 2, 3, 10, 11, 12, 7, 8, 9]
      52              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm5 = [7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6]
      53              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm6 = [7, 8, 9, 10, 11, 12, 4, 5, 6, 1, 2, 3]
      54              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm7 = [10, 11, 12, 7, 8, 9, 1, 2, 3, 4, 5, 6]
      55              :    INTEGER, DIMENSION(12), PARAMETER :: full_perm8 = [10, 11, 12, 7, 8, 9, 4, 5, 6, 1, 2, 3]
      56              : 
      57              : !***
      58              : 
      59              : CONTAINS
      60              : 
      61              : ! **************************************************************************************************
      62              : !> \brief Fill data structure used in libint
      63              : !> \param A ...
      64              : !> \param B ...
      65              : !> \param C ...
      66              : !> \param D ...
      67              : !> \param Zeta_A ...
      68              : !> \param Zeta_B ...
      69              : !> \param Zeta_C ...
      70              : !> \param Zeta_D ...
      71              : !> \param m_max ...
      72              : !> \param potential_parameter ...
      73              : !> \param libint ...
      74              : !> \param R11 ...
      75              : !> \param R22 ...
      76              : !> \par History
      77              : !>      03.2007 created [Manuel Guidon]
      78              : !> \author Manuel Guidon
      79              : ! **************************************************************************************************
      80     89121390 :    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     89121390 :       Zeta = Zeta_A + Zeta_B
      98     89121390 :       ZetaInv = 1.0_dp/Zeta
      99     89121390 :       Eta = Zeta_C + Zeta_D
     100     89121390 :       EtaInv = 1.0_dp/Eta
     101     89121390 :       ZetapEtaInv = Zeta + Eta
     102     89121390 :       ZetapEtaInv = 1.0_dp/ZetapEtaInv
     103     89121390 :       Rho = Zeta*Eta*ZetapEtaInv
     104     89121390 :       RhoInv = 1.0_dp/Rho
     105              : 
     106    356485560 :       DO i = 1, 3
     107    267364170 :          P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv
     108    267364170 :          Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv
     109    267364170 :          AB(i) = A(i) - B(i)
     110    267364170 :          CD(i) = C(i) - D(i)
     111    267364170 :          PQ(i) = P(i) - Q(i)
     112    356485560 :          W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv
     113              :       END DO
     114              : 
     115    356485560 :       AB2 = DOT_PRODUCT(AB, AB)
     116    356485560 :       CD2 = DOT_PRODUCT(CD, CD)
     117    356485560 :       PQ2 = DOT_PRODUCT(PQ, PQ)
     118              : 
     119     89121390 :       S1234 = EXP((-Zeta_A*Zeta_B*ZetaInv*AB2) + (-Zeta_C*Zeta_D*EtaInv*CD2))
     120     89121390 :       T = Rho*PQ2
     121              : 
     122    102496416 :       SELECT CASE (potential_parameter%potential_type)
     123              :       CASE (do_potential_truncated)
     124     13375026 :          R = potential_parameter%cutoff_radius*SQRT(rho)
     125     13375026 :          R1 = R11
     126     13375026 :          R2 = R22
     127     13375026 :          IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN
     128            0 :             RETURN
     129              :          END IF
     130     13375026 :          CALL t_c_g0_n(F(1), use_gamma, R, T, m_max)
     131     13375026 :          IF (use_gamma) CALL fgamma(m_max, T, F(1))
     132     13375026 :          factor = 2.0_dp*Pi*RhoInv
     133              :       CASE (do_potential_coulomb)
     134     68614956 :          CALL fgamma(m_max, T, F(1))
     135     68614956 :          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      1007172 :          omega2 = potential_parameter%omega**2
     157      1007172 :          omega_corr2 = omega2/(omega2 + Rho)
     158      1007172 :          omega_corr = SQRT(omega_corr2)
     159      1007172 :          T = T*omega_corr2
     160      1007172 :          CALL fgamma(m_max, T, F(1))
     161      1007172 :          tmp = omega_corr
     162      3623880 :          DO i = 1, m_max + 1
     163      2616708 :             F(i) = F(i)*tmp
     164      3623880 :             tmp = tmp*omega_corr2
     165              :          END DO
     166      1007172 :          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    103698114 :          IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
     256              :       END SELECT
     257              : 
     258     89121390 :       tmp = (Pi*ZetapEtaInv)**3
     259     89121390 :       factor = factor*S1234*SQRT(tmp)
     260              : 
     261    343477164 :       DO i = 1, m_max + 1
     262    343477164 :          F(i) = F(i)*factor
     263              :       END DO
     264              : 
     265     89121390 :       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    105499289 :    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    105499289 :                                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    105499289 :                                           ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F)
     313              : 
     314    105499289 :    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     34517570 :    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     34517570 :                                  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     34517570 :                                  sphi_a, sphi_b, sphi_c, sphi_d, &
     388     34517570 :                                  work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
     389     34517570 :                                  use_virial, work_virial, work2_virial, primitives_tmp_virial, &
     390     34517570 :                                  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     34517570 :       m_max = n_a + n_b + n_c + n_d
     433     34517570 :       m_max = m_max + 1
     434     34517570 :       mysize = ncoa*ncob*ncoc*ncod
     435     69035140 :       a_mysize = mysize
     436              : 
     437   4837785266 :       work = 0.0_dp
     438   9199218758 :       work2 = 0.0_dp
     439              : 
     440     34517570 :       IF (use_virial) THEN
     441    860371492 :          work_virial = 0.0_dp
     442   1484670640 :          work2_virial = 0.0_dp
     443              :       END IF
     444              : 
     445     34517570 :       perm_case = 1
     446     34517570 :       IF (n_a < n_b) THEN
     447      3349021 :          perm_case = perm_case + 1
     448              :       END IF
     449     34517570 :       IF (n_c < n_d) THEN
     450      5672323 :          perm_case = perm_case + 2
     451              :       END IF
     452     34517570 :       IF (n_a + n_b > n_c + n_d) THEN
     453      6395773 :          perm_case = perm_case + 4
     454              :       END IF
     455              : 
     456              :       SELECT CASE (perm_case)
     457              :       CASE (1)
     458     86406599 :          DO i = 1, nproducts
     459    258180064 :             A = pgf_product_list(i)%ra
     460    258180064 :             B = pgf_product_list(i)%rb
     461    258180064 :             C = pgf_product_list(i)%rc
     462    258180064 :             D = pgf_product_list(i)%rd
     463     64545016 :             Rho = pgf_product_list(i)%Rho
     464              :             RhoInv = pgf_product_list(i)%RhoInv
     465    258180064 :             P = pgf_product_list(i)%P
     466    258180064 :             Q = pgf_product_list(i)%Q
     467    258180064 :             W = pgf_product_list(i)%W
     468    258180064 :             AB = pgf_product_list(i)%AB
     469    258180064 :             CD = pgf_product_list(i)%CD
     470     64545016 :             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     64545016 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     476     64545016 :             CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize)
     477    258180064 :             DO k = 4, 6
     478   1990998319 :                DO j = 1, mysize
     479              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     480              :                                                work_forces(j, k + 3) + &
     481   1926453303 :                                                work_forces(j, k + 6))
     482              :                END DO
     483              :             END DO
     484    839085208 :             DO k = 1, 12
     485   7770358228 :                DO j = 1, mysize
     486   7705813212 :                   work(j, k) = work(j, k) + work_forces(j, k)
     487              :                END DO
     488              :             END DO
     489     64545016 :             neris = neris + 12*mysize
     490     86406599 :             IF (use_virial) THEN
     491     13768212 :                CALL real_to_scaled(scoord(1:3), A, cell)
     492     13768212 :                CALL real_to_scaled(scoord(4:6), B, cell)
     493     13768212 :                CALL real_to_scaled(scoord(7:9), C, cell)
     494     13768212 :                CALL real_to_scaled(scoord(10:12), D, cell)
     495    178986756 :                DO k = 1, 12
     496   1430151108 :                   DO j = 1, mysize
     497   5169875952 :                      DO m = 1, 3
     498   5004657408 :                         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    284200579 :          DO n = 1, 12
     506              :             tmp_max = 0.0_dp
     507   2356928976 :             DO i = 1, mysize
     508   2356928976 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     509              :             END DO
     510    262338996 :             tmp_max = tmp_max*max_contraction
     511    262338996 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     512              : 
     513    651524407 :             DO i = 1, ncoa
     514    367323828 :                p1 = (i - 1)*ncob
     515   1036733388 :                DO j = 1, ncob
     516    407070564 :                   p2 = (p1 + j - 1)*ncoc
     517   1890177420 :                   DO k = 1, ncoc
     518   1115783028 :                      p3 = (p2 + k - 1)*ncod
     519   3617443572 :                      DO l = 1, ncod
     520   2094589980 :                         p4 = p3 + l
     521   3210373008 :                         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     21861583 :          IF (use_virial) THEN
     528      6432582 :             DO n = 1, 12
     529              :                tmp_max_virial = 0.0_dp
     530    114644292 :                DO i = 1, mysize
     531              :                   tmp_max_virial = MAX(tmp_max_virial, &
     532    114644292 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     533              :                END DO
     534      5937768 :                tmp_max_virial = tmp_max_virial*max_contraction
     535      5937768 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     536              : 
     537     17003598 :                DO i = 1, ncoa
     538     10571016 :                   p1 = (i - 1)*ncob
     539     30755184 :                   DO j = 1, ncob
     540     14246400 :                      p2 = (p1 + j - 1)*ncoc
     541     71065380 :                      DO k = 1, ncoc
     542     46247964 :                         p3 = (p2 + k - 1)*ncod
     543    169200888 :                         DO l = 1, ncod
     544    108706524 :                            p4 = p3 + l
     545    481074060 :                            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      6176104 :          DO i = 1, nproducts
     554     19354916 :             A = pgf_product_list(i)%ra
     555     19354916 :             B = pgf_product_list(i)%rb
     556     19354916 :             C = pgf_product_list(i)%rc
     557     19354916 :             D = pgf_product_list(i)%rd
     558      4838729 :             Rho = pgf_product_list(i)%Rho
     559              :             RhoInv = pgf_product_list(i)%RhoInv
     560     19354916 :             P = pgf_product_list(i)%P
     561     19354916 :             Q = pgf_product_list(i)%Q
     562     19354916 :             W = pgf_product_list(i)%W
     563     19354916 :             AB = pgf_product_list(i)%AB
     564     19354916 :             CD = pgf_product_list(i)%CD
     565      4838729 :             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      4838729 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     571      4838729 :             CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize)
     572     19354916 :             DO k = 4, 6
     573    347512943 :                DO j = 1, mysize
     574              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     575              :                                                work_forces(j, k + 3) + &
     576    342674214 :                                                work_forces(j, k + 6))
     577              :                END DO
     578              :             END DO
     579     62903477 :             DO k = 1, 12
     580   1375535585 :                DO j = 1, mysize
     581   1370696856 :                   work(j, k) = work(j, k) + work_forces(j, k)
     582              :                END DO
     583              :             END DO
     584      4838729 :             neris = neris + 12*mysize
     585      6176104 :             IF (use_virial) THEN
     586       802687 :                CALL real_to_scaled(scoord(1:3), B, cell)
     587       802687 :                CALL real_to_scaled(scoord(4:6), A, cell)
     588       802687 :                CALL real_to_scaled(scoord(7:9), C, cell)
     589       802687 :                CALL real_to_scaled(scoord(10:12), D, cell)
     590     10434931 :                DO k = 1, 12
     591    209737051 :                   DO j = 1, mysize
     592    806840724 :                      DO m = 1, 3
     593    797208480 :                         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     17385875 :          DO n = 1, 12
     601              :             tmp_max = 0.0_dp
     602    433075488 :             DO i = 1, mysize
     603    433075488 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     604              :             END DO
     605     16048500 :             tmp_max = tmp_max*max_contraction
     606     16048500 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     607              : 
     608     68618735 :             DO j = 1, ncob
     609     51232860 :                p1 = (j - 1)*ncoa
     610    120390396 :                DO i = 1, ncoa
     611     53109036 :                   p2 = (p1 + i - 1)*ncoc
     612    297486396 :                   DO k = 1, ncoc
     613    193144500 :                      p3 = (p2 + k - 1)*ncod
     614    663280524 :                      DO l = 1, ncod
     615    417026988 :                         p4 = p3 + l
     616    610171488 :                         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      1337375 :          IF (use_virial) THEN
     623      1292876 :             DO n = 1, 12
     624              :                tmp_max_virial = 0.0_dp
     625     35014056 :                DO i = 1, mysize
     626              :                   tmp_max_virial = MAX(tmp_max_virial, &
     627     35014056 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     628              :                END DO
     629      1193424 :                tmp_max_virial = tmp_max_virial*max_contraction
     630      1193424 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     631              : 
     632      5149628 :                DO j = 1, ncob
     633      3856752 :                   p1 = (j - 1)*ncoa
     634      9091248 :                   DO i = 1, ncoa
     635      4041072 :                      p2 = (p1 + i - 1)*ncoc
     636     22533984 :                      DO k = 1, ncoc
     637     14636160 :                         p3 = (p2 + k - 1)*ncod
     638     52497864 :                         DO l = 1, ncod
     639     33820632 :                            p4 = p3 + l
     640    149918688 :                            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     17793779 :          DO i = 1, nproducts
     649     54419172 :             A = pgf_product_list(i)%ra
     650     54419172 :             B = pgf_product_list(i)%rb
     651     54419172 :             C = pgf_product_list(i)%rc
     652     54419172 :             D = pgf_product_list(i)%rd
     653     13604793 :             Rho = pgf_product_list(i)%Rho
     654              :             RhoInv = pgf_product_list(i)%RhoInv
     655     54419172 :             P = pgf_product_list(i)%P
     656     54419172 :             Q = pgf_product_list(i)%Q
     657     54419172 :             W = pgf_product_list(i)%W
     658     54419172 :             AB = pgf_product_list(i)%AB
     659     54419172 :             CD = pgf_product_list(i)%CD
     660     13604793 :             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     13604793 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     666     13604793 :             CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize)
     667     54419172 :             DO k = 4, 6
     668    483826263 :                DO j = 1, mysize
     669              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     670              :                                                work_forces(j, k + 3) + &
     671    470221470 :                                                work_forces(j, k + 6))
     672              :                END DO
     673              :             END DO
     674    176862309 :             DO k = 1, 12
     675   1894490673 :                DO j = 1, mysize
     676   1880885880 :                   work(j, k) = work(j, k) + work_forces(j, k)
     677              :                END DO
     678              :             END DO
     679     13604793 :             neris = neris + 12*mysize
     680     17793779 :             IF (use_virial) THEN
     681      2211138 :                CALL real_to_scaled(scoord(1:3), A, cell)
     682      2211138 :                CALL real_to_scaled(scoord(4:6), B, cell)
     683      2211138 :                CALL real_to_scaled(scoord(7:9), D, cell)
     684      2211138 :                CALL real_to_scaled(scoord(10:12), C, cell)
     685     28744794 :                DO k = 1, 12
     686    234139914 :                   DO j = 1, mysize
     687    848114136 :                      DO m = 1, 3
     688    821580480 :                         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     54456818 :          DO n = 1, 12
     696              :             tmp_max = 0.0_dp
     697    615805080 :             DO i = 1, mysize
     698    615805080 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     699              :             END DO
     700     50267832 :             tmp_max = tmp_max*max_contraction
     701     50267832 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     702              : 
     703    142824806 :             DO i = 1, ncoa
     704     88367988 :                p1 = (i - 1)*ncob
     705    236060760 :                DO j = 1, ncob
     706     97424940 :                   p2 = (p1 + j - 1)*ncod
     707    568204656 :                   DO l = 1, ncod
     708    382411728 :                      p3 = (p2 + l - 1)*ncoc
     709   1045373916 :                      DO k = 1, ncoc
     710    565537248 :                         p4 = p3 + k
     711    947948976 :                         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      4188986 :          IF (use_virial) THEN
     718      1870583 :             DO n = 1, 12
     719              :                tmp_max_virial = 0.0_dp
     720     32305632 :                DO i = 1, mysize
     721              :                   tmp_max_virial = MAX(tmp_max_virial, &
     722     32305632 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     723              :                END DO
     724      1726692 :                tmp_max_virial = tmp_max_virial*max_contraction
     725      1726692 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     726              : 
     727      5393483 :                DO i = 1, ncoa
     728      3522900 :                   p1 = (i - 1)*ncob
     729      9544332 :                   DO j = 1, ncob
     730      4294740 :                      p2 = (p1 + j - 1)*ncod
     731     26142180 :                      DO l = 1, ncod
     732     18324540 :                         p3 = (p2 + l - 1)*ncoc
     733     53198220 :                         DO k = 1, ncoc
     734     30578940 :                            p4 = p3 + k
     735    140640300 :                            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      3127635 :          DO i = 1, nproducts
     744      9575128 :             A = pgf_product_list(i)%ra
     745      9575128 :             B = pgf_product_list(i)%rb
     746      9575128 :             C = pgf_product_list(i)%rc
     747      9575128 :             D = pgf_product_list(i)%rd
     748      2393782 :             Rho = pgf_product_list(i)%Rho
     749              :             RhoInv = pgf_product_list(i)%RhoInv
     750      9575128 :             P = pgf_product_list(i)%P
     751      9575128 :             Q = pgf_product_list(i)%Q
     752      9575128 :             W = pgf_product_list(i)%W
     753      9575128 :             AB = pgf_product_list(i)%AB
     754      9575128 :             CD = pgf_product_list(i)%CD
     755      2393782 :             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      2393782 :                                   P, Q, W, m_max, pgf_product_list(i)%Fm)
     760      2393782 :             CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize)
     761      9575128 :             DO k = 4, 6
     762    132186394 :                DO j = 1, mysize
     763              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     764              :                                                work_forces(j, k + 3) + &
     765    129792612 :                                                work_forces(j, k + 6))
     766              :                END DO
     767              :             END DO
     768     31119166 :             DO k = 1, 12
     769    521564230 :                DO j = 1, mysize
     770    519170448 :                   work(j, k) = work(j, k) + work_forces(j, k)
     771              :                END DO
     772              :             END DO
     773      2393782 :             neris = neris + 12*mysize
     774      3127635 :             IF (use_virial) THEN
     775       392668 :                CALL real_to_scaled(scoord(1:3), B, cell)
     776       392668 :                CALL real_to_scaled(scoord(4:6), A, cell)
     777       392668 :                CALL real_to_scaled(scoord(7:9), D, cell)
     778       392668 :                CALL real_to_scaled(scoord(10:12), C, cell)
     779      5104684 :                DO k = 1, 12
     780     67425004 :                   DO j = 1, mysize
     781    253993296 :                      DO m = 1, 3
     782    249281280 :                         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      9540089 :          DO n = 1, 12
     790              :             tmp_max = 0.0_dp
     791    187823904 :             DO i = 1, mysize
     792    187823904 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     793              :             END DO
     794      8806236 :             tmp_max = tmp_max*max_contraction
     795      8806236 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     796              : 
     797     37300877 :             DO j = 1, ncob
     798     27760788 :                p1 = (j - 1)*ncoa
     799     65874084 :                DO i = 1, ncoa
     800     29307060 :                   p2 = (p1 + i - 1)*ncod
     801    172147788 :                   DO l = 1, ncod
     802    115079940 :                      p3 = (p2 + l - 1)*ncoc
     803    323404668 :                      DO k = 1, ncoc
     804    179017668 :                         p4 = p3 + k
     805    294097608 :                         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       733853 :          IF (use_virial) THEN
     812       671944 :             DO n = 1, 12
     813              :                tmp_max_virial = 0.0_dp
     814     14524608 :                DO i = 1, mysize
     815              :                   tmp_max_virial = MAX(tmp_max_virial, &
     816     14524608 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     817              :                END DO
     818       620256 :                tmp_max_virial = tmp_max_virial*max_contraction
     819       620256 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     820              : 
     821      2643304 :                DO j = 1, ncob
     822      1971360 :                   p1 = (j - 1)*ncoa
     823      4710432 :                   DO i = 1, ncoa
     824      2118816 :                      p2 = (p1 + i - 1)*ncod
     825     12520224 :                      DO l = 1, ncod
     826      8430048 :                         p3 = (p2 + l - 1)*ncoc
     827     24453216 :                         DO k = 1, ncoc
     828     13904352 :                            p4 = p3 + k
     829     64047456 :                            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     18568660 :          DO i = 1, nproducts
     838     56362920 :             A = pgf_product_list(i)%ra
     839     56362920 :             B = pgf_product_list(i)%rb
     840     56362920 :             C = pgf_product_list(i)%rc
     841     56362920 :             D = pgf_product_list(i)%rd
     842     14090730 :             Rho = pgf_product_list(i)%Rho
     843              :             RhoInv = pgf_product_list(i)%RhoInv
     844     56362920 :             P = pgf_product_list(i)%P
     845     56362920 :             Q = pgf_product_list(i)%Q
     846     56362920 :             W = pgf_product_list(i)%W
     847     56362920 :             AB = pgf_product_list(i)%AB
     848     56362920 :             CD = pgf_product_list(i)%CD
     849     14090730 :             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     14090730 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     854     14090730 :             CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize)
     855     56362920 :             DO k = 4, 6
     856    542730714 :                DO j = 1, mysize
     857              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     858              :                                                work_forces(j, k + 3) + &
     859    528639984 :                                                work_forces(j, k + 6))
     860              :                END DO
     861              :             END DO
     862    183179490 :             DO k = 1, 12
     863   2128650666 :                DO j = 1, mysize
     864   2114559936 :                   work(j, k) = work(j, k) + work_forces(j, k)
     865              :                END DO
     866              :             END DO
     867     14090730 :             neris = neris + 12*mysize
     868     18568660 :             IF (use_virial) THEN
     869      2459890 :                CALL real_to_scaled(scoord(1:3), C, cell)
     870      2459890 :                CALL real_to_scaled(scoord(4:6), D, cell)
     871      2459890 :                CALL real_to_scaled(scoord(7:9), A, cell)
     872      2459890 :                CALL real_to_scaled(scoord(10:12), B, cell)
     873     31978570 :                DO k = 1, 12
     874    313018906 :                   DO j = 1, mysize
     875   1153680024 :                      DO m = 1, 3
     876   1124161344 :                         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     58213090 :          DO n = 1, 12
     884              :             tmp_max = 0.0_dp
     885    681929256 :             DO i = 1, mysize
     886    681929256 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     887              :             END DO
     888     53735160 :             tmp_max = tmp_max*max_contraction
     889     53735160 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     890              : 
     891    131661598 :             DO k = 1, ncoc
     892     73448508 :                p1 = (k - 1)*ncod
     893    205279056 :                DO l = 1, ncod
     894     78095388 :                   p2 = (p1 + l - 1)*ncoa
     895    449858472 :                   DO i = 1, ncoa
     896    298314576 :                      p3 = (p2 + i - 1)*ncob
     897   1004604060 :                      DO j = 1, ncob
     898    628194096 :                         p4 = p3 + j
     899    926508672 :                         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      4477930 :          IF (use_virial) THEN
     906      2299388 :             DO n = 1, 12
     907              :                tmp_max_virial = 0.0_dp
     908     45337308 :                DO i = 1, mysize
     909              :                   tmp_max_virial = MAX(tmp_max_virial, &
     910     45337308 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
     911              :                END DO
     912      2122512 :                tmp_max_virial = tmp_max_virial*max_contraction
     913      2122512 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
     914              : 
     915      5855948 :                DO k = 1, ncoc
     916      3556560 :                   p1 = (k - 1)*ncod
     917      9656256 :                   DO l = 1, ncod
     918      3977184 :                      p2 = (p1 + l - 1)*ncoa
     919     23284284 :                      DO i = 1, ncoa
     920     15750540 :                         p3 = (p2 + i - 1)*ncob
     921     62942520 :                         DO j = 1, ncob
     922     43214796 :                            p4 = p3 + j
     923    188609724 :                            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      4753861 :          DO i = 1, nproducts
     932     14342008 :             A = pgf_product_list(i)%ra
     933     14342008 :             B = pgf_product_list(i)%rb
     934     14342008 :             C = pgf_product_list(i)%rc
     935     14342008 :             D = pgf_product_list(i)%rd
     936      3585502 :             Rho = pgf_product_list(i)%Rho
     937              :             RhoInv = pgf_product_list(i)%RhoInv
     938     14342008 :             P = pgf_product_list(i)%P
     939     14342008 :             Q = pgf_product_list(i)%Q
     940     14342008 :             W = pgf_product_list(i)%W
     941     14342008 :             AB = pgf_product_list(i)%AB
     942     14342008 :             CD = pgf_product_list(i)%CD
     943      3585502 :             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      3585502 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
     949      3585502 :             CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize)
     950     14342008 :             DO k = 4, 6
     951    137771320 :                DO j = 1, mysize
     952              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
     953              :                                                work_forces(j, k + 3) + &
     954    134185818 :                                                work_forces(j, k + 6))
     955              :                END DO
     956              :             END DO
     957     46611526 :             DO k = 1, 12
     958    540328774 :                DO j = 1, mysize
     959    536743272 :                   work(j, k) = work(j, k) + work_forces(j, k)
     960              :                END DO
     961              :             END DO
     962      3585502 :             neris = neris + 12*mysize
     963      4753861 :             IF (use_virial) THEN
     964       442519 :                CALL real_to_scaled(scoord(1:3), C, cell)
     965       442519 :                CALL real_to_scaled(scoord(4:6), D, cell)
     966       442519 :                CALL real_to_scaled(scoord(7:9), B, cell)
     967       442519 :                CALL real_to_scaled(scoord(10:12), A, cell)
     968      5752747 :                DO k = 1, 12
     969     54042139 :                   DO j = 1, mysize
     970    198467796 :                      DO m = 1, 3
     971    193157568 :                         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     15188667 :          DO n = 1, 12
     979              :             tmp_max = 0.0_dp
     980    196283664 :             DO i = 1, mysize
     981    196283664 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
     982              :             END DO
     983     14020308 :             tmp_max = tmp_max*max_contraction
     984     14020308 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
     985              : 
     986     33207435 :             DO k = 1, ncoc
     987     18018768 :                p1 = (k - 1)*ncod
     988     52958148 :                DO l = 1, ncod
     989     20919072 :                   p2 = (p1 + l - 1)*ncob
     990    132426060 :                   DO j = 1, ncob
     991     93488220 :                      p3 = (p2 + j - 1)*ncoa
     992    296670648 :                      DO i = 1, ncoa
     993    182263356 :                         p4 = p3 + i
     994    275751576 :                         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      1168359 :          IF (use_virial) THEN
    1001       846261 :             DO n = 1, 12
    1002              :                tmp_max_virial = 0.0_dp
    1003     16358832 :                DO i = 1, mysize
    1004              :                   tmp_max_virial = MAX(tmp_max_virial, &
    1005     16358832 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1006              :                END DO
    1007       781164 :                tmp_max_virial = tmp_max_virial*max_contraction
    1008       781164 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1009              : 
    1010      1983777 :                DO k = 1, ncoc
    1011      1137516 :                   p1 = (k - 1)*ncod
    1012      3351108 :                   DO l = 1, ncod
    1013      1432428 :                      p2 = (p1 + l - 1)*ncob
    1014      9595164 :                      DO j = 1, ncob
    1015      7025220 :                         p3 = (p2 + j - 1)*ncoa
    1016     24035316 :                         DO i = 1, ncoa
    1017     15577668 :                            p4 = p3 + i
    1018     69335892 :                            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      2815690 :          DO i = 1, nproducts
    1027      8702560 :             A = pgf_product_list(i)%ra
    1028      8702560 :             B = pgf_product_list(i)%rb
    1029      8702560 :             C = pgf_product_list(i)%rc
    1030      8702560 :             D = pgf_product_list(i)%rd
    1031      2175640 :             Rho = pgf_product_list(i)%Rho
    1032              :             RhoInv = pgf_product_list(i)%RhoInv
    1033      8702560 :             P = pgf_product_list(i)%P
    1034      8702560 :             Q = pgf_product_list(i)%Q
    1035      8702560 :             W = pgf_product_list(i)%W
    1036      8702560 :             AB = pgf_product_list(i)%AB
    1037      8702560 :             CD = pgf_product_list(i)%CD
    1038      2175640 :             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      2175640 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1044      2175640 :             CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize)
    1045      8702560 :             DO k = 4, 6
    1046    212266468 :                DO j = 1, mysize
    1047              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1048              :                                                work_forces(j, k + 3) + &
    1049    210090828 :                                                work_forces(j, k + 6))
    1050              :                END DO
    1051              :             END DO
    1052     28283320 :             DO k = 1, 12
    1053    842538952 :                DO j = 1, mysize
    1054    840363312 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1055              :                END DO
    1056              :             END DO
    1057      2175640 :             neris = neris + 12*mysize
    1058      2815690 :             IF (use_virial) THEN
    1059       389649 :                CALL real_to_scaled(scoord(1:3), D, cell)
    1060       389649 :                CALL real_to_scaled(scoord(4:6), C, cell)
    1061       389649 :                CALL real_to_scaled(scoord(7:9), A, cell)
    1062       389649 :                CALL real_to_scaled(scoord(10:12), B, cell)
    1063      5065437 :                DO k = 1, 12
    1064    141227301 :                   DO j = 1, mysize
    1065    549323244 :                      DO m = 1, 3
    1066    544647456 :                         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      8320650 :          DO n = 1, 12
    1074              :             tmp_max = 0.0_dp
    1075    274647960 :             DO i = 1, mysize
    1076    274647960 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1077              :             END DO
    1078      7680600 :             tmp_max = tmp_max*max_contraction
    1079      7680600 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1080              : 
    1081     32082522 :             DO l = 1, ncod
    1082     23761872 :                p1 = (l - 1)*ncoc
    1083     55646424 :                DO k = 1, ncoc
    1084     24203952 :                   p2 = (p1 + k - 1)*ncoa
    1085    146772000 :                   DO i = 1, ncoa
    1086     98806176 :                      p3 = (p2 + i - 1)*ncob
    1087    389977488 :                      DO j = 1, ncob
    1088    266967360 :                         p4 = p3 + j
    1089    365773536 :                         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       640050 :          IF (use_virial) THEN
    1096       650013 :             DO n = 1, 12
    1097              :                tmp_max_virial = 0.0_dp
    1098     22174416 :                DO i = 1, mysize
    1099              :                   tmp_max_virial = MAX(tmp_max_virial, &
    1100     22174416 :                                        ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3)))
    1101              :                END DO
    1102       600012 :                tmp_max_virial = tmp_max_virial*max_contraction
    1103       600012 :                tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial)
    1104              : 
    1105      2506641 :                DO l = 1, ncod
    1106      1856628 :                   p1 = (l - 1)*ncoc
    1107      4350132 :                   DO k = 1, ncoc
    1108      1893492 :                      p2 = (p1 + k - 1)*ncoa
    1109     10999188 :                      DO i = 1, ncoa
    1110      7249068 :                         p3 = (p2 + i - 1)*ncob
    1111     30716964 :                         DO j = 1, ncob
    1112     21574404 :                            p4 = p3 + j
    1113     93546684 :                            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       374531 :          DO i = 1, nproducts
    1122      1060388 :             A = pgf_product_list(i)%ra
    1123      1060388 :             B = pgf_product_list(i)%rb
    1124      1060388 :             C = pgf_product_list(i)%rc
    1125      1060388 :             D = pgf_product_list(i)%rd
    1126       265097 :             Rho = pgf_product_list(i)%Rho
    1127              :             RhoInv = pgf_product_list(i)%RhoInv
    1128      1060388 :             P = pgf_product_list(i)%P
    1129      1060388 :             Q = pgf_product_list(i)%Q
    1130      1060388 :             W = pgf_product_list(i)%W
    1131      1060388 :             AB = pgf_product_list(i)%AB
    1132      1060388 :             CD = pgf_product_list(i)%CD
    1133       265097 :             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       265097 :                                   Q, P, W, m_max, pgf_product_list(i)%Fm)
    1139       265097 :             CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize)
    1140      1060388 :             DO k = 4, 6
    1141     34002656 :                DO j = 1, mysize
    1142              :                   work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
    1143              :                                                work_forces(j, k + 3) + &
    1144     33737559 :                                                work_forces(j, k + 6))
    1145              :                END DO
    1146              :             END DO
    1147      3446261 :             DO k = 1, 12
    1148    135215333 :                DO j = 1, mysize
    1149    134950236 :                   work(j, k) = work(j, k) + work_forces(j, k)
    1150              :                END DO
    1151              :             END DO
    1152       265097 :             neris = neris + 12*mysize
    1153       374531 :             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      1422642 :          DO n = 1, 12
    1169              :             tmp_max = 0.0_dp
    1170     56773368 :             DO i = 1, mysize
    1171     56773368 :                tmp_max = MAX(tmp_max, ABS(work(i, n)))
    1172              :             END DO
    1173      1313208 :             tmp_max = tmp_max*max_contraction
    1174      1313208 :             tmp_max_all = MAX(tmp_max_all, tmp_max)
    1175              : 
    1176      5788146 :             DO l = 1, ncod
    1177      4365504 :                p1 = (l - 1)*ncoc
    1178     10044216 :                DO k = 1, ncoc
    1179      4365504 :                   p2 = (p1 + k - 1)*ncob
    1180     34924032 :                   DO j = 1, ncob
    1181     26193024 :                      p3 = (p2 + j - 1)*ncoa
    1182     86018688 :                      DO i = 1, ncoa
    1183     55460160 :                         p4 = p3 + i
    1184     81653184 :                         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     34627004 :          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     34517570 :       IF (.NOT. use_virial) THEN
    1218     33426535 :          IF (tmp_max_all < eps_schwarz) RETURN
    1219              :       END IF
    1220              : 
    1221     28795269 :       IF (tmp_max_all >= eps_schwarz) THEN
    1222    372127457 :          DO i = 1, 12
    1223  20291234316 :             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    343502268 :                           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  20319859505 :                           s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
    1241              :          END DO
    1242              :       END IF
    1243              : 
    1244     28795269 :       IF (use_virial .AND. tmp_max_all_virial >= eps_schwarz) THEN
    1245     12130053 :          DO i = 1, 12
    1246     45720969 :             DO m = 1, 3
    1247   4955236956 :                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     33590916 :                              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   4966433928 :                                     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     89121390 :    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     89121390 :       m_max = n_a + n_b + n_c + n_d
    1313     89121390 :       mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
    1314    178242780 :       a_mysize = mysize
    1315              : 
    1316     89121390 :       IF (m_max /= 0) THEN
    1317     55805328 :          perm_case = 1
    1318     55805328 :          IF (n_a < n_b) THEN
    1319     22652280 :             perm_case = perm_case + 1
    1320              :          END IF
    1321     55805328 :          IF (n_c < n_d) THEN
    1322     22652280 :             perm_case = perm_case + 2
    1323              :          END IF
    1324     55805328 :          IF (n_a + n_b > n_c + n_d) THEN
    1325            0 :             perm_case = perm_case + 4
    1326              :          END IF
    1327              : 
    1328     33153048 :          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     33153048 :                                            potential_parameter, libint, R1, R2)
    1332              : 
    1333     33153048 :             CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1334   2817475800 :             DO i = 1, mysize
    1335   2817475800 :                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     22652280 :                                            potential_parameter, libint, R1, R2)
    1356              : 
    1357     22652280 :             CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1358   1010667408 :             DO i = 1, mysize
    1359   1010667408 :                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     55805328 :             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     33316062 :                                         potential_parameter, libint, R1, R2)
    1397     33316062 :          max_val = ABS(get_ssss_f_val(libint))
    1398              :       END IF
    1399              : 
    1400     89121390 :    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    181862510 :    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    181862510 :                            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    181862510 :                            sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
    1463    181862510 :                            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    181862510 :       m_max = n_a + n_b + n_c + n_d
    1496    181862510 :       mysize = ncoa*ncob*ncoc*ncod
    1497    363725020 :       a_mysize = mysize
    1498              : 
    1499   3640419540 :       work = 0.0_dp
    1500    181862510 :       IF (m_max /= 0) THEN
    1501    142507621 :          perm_case = 1
    1502    142507621 :          IF (n_a < n_b) THEN
    1503     33879249 :             perm_case = perm_case + 1
    1504              :          END IF
    1505    142507621 :          IF (n_c < n_d) THEN
    1506     40971291 :             perm_case = perm_case + 2
    1507              :          END IF
    1508    142507621 :          IF (n_a + n_b > n_c + n_d) THEN
    1509     49782429 :             perm_case = perm_case + 4
    1510              :          END IF
    1511              :          SELECT CASE (perm_case)
    1512              :          CASE (1)
    1513    214064621 :             DO i = 1, nproducts
    1514    658767120 :                A = pgf_product_list(i)%ra
    1515    658767120 :                B = pgf_product_list(i)%rb
    1516    658767120 :                C = pgf_product_list(i)%rc
    1517    658767120 :                D = pgf_product_list(i)%rd
    1518    164691780 :                Rho = pgf_product_list(i)%Rho
    1519    164691780 :                RhoInv = pgf_product_list(i)%RhoInv
    1520    658767120 :                P = pgf_product_list(i)%P
    1521    658767120 :                Q = pgf_product_list(i)%Q
    1522    658767120 :                W = pgf_product_list(i)%W
    1523    658767120 :                AB = pgf_product_list(i)%AB
    1524    658767120 :                CD = pgf_product_list(i)%CD
    1525    164691780 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1526    648999318 :                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    164691780 :                                        P, Q, W, m_max, F)
    1530    164691780 :                CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
    1531   2899546281 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1532    214064621 :                neris = neris + mysize
    1533              :             END DO
    1534   1150404592 :             DO i = 1, mysize
    1535   1150404592 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1536              :             END DO
    1537     49372841 :             tmp_max = tmp_max*max_contraction
    1538     49372841 :             IF (tmp_max < eps_schwarz) THEN
    1539     45420780 :                RETURN
    1540              :             END IF
    1541              : 
    1542    112745154 :             DO i = 1, ncoa
    1543     73893209 :                p1 = (i - 1)*ncob
    1544    206276501 :                DO j = 1, ncob
    1545     93531347 :                   p2 = (p1 + j - 1)*ncoc
    1546    540304759 :                   DO k = 1, ncoc
    1547    372880203 :                      p3 = (p2 + k - 1)*ncod
    1548   1407424567 :                      DO l = 1, ncod
    1549    941013017 :                         p4 = p3 + l
    1550   1313893220 :                         work2(i, j, k, l) = work(p4)
    1551              :                      END DO
    1552              :                   END DO
    1553              :                END DO
    1554              :             END DO
    1555              :          CASE (2)
    1556     35517685 :             DO i = 1, nproducts
    1557    101910416 :                A = pgf_product_list(i)%ra
    1558    101910416 :                B = pgf_product_list(i)%rb
    1559    101910416 :                C = pgf_product_list(i)%rc
    1560    101910416 :                D = pgf_product_list(i)%rd
    1561     25477604 :                Rho = pgf_product_list(i)%Rho
    1562     25477604 :                RhoInv = pgf_product_list(i)%RhoInv
    1563    101910416 :                P = pgf_product_list(i)%P
    1564    101910416 :                Q = pgf_product_list(i)%Q
    1565    101910416 :                W = pgf_product_list(i)%W
    1566    101910416 :                AB = pgf_product_list(i)%AB
    1567    101910416 :                CD = pgf_product_list(i)%CD
    1568     25477604 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1569    121276770 :                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     25477604 :                                        P, Q, W, m_max, F)
    1574              : 
    1575     25477604 :                CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
    1576    789548928 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1577     35517685 :                neris = neris + mysize
    1578              :             END DO
    1579    421813328 :             DO i = 1, mysize
    1580    421813328 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1581              :             END DO
    1582     10040081 :             tmp_max = tmp_max*max_contraction
    1583     10040081 :             IF (tmp_max < eps_schwarz) THEN
    1584              :                RETURN
    1585              :             END IF
    1586              : 
    1587     32224549 :             DO j = 1, ncob
    1588     25118675 :                p1 = (j - 1)*ncoa
    1589     60330082 :                DO i = 1, ncoa
    1590     28105533 :                   p2 = (p1 + i - 1)*ncoc
    1591    173512131 :                   DO k = 1, ncoc
    1592    120287923 :                      p3 = (p2 + k - 1)*ncod
    1593    473674081 :                      DO l = 1, ncod
    1594    325280625 :                         p4 = p3 + l
    1595    445568548 :                         work2(i, j, k, l) = work(p4)
    1596              :                      END DO
    1597              :                   END DO
    1598              :                END DO
    1599              :             END DO
    1600              :          CASE (3)
    1601     92620978 :             DO i = 1, nproducts
    1602    265969384 :                A = pgf_product_list(i)%ra
    1603    265969384 :                B = pgf_product_list(i)%rb
    1604    265969384 :                C = pgf_product_list(i)%rc
    1605    265969384 :                D = pgf_product_list(i)%rd
    1606     66492346 :                Rho = pgf_product_list(i)%Rho
    1607     66492346 :                RhoInv = pgf_product_list(i)%RhoInv
    1608    265969384 :                P = pgf_product_list(i)%P
    1609    265969384 :                Q = pgf_product_list(i)%Q
    1610    265969384 :                W = pgf_product_list(i)%W
    1611    265969384 :                AB = pgf_product_list(i)%AB
    1612    265969384 :                CD = pgf_product_list(i)%CD
    1613     66492346 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1614    250437482 :                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     66492346 :                                        P, Q, W, m_max, F)
    1619              : 
    1620     66492346 :                CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
    1621    953888587 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1622     92620978 :                neris = neris + mysize
    1623              :             END DO
    1624    463108112 :             DO i = 1, mysize
    1625    463108112 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1626              :             END DO
    1627     26128632 :             tmp_max = tmp_max*max_contraction
    1628     26128632 :             IF (tmp_max < eps_schwarz) THEN
    1629              :                RETURN
    1630              :             END IF
    1631              : 
    1632     55251439 :             DO i = 1, ncoa
    1633     36733726 :                p1 = (i - 1)*ncob
    1634     99636593 :                DO j = 1, ncob
    1635     44385154 :                   p2 = (p1 + j - 1)*ncod
    1636    285301434 :                   DO l = 1, ncod
    1637    204182554 :                      p3 = (p2 + l - 1)*ncoc
    1638    595180492 :                      DO k = 1, ncoc
    1639    346612784 :                         p4 = p3 + k
    1640    550795338 :                         work2(i, j, k, l) = work(p4)
    1641              :                      END DO
    1642              :                   END DO
    1643              :                END DO
    1644              :             END DO
    1645              :          CASE (4)
    1646     21908309 :             DO i = 1, nproducts
    1647     58898684 :                A = pgf_product_list(i)%ra
    1648     58898684 :                B = pgf_product_list(i)%rb
    1649     58898684 :                C = pgf_product_list(i)%rc
    1650     58898684 :                D = pgf_product_list(i)%rd
    1651     14724671 :                Rho = pgf_product_list(i)%Rho
    1652     14724671 :                RhoInv = pgf_product_list(i)%RhoInv
    1653     58898684 :                P = pgf_product_list(i)%P
    1654     58898684 :                Q = pgf_product_list(i)%Q
    1655     58898684 :                W = pgf_product_list(i)%W
    1656     58898684 :                AB = pgf_product_list(i)%AB
    1657     58898684 :                CD = pgf_product_list(i)%CD
    1658     14724671 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1659     66975971 :                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     14724671 :                                        P, Q, W, m_max, F)
    1664              : 
    1665     14724671 :                CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
    1666    390905968 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1667     21908309 :                neris = neris + mysize
    1668              :             END DO
    1669    246609747 :             DO i = 1, mysize
    1670    246609747 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1671              :             END DO
    1672      7183638 :             tmp_max = tmp_max*max_contraction
    1673      7183638 :             IF (tmp_max < eps_schwarz) THEN
    1674              :                RETURN
    1675              :             END IF
    1676              : 
    1677     23815341 :             DO j = 1, ncob
    1678     18531172 :                p1 = (j - 1)*ncoa
    1679     45194519 :                DO i = 1, ncoa
    1680     21379178 :                   p2 = (p1 + i - 1)*ncod
    1681    141180295 :                   DO l = 1, ncod
    1682    101269945 :                      p3 = (p2 + l - 1)*ncoc
    1683    317271908 :                      DO k = 1, ncoc
    1684    194622785 :                         p4 = p3 + k
    1685    295892730 :                         work2(i, j, k, l) = work(p4)
    1686              :                      END DO
    1687              :                   END DO
    1688              :                END DO
    1689              :             END DO
    1690              :          CASE (5)
    1691     98837783 :             DO i = 1, nproducts
    1692    284668412 :                A = pgf_product_list(i)%ra
    1693    284668412 :                B = pgf_product_list(i)%rb
    1694    284668412 :                C = pgf_product_list(i)%rc
    1695    284668412 :                D = pgf_product_list(i)%rd
    1696     71167103 :                Rho = pgf_product_list(i)%Rho
    1697     71167103 :                RhoInv = pgf_product_list(i)%RhoInv
    1698    284668412 :                P = pgf_product_list(i)%P
    1699    284668412 :                Q = pgf_product_list(i)%Q
    1700    284668412 :                W = pgf_product_list(i)%W
    1701    284668412 :                AB = pgf_product_list(i)%AB
    1702    284668412 :                CD = pgf_product_list(i)%CD
    1703     71167103 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1704    269550346 :                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     71167103 :                                        Q, P, W, m_max, F)
    1709              : 
    1710     71167103 :                CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
    1711   1138550087 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1712     98837783 :                neris = neris + mysize
    1713              :             END DO
    1714    583394454 :             DO i = 1, mysize
    1715    583394454 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1716              :             END DO
    1717     27670680 :             tmp_max = tmp_max*max_contraction
    1718     27670680 :             IF (tmp_max < eps_schwarz) THEN
    1719              :                RETURN
    1720              :             END IF
    1721              : 
    1722     49191261 :             DO k = 1, ncoc
    1723     29172163 :                p1 = (k - 1)*ncod
    1724     82683404 :                DO l = 1, ncod
    1725     33492143 :                   p2 = (p1 + l - 1)*ncoa
    1726    214511650 :                   DO i = 1, ncoa
    1727    151847344 :                      p3 = (p2 + i - 1)*ncob
    1728    627864971 :                      DO j = 1, ncob
    1729    442525484 :                         p4 = p3 + j
    1730    594372828 :                         work2(i, j, k, l) = work(p4)
    1731              :                      END DO
    1732              :                   END DO
    1733              :                END DO
    1734              :             END DO
    1735              :          CASE (6)
    1736     41106485 :             DO i = 1, nproducts
    1737    106615028 :                A = pgf_product_list(i)%ra
    1738    106615028 :                B = pgf_product_list(i)%rb
    1739    106615028 :                C = pgf_product_list(i)%rc
    1740    106615028 :                D = pgf_product_list(i)%rd
    1741     26653757 :                Rho = pgf_product_list(i)%Rho
    1742     26653757 :                RhoInv = pgf_product_list(i)%RhoInv
    1743    106615028 :                P = pgf_product_list(i)%P
    1744    106615028 :                Q = pgf_product_list(i)%Q
    1745    106615028 :                W = pgf_product_list(i)%W
    1746    106615028 :                AB = pgf_product_list(i)%AB
    1747    106615028 :                CD = pgf_product_list(i)%CD
    1748     26653757 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1749     96081767 :                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     26653757 :                                        Q, P, W, m_max, F)
    1754              : 
    1755     26653757 :                CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
    1756    380327605 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1757     41106485 :                neris = neris + mysize
    1758              :             END DO
    1759    228689014 :             DO i = 1, mysize
    1760    228689014 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1761              :             END DO
    1762     14452728 :             tmp_max = tmp_max*max_contraction
    1763     14452728 :             IF (tmp_max < eps_schwarz) THEN
    1764              :                RETURN
    1765              :             END IF
    1766              : 
    1767     21875791 :             DO k = 1, ncoc
    1768     12749905 :                p1 = (k - 1)*ncod
    1769     36755966 :                DO l = 1, ncod
    1770     14880175 :                   p2 = (p1 + l - 1)*ncob
    1771    102059795 :                   DO j = 1, ncob
    1772     74429715 :                      p3 = (p2 + j - 1)*ncoa
    1773    244974513 :                      DO i = 1, ncoa
    1774    155664623 :                         p4 = p3 + i
    1775    230094338 :                         work2(i, j, k, l) = work(p4)
    1776              :                      END DO
    1777              :                   END DO
    1778              :                END DO
    1779              :             END DO
    1780              :          CASE (7)
    1781     17495552 :             DO i = 1, nproducts
    1782     48157332 :                A = pgf_product_list(i)%ra
    1783     48157332 :                B = pgf_product_list(i)%rb
    1784     48157332 :                C = pgf_product_list(i)%rc
    1785     48157332 :                D = pgf_product_list(i)%rd
    1786     12039333 :                Rho = pgf_product_list(i)%Rho
    1787     12039333 :                RhoInv = pgf_product_list(i)%RhoInv
    1788     48157332 :                P = pgf_product_list(i)%P
    1789     48157332 :                Q = pgf_product_list(i)%Q
    1790     48157332 :                W = pgf_product_list(i)%W
    1791     48157332 :                AB = pgf_product_list(i)%AB
    1792     48157332 :                CD = pgf_product_list(i)%CD
    1793     12039333 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1794     64453034 :                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     12039333 :                                        Q, P, W, m_max, F)
    1799              : 
    1800     12039333 :                CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
    1801    571177130 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1802     17495552 :                neris = neris + mysize
    1803              :             END DO
    1804    354606183 :             DO i = 1, mysize
    1805    354606183 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1806              :             END DO
    1807      5456219 :             tmp_max = tmp_max*max_contraction
    1808      5456219 :             IF (tmp_max < eps_schwarz) THEN
    1809              :                RETURN
    1810              :             END IF
    1811              : 
    1812     16712919 :             DO l = 1, ncod
    1813     12957919 :                p1 = (l - 1)*ncoc
    1814     31349594 :                DO k = 1, ncoc
    1815     14636675 :                   p2 = (p1 + k - 1)*ncoa
    1816    103614759 :                   DO i = 1, ncoa
    1817     76020165 :                      p3 = (p2 + i - 1)*ncob
    1818    368742393 :                      DO j = 1, ncob
    1819    278085553 :                         p4 = p3 + j
    1820    354105718 :                         work2(i, j, k, l) = work(p4)
    1821              :                      END DO
    1822              :                   END DO
    1823              :                END DO
    1824              :             END DO
    1825              :          CASE (8)
    1826      5370230 :             DO i = 1, nproducts
    1827     12669712 :                A = pgf_product_list(i)%ra
    1828     12669712 :                B = pgf_product_list(i)%rb
    1829     12669712 :                C = pgf_product_list(i)%rc
    1830     12669712 :                D = pgf_product_list(i)%rd
    1831      3167428 :                Rho = pgf_product_list(i)%Rho
    1832      3167428 :                RhoInv = pgf_product_list(i)%RhoInv
    1833     12669712 :                P = pgf_product_list(i)%P
    1834     12669712 :                Q = pgf_product_list(i)%Q
    1835     12669712 :                W = pgf_product_list(i)%W
    1836     12669712 :                AB = pgf_product_list(i)%AB
    1837     12669712 :                CD = pgf_product_list(i)%CD
    1838      3167428 :                ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1839     17995709 :                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      3167428 :                                        Q, P, W, m_max, F)
    1844              : 
    1845      3167428 :                CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
    1846    153381898 :                work(1:mysize) = work(1:mysize) + p_work(1:mysize)
    1847      5370230 :                neris = neris + mysize
    1848              :             END DO
    1849    113084332 :             DO i = 1, mysize
    1850    113084332 :                tmp_max = MAX(tmp_max, ABS(work(i)))
    1851              :             END DO
    1852      2202802 :             tmp_max = tmp_max*max_contraction
    1853      2202802 :             IF (tmp_max < eps_schwarz) THEN
    1854              :                RETURN
    1855              :             END IF
    1856              : 
    1857    149996322 :             DO l = 1, ncod
    1858      5844730 :                p1 = (l - 1)*ncoc
    1859     13393751 :                DO k = 1, ncoc
    1860      5905050 :                   p2 = (p1 + k - 1)*ncob
    1861     48249964 :                   DO j = 1, ncob
    1862     36500184 :                      p3 = (p2 + j - 1)*ncoa
    1863    128018364 :                      DO i = 1, ncoa
    1864     85613130 :                         p4 = p3 + i
    1865    122113314 :                         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    162422585 :          DO i = 1, nproducts
    1873    492270784 :             A = pgf_product_list(i)%ra
    1874    492270784 :             B = pgf_product_list(i)%rb
    1875    492270784 :             C = pgf_product_list(i)%rc
    1876    492270784 :             D = pgf_product_list(i)%rd
    1877    123067696 :             Rho = pgf_product_list(i)%Rho
    1878    123067696 :             RhoInv = pgf_product_list(i)%RhoInv
    1879    492270784 :             P = pgf_product_list(i)%P
    1880    492270784 :             Q = pgf_product_list(i)%Q
    1881    492270784 :             W = pgf_product_list(i)%W
    1882    123067696 :             ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv
    1883    246135392 :             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    123067696 :                                     P, Q, W, m_max, F)
    1888    123067696 :             work(1) = work(1) + F(1)
    1889    162422585 :             neris = neris + mysize
    1890              :          END DO
    1891     39354889 :          work2(1, 1, 1, 1) = work(1)
    1892     39354889 :          tmp_max = max_contraction*ABS(work(1))
    1893     39354889 :          IF (tmp_max < eps_schwarz) RETURN
    1894              :       END IF
    1895              : 
    1896    136441730 :       IF (tmp_max < eps_schwarz) RETURN
    1897   8376307702 :       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    136441730 :                     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   8376307702 :                     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    507481718 :    SUBROUTINE build_quartet_data(libint, A, B, C, D, &
    1940              :                                  ZetaInv, EtaInv, ZetapEtaInv, Rho, &
    1941    507481718 :                                  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    507481718 :       CALL cp_libint_set_params_eri(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, P, Q, W, m_max, F)
    1951    507481718 :    END SUBROUTINE build_quartet_data
    1952              : 
    1953              : END MODULE hfx_libint_interface
        

Generated by: LCOV version 2.0-1