LCOV - code coverage report
Current view: top level - src - hfx_pair_list_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.7 % 337 326
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 5 5

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Routines for optimizing load balance between processes in HFX calculations
      10              : !> \par History
      11              : !>      04.2008 created [Manuel Guidon]
      12              : !>      11.2019 fixed initial value for potential_id (A. Bussy)
      13              : !> \author Manuel Guidon
      14              : ! **************************************************************************************************
      15              : MODULE hfx_pair_list_methods
      16              :    USE cell_types,                      ONLY: cell_type,&
      17              :                                               pbc
      18              :    USE gamma,                           ONLY: fgamma => fgamma_0
      19              :    USE hfx_types,                       ONLY: &
      20              :         hfx_basis_type, hfx_block_range_type, hfx_cell_type, hfx_pgf_list, hfx_pgf_product_list, &
      21              :         hfx_potential_type, hfx_screen_coeff_type, pair_list_type, pair_set_list_type
      22              :    USE input_constants,                 ONLY: &
      23              :         do_potential_TShPSC, do_potential_coulomb, do_potential_gaussian, do_potential_id, &
      24              :         do_potential_long, do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, &
      25              :         do_potential_short, do_potential_truncated
      26              :    USE kinds,                           ONLY: dp
      27              :    USE libint_wrapper,                  ONLY: prim_data_f_size
      28              :    USE mathconstants,                   ONLY: pi
      29              :    USE mp2_types,                       ONLY: pair_list_type_mp2
      30              :    USE particle_types,                  ONLY: particle_type
      31              :    USE t_c_g0,                          ONLY: t_c_g0_n
      32              :    USE t_sh_p_s_c,                      ONLY: trunc_CS_poly_n20
      33              : #include "./base/base_uses.f90"
      34              : 
      35              :    IMPLICIT NONE
      36              :    PRIVATE
      37              : 
      38              :    PUBLIC :: build_pair_list, &
      39              :              build_pair_list_mp2, &
      40              :              build_pair_list_pgf, &
      41              :              build_pgf_product_list, &
      42              :              build_atomic_pair_list, &
      43              :              pgf_product_list_size
      44              : 
      45              :    ! an initial estimate for the size of the product list
      46              :    INTEGER, SAVE :: pgf_product_list_size = 128
      47              : 
      48              : !***
      49              : 
      50              : CONTAINS
      51              : 
      52              : ! **************************************************************************************************
      53              : !> \brief ...
      54              : !> \param list1 ...
      55              : !> \param list2 ...
      56              : !> \param product_list ...
      57              : !> \param nproducts ...
      58              : !> \param log10_pmax ...
      59              : !> \param log10_eps_schwarz ...
      60              : !> \param neighbor_cells ...
      61              : !> \param cell ...
      62              : !> \param potential_parameter ...
      63              : !> \param m_max ...
      64              : !> \param do_periodic ...
      65              : ! **************************************************************************************************
      66     84555358 :    SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
      67              :                                      log10_pmax, log10_eps_schwarz, neighbor_cells, &
      68              :                                      cell, potential_parameter, m_max, do_periodic)
      69              : 
      70              :       TYPE(hfx_pgf_list)                                 :: list1, list2
      71              :       TYPE(hfx_pgf_product_list), ALLOCATABLE, &
      72              :          DIMENSION(:), INTENT(INOUT)                     :: product_list
      73              :       INTEGER, INTENT(OUT)                               :: nproducts
      74              :       REAL(dp), INTENT(IN)                               :: log10_pmax, log10_eps_schwarz
      75              :       TYPE(hfx_cell_type), DIMENSION(:), POINTER         :: neighbor_cells
      76              :       TYPE(cell_type), POINTER                           :: cell
      77              :       TYPE(hfx_potential_type)                           :: potential_parameter
      78              :       INTEGER, INTENT(IN)                                :: m_max
      79              :       LOGICAL, INTENT(IN)                                :: do_periodic
      80              : 
      81              :       INTEGER                                            :: i, j, k, l, nimages1, nimages2, tmp_i4
      82              :       LOGICAL                                            :: use_gamma
      83              :       REAL(dp) :: C11(3), den, Eta, EtaInv, factor, Fm(prim_data_f_size), G(3), num, omega2, &
      84              :          omega_corr, omega_corr2, P(3), pgf_max_1, pgf_max_2, PQ(3), Q(3), R, R1, R2, ra(3), &
      85              :          rb(3), rc(3), rd(3), Rho, RhoInv, rpq2, S1234, S1234a, S1234b, shift(3), ssss, T, &
      86              :          temp(3), temp_CC(3), temp_DD(3), tmp, tmp_D(3), W(3), Zeta1, Zeta_C, Zeta_D, ZetapEtaInv
      87              :       TYPE(hfx_pgf_product_list), ALLOCATABLE, &
      88     84555358 :          DIMENSION(:)                                    :: tmp_product_list
      89              : 
      90     84555358 :       nimages1 = list1%nimages
      91     84555358 :       nimages2 = list2%nimages
      92     84555358 :       nproducts = 0
      93     84555358 :       Zeta1 = list1%zetapzetb
      94     84555358 :       Eta = list2%zetapzetb
      95     84555358 :       EtaInv = list2%ZetaInv
      96     84555358 :       Zeta_C = list2%zeta
      97     84555358 :       Zeta_D = list2%zetb
      98     84555358 :       temp_CC = 0.0_dp
      99     84555358 :       temp_DD = 0.0_dp
     100    179500930 :       DO i = 1, nimages1
     101    379782288 :          P = list1%image_list(i)%P
     102     94945572 :          R1 = list1%image_list(i)%R
     103     94945572 :          S1234a = list1%image_list(i)%S1234
     104     94945572 :          pgf_max_1 = list1%image_list(i)%pgf_max
     105    379782288 :          ra = list1%image_list(i)%ra
     106    379782288 :          rb = list1%image_list(i)%rb
     107    314960010 :          DO j = 1, nimages2
     108    135459080 :             pgf_max_2 = list2%image_list(j)%pgf_max
     109    135459080 :             IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE
     110    373191368 :             Q = list2%image_list(j)%P
     111     93297842 :             R2 = list2%image_list(j)%R
     112     93297842 :             S1234b = list2%image_list(j)%S1234
     113    373191368 :             rc = list2%image_list(j)%ra
     114    373191368 :             rd = list2%image_list(j)%rb
     115              : 
     116     93297842 :             ZetapEtaInv = Zeta1 + Eta
     117     93297842 :             ZetapEtaInv = 1.0_dp/ZetapEtaInv
     118     93297842 :             Rho = Zeta1*Eta*ZetapEtaInv
     119     93297842 :             RhoInv = 1.0_dp/Rho
     120     93297842 :             S1234 = EXP(S1234a + S1234b)
     121     93297842 :             IF (do_periodic) THEN
     122    185331416 :                temp = P - Q
     123     46332854 :                PQ = pbc(temp, cell)
     124    185331416 :                shift = -PQ + temp
     125    185331416 :                temp_CC = rc + shift
     126    185331416 :                temp_DD = rd + shift
     127              :             END IF
     128              : 
     129   2456631312 :             DO k = 1, SIZE(neighbor_cells)
     130   2268387898 :                IF (do_periodic) THEN
     131   8885691640 :                   C11 = temp_CC + neighbor_cells(k)%cell_r(:)
     132   8885691640 :                   tmp_D = temp_DD + neighbor_cells(k)%cell_r(:)
     133              :                ELSE
     134     46964988 :                   C11 = rc
     135     46964988 :                   tmp_D = rd
     136              :                END IF
     137   9073551592 :                Q = (Zeta_C*C11 + Zeta_D*tmp_D)*EtaInv
     138   2268387898 :                rpq2 = (P(1) - Q(1))**2 + (P(2) - Q(2))**2 + (P(3) - Q(3))**2
     139              :                IF (potential_parameter%potential_type == do_potential_truncated .OR. &
     140   2268387898 :                    potential_parameter%potential_type == do_potential_short .OR. &
     141              :                    potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
     142   2203076203 :                   IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) CYCLE
     143              :                END IF
     144    196521231 :                IF (potential_parameter%potential_type == do_potential_TShPSC) THEN
     145      1160608 :                   IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius*2.0_dp)**2) CYCLE
     146              :                END IF
     147    196521231 :                nproducts = nproducts + 1
     148              : 
     149              :                ! allocate size as needed,
     150              :                ! updating the global size estimate to make this a rare event in longer simulations
     151    196521231 :                IF (nproducts > SIZE(product_list)) THEN
     152         1333 : !$OMP              ATOMIC READ
     153              :                   tmp_i4 = pgf_product_list_size
     154         1333 :                   tmp_i4 = MAX(pgf_product_list_size, (3*nproducts + 1)/2)
     155         1333 : !$OMP              ATOMIC WRITE
     156              :                   pgf_product_list_size = tmp_i4
     157      2718992 :                   ALLOCATE (tmp_product_list(SIZE(product_list)))
     158      2652342 :                   tmp_product_list(:) = product_list
     159         1333 :                   DEALLOCATE (product_list)
     160      4046867 :                   ALLOCATE (product_list(tmp_i4))
     161      2652342 :                   product_list(1:SIZE(tmp_product_list)) = tmp_product_list
     162         1333 :                   DEALLOCATE (tmp_product_list)
     163              :                END IF
     164              : 
     165    196521231 :                T = Rho*rpq2
     166    315250910 :                SELECT CASE (potential_parameter%potential_type)
     167              :                CASE (do_potential_truncated)
     168    118729679 :                   R = potential_parameter%cutoff_radius*SQRT(Rho)
     169    118729679 :                   CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
     170    118729679 :                   IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
     171    118729679 :                   factor = 2.0_dp*Pi*RhoInv
     172              :                CASE (do_potential_TShPSC)
     173      1160608 :                   R = potential_parameter%cutoff_radius*SQRT(Rho)
     174     25533376 :                   product_list(nproducts)%Fm = 0.0_dp
     175      1160608 :                   CALL trunc_CS_poly_n20(product_list(nproducts)%Fm(1), R, T, m_max)
     176      1160608 :                   factor = 2.0_dp*Pi*RhoInv
     177              :                CASE (do_potential_coulomb)
     178     62150812 :                   CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
     179     62150812 :                   factor = 2.0_dp*Pi*RhoInv
     180              :                CASE (do_potential_short)
     181      6539793 :                   CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
     182      6539793 :                   omega2 = potential_parameter%omega**2
     183      6539793 :                   omega_corr2 = omega2/(omega2 + Rho)
     184      6539793 :                   omega_corr = SQRT(omega_corr2)
     185      6539793 :                   T = T*omega_corr2
     186      6539793 :                   CALL fgamma(m_max, T, Fm)
     187      6539793 :                   tmp = -omega_corr
     188     29781474 :                   DO l = 1, m_max + 1
     189     23241681 :                      product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp
     190     29781474 :                      tmp = tmp*omega_corr2
     191              :                   END DO
     192      6539793 :                   factor = 2.0_dp*Pi*RhoInv
     193              :                CASE (do_potential_long)
     194      1087985 :                   omega2 = potential_parameter%omega**2
     195      1087985 :                   omega_corr2 = omega2/(omega2 + Rho)
     196      1087985 :                   omega_corr = SQRT(omega_corr2)
     197      1087985 :                   T = T*omega_corr2
     198      1087985 :                   CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
     199      1087985 :                   tmp = omega_corr
     200      3837019 :                   DO l = 1, m_max + 1
     201      2749034 :                      product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*tmp
     202      3837019 :                      tmp = tmp*omega_corr2
     203              :                   END DO
     204      1087985 :                   factor = 2.0_dp*Pi*RhoInv
     205              :                CASE (do_potential_mix_cl)
     206       830676 :                   CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
     207       830676 :                   omega2 = potential_parameter%omega**2
     208       830676 :                   omega_corr2 = omega2/(omega2 + Rho)
     209       830676 :                   omega_corr = SQRT(omega_corr2)
     210       830676 :                   T = T*omega_corr2
     211       830676 :                   CALL fgamma(m_max, T, Fm)
     212       830676 :                   tmp = omega_corr
     213      3047709 :                   DO l = 1, m_max + 1
     214              :                      product_list(nproducts)%Fm(l) = &
     215              :                         product_list(nproducts)%Fm(l)*potential_parameter%scale_coulomb &
     216      2217033 :                         + Fm(l)*tmp*potential_parameter%scale_longrange
     217      3047709 :                      tmp = tmp*omega_corr2
     218              :                   END DO
     219       830676 :                   factor = 2.0_dp*Pi*RhoInv
     220              :                CASE (do_potential_mix_cl_trunc)
     221              : 
     222              :                   ! truncated
     223      5940064 :                   R = potential_parameter%cutoff_radius*SQRT(rho)
     224      5940064 :                   CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
     225      5940064 :                   IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
     226              : 
     227              :                   ! Coulomb
     228      5940064 :                   CALL fgamma(m_max, T, Fm)
     229              : 
     230     23020893 :                   DO l = 1, m_max + 1
     231              :                      product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)* &
     232              :                                                      (potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
     233     23020893 :                                                      Fm(l)*potential_parameter%scale_longrange
     234              :                   END DO
     235              : 
     236              :                   ! longrange
     237      5940064 :                   omega2 = potential_parameter%omega**2
     238      5940064 :                   omega_corr2 = omega2/(omega2 + Rho)
     239      5940064 :                   omega_corr = SQRT(omega_corr2)
     240      5940064 :                   T = T*omega_corr2
     241      5940064 :                   CALL fgamma(m_max, T, Fm)
     242      5940064 :                   tmp = omega_corr
     243     23020893 :                   DO l = 1, m_max + 1
     244     17080829 :                      product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
     245     23020893 :                      tmp = tmp*omega_corr2
     246              :                   END DO
     247      5940064 :                   factor = 2.0_dp*Pi*RhoInv
     248              : 
     249              :                CASE (do_potential_gaussian)
     250            0 :                   omega2 = potential_parameter%omega**2
     251            0 :                   T = -omega2*T/(Rho + omega2)
     252            0 :                   tmp = 1.0_dp
     253            0 :                   DO l = 1, m_max + 1
     254            0 :                      product_list(nproducts)%Fm(l) = EXP(T)*tmp
     255            0 :                      tmp = tmp*omega2/(Rho + omega2)
     256              :                   END DO
     257            0 :                   factor = (Pi/(Rho + omega2))**(1.5_dp)
     258              :                CASE (do_potential_mix_lg)
     259        29282 :                   omega2 = potential_parameter%omega**2
     260        29282 :                   omega_corr2 = omega2/(omega2 + Rho)
     261        29282 :                   omega_corr = SQRT(omega_corr2)
     262        29282 :                   T = T*omega_corr2
     263        29282 :                   CALL fgamma(m_max, T, Fm)
     264        29282 :                   tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange
     265       122452 :                   DO l = 1, m_max + 1
     266        93170 :                      Fm(l) = Fm(l)*tmp
     267       122452 :                      tmp = tmp*omega_corr2
     268              :                   END DO
     269              :                   T = Rho*rpq2
     270        29282 :                   T = -omega2*T/(Rho + omega2)
     271        29282 :                   tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
     272       122452 :                   DO l = 1, m_max + 1
     273        93170 :                      product_list(nproducts)%Fm(l) = EXP(T)*tmp + Fm(l)
     274       122452 :                      tmp = tmp*omega2/(Rho + omega2)
     275              :                   END DO
     276        52332 :                   factor = 1.0_dp
     277              :                CASE (do_potential_id)
     278        52332 :                   num = list1%zeta*list1%zetb
     279        52332 :                   den = list1%zeta + list1%zetb
     280       209328 :                   ssss = -num/den*SUM((ra - rb)**2)
     281              : 
     282        52332 :                   num = den*Zeta_C
     283        52332 :                   den = den + Zeta_C
     284       209328 :                   ssss = ssss - num/den*SUM((P - rc)**2)
     285              : 
     286       209328 :                   G(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + Zeta_C*rc(:))/den
     287        52332 :                   num = den*Zeta_D
     288        52332 :                   den = den + Zeta_D
     289       209328 :                   ssss = ssss - num/den*SUM((G - rd)**2)
     290              : 
     291      1151304 :                   product_list(nproducts)%Fm(:) = EXP(ssss)
     292        52332 :                   factor = 1.0_dp
     293    196573563 :                   IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
     294              :                END SELECT
     295              : 
     296    196521231 :                tmp = (Pi*ZetapEtaInv)**3
     297    196521231 :                factor = factor*S1234*SQRT(tmp)
     298              : 
     299    744588276 :                DO l = 1, m_max + 1
     300    744588276 :                   product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*factor
     301              :                END DO
     302              : 
     303    786084924 :                W = (Zeta1*P + Eta*Q)*ZetapEtaInv
     304    786084924 :                product_list(nproducts)%ra = ra
     305    786084924 :                product_list(nproducts)%rb = rb
     306    786084924 :                product_list(nproducts)%rc = C11
     307    786084924 :                product_list(nproducts)%rd = tmp_D
     308    196521231 :                product_list(nproducts)%ZetapEtaInv = ZetapEtaInv
     309    196521231 :                product_list(nproducts)%Rho = Rho
     310    196521231 :                product_list(nproducts)%RhoInv = RhoInv
     311    786084924 :                product_list(nproducts)%P = P
     312    786084924 :                product_list(nproducts)%Q = Q
     313    786084924 :                product_list(nproducts)%W = W
     314    786084924 :                product_list(nproducts)%AB = ra - rb
     315   2993410671 :                product_list(nproducts)%CD = C11 - tmp_D
     316              :             END DO
     317              :          END DO
     318              :       END DO
     319              : 
     320     84555358 :    END SUBROUTINE build_pgf_product_list
     321              : 
     322              : ! **************************************************************************************************
     323              : !> \brief ...
     324              : !> \param npgfa ...
     325              : !> \param npgfb ...
     326              : !> \param list ...
     327              : !> \param zeta ...
     328              : !> \param zetb ...
     329              : !> \param screen1 ...
     330              : !> \param screen2 ...
     331              : !> \param pgf ...
     332              : !> \param R_pgf ...
     333              : !> \param log10_pmax ...
     334              : !> \param log10_eps_schwarz ...
     335              : !> \param ra ...
     336              : !> \param rb ...
     337              : !> \param nelements ...
     338              : !> \param neighbor_cells ...
     339              : !> \param nimages ...
     340              : !> \param do_periodic ...
     341              : ! **************************************************************************************************
     342     19406500 :    SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, &
     343              :                                   log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
     344     19406500 :                                   neighbor_cells, nimages, do_periodic)
     345              :       INTEGER, INTENT(IN)                                :: npgfa, npgfb
     346              :       TYPE(hfx_pgf_list), DIMENSION(npgfa*npgfb)         :: list
     347              :       REAL(dp), DIMENSION(1:npgfa), INTENT(IN)           :: zeta
     348              :       REAL(dp), DIMENSION(1:npgfb), INTENT(IN)           :: zetb
     349              :       REAL(dp), INTENT(IN)                               :: screen1(2), screen2(2)
     350              :       TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
     351              :          POINTER                                         :: pgf, R_pgf
     352              :       REAL(dp), INTENT(IN)                               :: log10_pmax, log10_eps_schwarz, ra(3), &
     353              :                                                             rb(3)
     354              :       INTEGER, INTENT(OUT)                               :: nelements
     355              :       TYPE(hfx_cell_type), DIMENSION(:), POINTER         :: neighbor_cells
     356              :       INTEGER                                            :: nimages(npgfa*npgfb)
     357              :       LOGICAL, INTENT(IN)                                :: do_periodic
     358              : 
     359              :       INTEGER                                            :: element_counter, i, ipgf, j, jpgf
     360              :       REAL(dp)                                           :: AB(3), im_B(3), pgf_max, rab2, Zeta1, &
     361              :                                                             Zeta_A, Zeta_B, ZetaInv
     362              : 
     363     72303266 :       nimages = 0
     364              :       ! ** inner loop may never be reached
     365     19406500 :       nelements = npgfa*npgfb
     366    144634996 :       DO i = 1, SIZE(neighbor_cells)
     367    125228496 :          IF (do_periodic) THEN
     368    454906096 :             im_B = rb + neighbor_cells(i)%cell_r(:)
     369              :          ELSE
     370     11501972 :             im_B = rb
     371              :          END IF
     372    500913984 :          AB = ra - im_B
     373    125228496 :          rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
     374    125228496 :          IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
     375              :          element_counter = 0
     376     63687162 :          DO ipgf = 1, npgfa
     377    127146622 :             DO jpgf = 1, npgfb
     378     63459460 :                element_counter = element_counter + 1
     379     63459460 :                pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
     380     63459460 :                IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
     381              :                   CYCLE
     382              :                END IF
     383     53015469 :                nimages(element_counter) = nimages(element_counter) + 1
     384     53015469 :                list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
     385    212061876 :                list(element_counter)%image_list(nimages(element_counter))%ra = ra
     386    212061876 :                list(element_counter)%image_list(nimages(element_counter))%rb = im_B
     387     53015469 :                list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
     388              : 
     389     53015469 :                Zeta_A = zeta(ipgf)
     390     53015469 :                Zeta_B = zetb(jpgf)
     391     53015469 :                Zeta1 = Zeta_A + Zeta_B
     392     53015469 :                ZetaInv = 1.0_dp/Zeta1
     393              : 
     394     53015469 :                IF (nimages(element_counter) == 1) THEN
     395     46606517 :                   list(element_counter)%ipgf = ipgf
     396     46606517 :                   list(element_counter)%jpgf = jpgf
     397     46606517 :                   list(element_counter)%zetaInv = ZetaInv
     398     46606517 :                   list(element_counter)%zetapzetb = Zeta1
     399     46606517 :                   list(element_counter)%zeta = Zeta_A
     400     46606517 :                   list(element_counter)%zetb = Zeta_B
     401              :                END IF
     402              : 
     403     53015469 :                list(element_counter)%image_list(nimages(element_counter))%S1234 = (-Zeta_A*Zeta_B*ZetaInv*rab2)
     404    212061876 :                list(element_counter)%image_list(nimages(element_counter))%P = (Zeta_A*ra + Zeta_B*im_B)*ZetaInv
     405              :                list(element_counter)%image_list(nimages(element_counter))%R = &
     406     53015469 :                   MAX(0.0_dp, R_pgf(jpgf, ipgf)%x(1)*rab2 + R_pgf(jpgf, ipgf)%x(2))
     407    212061876 :                list(element_counter)%image_list(nimages(element_counter))%ra = ra
     408    212061876 :                list(element_counter)%image_list(nimages(element_counter))%rb = im_B
     409     53015469 :                list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
     410    261738088 :                list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
     411              :             END DO
     412              :          END DO
     413    144634996 :          nelements = MAX(nelements, element_counter)
     414              :       END DO
     415     72303266 :       DO j = 1, nelements
     416     72303266 :          list(j)%nimages = nimages(j)
     417              :       END DO
     418              :       ! ** Remove unused elements
     419              : 
     420              :       element_counter = 0
     421     72303266 :       DO j = 1, nelements
     422     52896766 :          IF (list(j)%nimages == 0) CYCLE
     423     46606517 :          element_counter = element_counter + 1
     424     46606517 :          list(element_counter)%nimages = list(j)%nimages
     425     46606517 :          list(element_counter)%zetapzetb = list(j)%zetapzetb
     426     46606517 :          list(element_counter)%ZetaInv = list(j)%ZetaInv
     427     46606517 :          list(element_counter)%zeta = list(j)%zeta
     428     46606517 :          list(element_counter)%zetb = list(j)%zetb
     429     46606517 :          list(element_counter)%ipgf = list(j)%ipgf
     430     46606517 :          list(element_counter)%jpgf = list(j)%jpgf
     431    119028486 :          DO i = 1, list(j)%nimages
     432    105912235 :             list(element_counter)%image_list(i) = list(j)%image_list(i)
     433              :          END DO
     434              :       END DO
     435              : 
     436     19406500 :       nelements = element_counter
     437              : 
     438     19406500 :    END SUBROUTINE build_pair_list_pgf
     439              : 
     440              : ! **************************************************************************************************
     441              : !> \brief ...
     442              : !> \param natom ...
     443              : !> \param list ...
     444              : !> \param set_list ...
     445              : !> \param i_start ...
     446              : !> \param i_end ...
     447              : !> \param j_start ...
     448              : !> \param j_end ...
     449              : !> \param kind_of ...
     450              : !> \param basis_parameter ...
     451              : !> \param particle_set ...
     452              : !> \param do_periodic ...
     453              : !> \param coeffs_set ...
     454              : !> \param coeffs_kind ...
     455              : !> \param coeffs_kind_max0 ...
     456              : !> \param log10_eps_schwarz ...
     457              : !> \param cell ...
     458              : !> \param pmax_blocks ...
     459              : !> \param atomic_pair_list ...
     460              : ! **************************************************************************************************
     461   3353863952 :    SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
     462      3178504 :                               do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
     463      3178504 :                               pmax_blocks, atomic_pair_list)
     464              : 
     465              :       INTEGER, INTENT(IN)                                :: natom
     466              :       TYPE(pair_list_type), INTENT(OUT)                  :: list
     467              :       TYPE(pair_set_list_type), DIMENSION(:), &
     468              :          INTENT(OUT)                                     :: set_list
     469              :       INTEGER, INTENT(IN)                                :: i_start, i_end, j_start, j_end, &
     470              :                                                             kind_of(*)
     471              :       TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
     472              :          POINTER                                         :: basis_parameter
     473              :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
     474              :          POINTER                                         :: particle_set
     475              :       LOGICAL, INTENT(IN)                                :: do_periodic
     476              :       TYPE(hfx_screen_coeff_type), &
     477              :          DIMENSION(:, :, :, :), INTENT(IN), POINTER      :: coeffs_set
     478              :       TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
     479              :          INTENT(IN)                                      :: coeffs_kind
     480              :       REAL(KIND=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
     481              :       TYPE(cell_type), POINTER                           :: cell
     482              :       REAL(dp), INTENT(IN)                               :: pmax_blocks
     483              :       LOGICAL, DIMENSION(natom, natom), INTENT(IN)       :: atomic_pair_list
     484              : 
     485              :       INTEGER                                            :: iatom, ikind, iset, jatom, jkind, jset, &
     486              :                                                             n_element, nset_ij, nseta, nsetb
     487              :       REAL(KIND=dp)                                      :: rab2
     488              :       REAL(KIND=dp), DIMENSION(3)                        :: B11, pbc_B, ra, rb, temp
     489              : 
     490      3178504 :       n_element = 0
     491      3178504 :       nset_ij = 0
     492              : 
     493      6357008 :       DO iatom = i_start, i_end
     494      9535512 :          DO jatom = j_start, j_end
     495      3178504 :             IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
     496              : 
     497      3017686 :             ikind = kind_of(iatom)
     498      3017686 :             nseta = basis_parameter(ikind)%nset
     499     12070744 :             ra = particle_set(iatom)%r(:)
     500              : 
     501      3017686 :             IF (jatom < iatom) CYCLE
     502      3017686 :             jkind = kind_of(jatom)
     503      3017686 :             nsetb = basis_parameter(jkind)%nset
     504     12070744 :             rb = particle_set(jatom)%r(:)
     505              : 
     506      3017686 :             IF (do_periodic) THEN
     507      5129248 :                temp = rb - ra
     508      1282312 :                pbc_B = pbc(temp, cell)
     509      5129248 :                B11 = ra + pbc_B
     510      1282312 :                rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
     511              :             ELSE
     512      1735374 :                rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
     513      1735374 :                B11 = rb ! ra - rb
     514              :             END IF
     515      3017686 :             IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
     516              :                  coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
     517              : 
     518      3007926 :             n_element = n_element + 1
     519      9023778 :             list%elements(n_element)%pair = (/iatom, jatom/)
     520      9023778 :             list%elements(n_element)%kind_pair = (/ikind, jkind/)
     521     12031704 :             list%elements(n_element)%r1 = ra
     522     12031704 :             list%elements(n_element)%r2 = B11
     523      3007926 :             list%elements(n_element)%dist2 = rab2
     524              :             ! build a list of guaranteed overlapping sets
     525      3007926 :             list%elements(n_element)%set_bounds(1) = nset_ij + 1
     526     10779724 :             DO iset = 1, nseta
     527     33559866 :                DO jset = 1, nsetb
     528     22780142 :                   IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
     529              :                       coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
     530     21409434 :                   nset_ij = nset_ij + 1
     531     73370808 :                   set_list(nset_ij)%pair = (/iset, jset/)
     532              :                END DO
     533              :             END DO
     534      6357008 :             list%elements(n_element)%set_bounds(2) = nset_ij
     535              :          END DO
     536              :       END DO
     537              : 
     538      3178504 :       list%n_element = n_element
     539              : 
     540      3178504 :    END SUBROUTINE build_pair_list
     541              : 
     542              : ! **************************************************************************************************
     543              : !> \brief ...
     544              : !> \param natom ...
     545              : !> \param atomic_pair_list ...
     546              : !> \param kind_of ...
     547              : !> \param basis_parameter ...
     548              : !> \param particle_set ...
     549              : !> \param do_periodic ...
     550              : !> \param coeffs_kind ...
     551              : !> \param coeffs_kind_max0 ...
     552              : !> \param log10_eps_schwarz ...
     553              : !> \param cell ...
     554              : !> \param blocks ...
     555              : ! **************************************************************************************************
     556         5156 :    SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
     557         5156 :                                      do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
     558              :                                      blocks)
     559              :       INTEGER, INTENT(IN)                                :: natom
     560              :       LOGICAL, DIMENSION(natom, natom)                   :: atomic_pair_list
     561              :       INTEGER, INTENT(IN)                                :: kind_of(*)
     562              :       TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
     563              :          POINTER                                         :: basis_parameter
     564              :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
     565              :          POINTER                                         :: particle_set
     566              :       LOGICAL, INTENT(IN)                                :: do_periodic
     567              :       TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
     568              :          INTENT(IN)                                      :: coeffs_kind
     569              :       REAL(KIND=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
     570              :       TYPE(cell_type), POINTER                           :: cell
     571              :       TYPE(hfx_block_range_type), DIMENSION(:), &
     572              :          INTENT(IN), POINTER                             :: blocks
     573              : 
     574              :       INTEGER                                            :: iatom, iatom_end, iatom_start, iblock, &
     575              :                                                             ikind, jatom, jatom_end, jatom_start, &
     576              :                                                             jblock, jkind, nseta, nsetb
     577              :       REAL(KIND=dp)                                      :: rab2
     578              :       REAL(KIND=dp), DIMENSION(3)                        :: B11, pbc_B, ra, rb, temp
     579              : 
     580        77424 :       atomic_pair_list = .FALSE.
     581              : 
     582        20828 :       DO iblock = 1, SIZE(blocks)
     583        15672 :          iatom_start = blocks(iblock)%istart
     584        15672 :          iatom_end = blocks(iblock)%iend
     585        77424 :          DO jblock = 1, SIZE(blocks)
     586        56596 :             jatom_start = blocks(jblock)%istart
     587        56596 :             jatom_end = blocks(jblock)%iend
     588              : 
     589       128864 :             DO iatom = iatom_start, iatom_end
     590        56596 :                ikind = kind_of(iatom)
     591        56596 :                nseta = basis_parameter(ikind)%nset
     592       226384 :                ra = particle_set(iatom)%r(:)
     593       169788 :                DO jatom = jatom_start, jatom_end
     594        56596 :                   IF (jatom < iatom) CYCLE
     595        36134 :                   jkind = kind_of(jatom)
     596        36134 :                   nsetb = basis_parameter(jkind)%nset
     597       144536 :                   rb = particle_set(jatom)%r(:)
     598              : 
     599        36134 :                   IF (do_periodic) THEN
     600        47952 :                      temp = rb - ra
     601        11988 :                      pbc_B = pbc(temp, cell)
     602        47952 :                      B11 = ra + pbc_B
     603        11988 :                      rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
     604              :                   ELSE
     605        24146 :                      rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
     606        24146 :                      B11 = rb ! ra - rb
     607              :                   END IF
     608        36134 :                   IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
     609              :                        coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) CYCLE
     610              : 
     611        35676 :                   atomic_pair_list(jatom, iatom) = .TRUE.
     612       113192 :                   atomic_pair_list(iatom, jatom) = .TRUE.
     613              :                END DO
     614              :             END DO
     615              :          END DO
     616              :       END DO
     617              : 
     618         5156 :    END SUBROUTINE build_atomic_pair_list
     619              : 
     620              : ! **************************************************************************************************
     621              : !> \brief ...
     622              : !> \param natom ...
     623              : !> \param list ...
     624              : !> \param set_list ...
     625              : !> \param i_start ...
     626              : !> \param i_end ...
     627              : !> \param j_start ...
     628              : !> \param j_end ...
     629              : !> \param kind_of ...
     630              : !> \param basis_parameter ...
     631              : !> \param particle_set ...
     632              : !> \param do_periodic ...
     633              : !> \param coeffs_set ...
     634              : !> \param coeffs_kind ...
     635              : !> \param coeffs_kind_max0 ...
     636              : !> \param log10_eps_schwarz ...
     637              : !> \param cell ...
     638              : !> \param pmax_blocks ...
     639              : !> \param atomic_pair_list ...
     640              : !> \param skip_atom_symmetry ...
     641              : ! **************************************************************************************************
     642       595568 :    SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
     643         2994 :                                   do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
     644         2994 :                                   pmax_blocks, atomic_pair_list, skip_atom_symmetry)
     645              : 
     646              :       INTEGER, INTENT(IN)                                :: natom
     647              :       TYPE(pair_list_type_mp2)                           :: list
     648              :       TYPE(pair_set_list_type), DIMENSION(:), &
     649              :          INTENT(OUT)                                     :: set_list
     650              :       INTEGER, INTENT(IN)                                :: i_start, i_end, j_start, j_end, &
     651              :                                                             kind_of(*)
     652              :       TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
     653              :          POINTER                                         :: basis_parameter
     654              :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
     655              :          POINTER                                         :: particle_set
     656              :       LOGICAL, INTENT(IN)                                :: do_periodic
     657              :       TYPE(hfx_screen_coeff_type), &
     658              :          DIMENSION(:, :, :, :), INTENT(IN), POINTER      :: coeffs_set
     659              :       TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
     660              :          INTENT(IN)                                      :: coeffs_kind
     661              :       REAL(KIND=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
     662              :       TYPE(cell_type), POINTER                           :: cell
     663              :       REAL(dp), INTENT(IN)                               :: pmax_blocks
     664              :       LOGICAL, DIMENSION(natom, natom), INTENT(IN)       :: atomic_pair_list
     665              :       LOGICAL, INTENT(IN), OPTIONAL                      :: skip_atom_symmetry
     666              : 
     667              :       INTEGER                                            :: iatom, ikind, iset, jatom, jkind, jset, &
     668              :                                                             n_element, nset_ij, nseta, nsetb
     669              :       LOGICAL                                            :: my_skip_atom_symmetry
     670              :       REAL(KIND=dp)                                      :: rab2
     671              :       REAL(KIND=dp), DIMENSION(3)                        :: B11, pbc_B, ra, rb, temp
     672              : 
     673         2994 :       n_element = 0
     674         2994 :       nset_ij = 0
     675              : 
     676         2994 :       my_skip_atom_symmetry = .FALSE.
     677         2994 :       IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry
     678              : 
     679         6076 :       DO iatom = i_start, i_end
     680         9454 :          DO jatom = j_start, j_end
     681         3378 :             IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE
     682              : 
     683         3378 :             ikind = kind_of(iatom)
     684         3378 :             nseta = basis_parameter(ikind)%nset
     685        13512 :             ra = particle_set(iatom)%r(:)
     686              : 
     687         3378 :             IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) CYCLE
     688         3230 :             jkind = kind_of(jatom)
     689         3230 :             nsetb = basis_parameter(jkind)%nset
     690        12920 :             rb = particle_set(jatom)%r(:)
     691              : 
     692         3230 :             IF (do_periodic) THEN
     693            0 :                temp = rb - ra
     694            0 :                pbc_B = pbc(temp, cell)
     695            0 :                B11 = ra + pbc_B
     696            0 :                rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
     697              :             ELSE
     698         3230 :                rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
     699         3230 :                B11 = rb ! ra - rb
     700              :             END IF
     701         3230 :             IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
     702              :                  coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
     703              : 
     704         3230 :             n_element = n_element + 1
     705         9690 :             list%elements(n_element)%pair = (/iatom, jatom/)
     706         9690 :             list%elements(n_element)%kind_pair = (/ikind, jkind/)
     707        12920 :             list%elements(n_element)%r1 = ra
     708        12920 :             list%elements(n_element)%r2 = B11
     709         3230 :             list%elements(n_element)%dist2 = rab2
     710              :             ! build a list of guaranteed overlapping sets
     711         3230 :             list%elements(n_element)%set_bounds(1) = nset_ij + 1
     712        13988 :             DO iset = 1, nseta
     713        52586 :                DO jset = 1, nsetb
     714        38598 :                   IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
     715              :                       coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
     716        38526 :                   nset_ij = nset_ij + 1
     717       126408 :                   set_list(nset_ij)%pair = (/iset, jset/)
     718              :                END DO
     719              :             END DO
     720         6460 :             list%elements(n_element)%set_bounds(2) = nset_ij
     721              :          END DO
     722              :       END DO
     723              : 
     724         2994 :       list%n_element = n_element
     725              : 
     726         2994 :    END SUBROUTINE build_pair_list_mp2
     727              : 
     728              : END MODULE hfx_pair_list_methods
        

Generated by: LCOV version 2.0-1