!--------------------------------------------------------------------------------------------------! ! CP2K: A general program to perform molecular dynamics simulations ! ! Copyright (C) 2000 - 2019 CP2K developers group ! !--------------------------------------------------------------------------------------------------! ! ************************************************************************************************** !> \brief Rountines for optimizing load balance between processes in HFX calculations !> \par History !> 04.2008 created [Manuel Guidon] !> 11.2019 fixed initial value for potential_id (A. Bussy) !> \author Manuel Guidon ! ************************************************************************************************** MODULE hfx_pair_list_methods USE cell_types, ONLY: cell_type,& pbc USE gamma, ONLY: fgamma => fgamma_0 USE hfx_types, ONLY: & hfx_basis_type, hfx_block_range_type, hfx_cell_type, hfx_pgf_list, hfx_pgf_product_list, & hfx_potential_type, hfx_screen_coeff_type, pair_list_type, pair_set_list_type USE input_constants, ONLY: & do_potential_TShPSC, do_potential_coulomb, do_potential_gaussian, do_potential_id, & do_potential_long, do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, & do_potential_short, do_potential_truncated USE kinds, ONLY: dp USE libint_wrapper, ONLY: prim_data_f_size USE mathconstants, ONLY: pi USE mp2_types, ONLY: pair_list_type_mp2 USE particle_types, ONLY: particle_type USE t_c_g0, ONLY: t_c_g0_n USE t_sh_p_s_c, ONLY: trunc_CS_poly_n20 #include "./base/base_uses.f90" IMPLICIT NONE PRIVATE PUBLIC :: build_pair_list, & build_pair_list_mp2, & build_pair_list_pgf, & build_pgf_product_list, & build_atomic_pair_list, & pgf_product_list_size ! an initial estimate for the size of the product list INTEGER, SAVE :: pgf_product_list_size = 128 !*** CONTAINS ! ************************************************************************************************** !> \brief ... !> \param list1 ... !> \param list2 ... !> \param product_list ... !> \param nproducts ... !> \param log10_pmax ... !> \param log10_eps_schwarz ... !> \param neighbor_cells ... !> \param cell ... !> \param potential_parameter ... !> \param m_max ... !> \param do_periodic ... ! ************************************************************************************************** SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, & log10_pmax, log10_eps_schwarz, neighbor_cells, & cell, potential_parameter, m_max, do_periodic) TYPE(hfx_pgf_list) :: list1, list2 TYPE(hfx_pgf_product_list), ALLOCATABLE, & DIMENSION(:), INTENT(INOUT) :: product_list INTEGER, INTENT(OUT) :: nproducts REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells TYPE(cell_type), POINTER :: cell TYPE(hfx_potential_type) :: potential_parameter INTEGER, INTENT(IN) :: m_max LOGICAL, INTENT(IN) :: do_periodic INTEGER :: i, j, k, l, nimages1, nimages2, tmp_i4 LOGICAL :: use_gamma REAL(dp) :: C11(3), den, Eta, EtaInv, factor, Fm(prim_data_f_size), G(3), num, omega2, & omega_corr, omega_corr2, P(3), pgf_max_1, pgf_max_2, PQ(3), Q(3), R, R1, R2, ra(3), & rb(3), rc(3), rd(3), Rho, RhoInv, rpq2, S1234, S1234a, S1234b, shift(3), ssss, T, & temp(3), temp_CC(3), temp_DD(3), tmp, tmp_D(3), W(3), Zeta1, Zeta_C, Zeta_D, ZetapEtaInv TYPE(hfx_pgf_product_list), ALLOCATABLE, & DIMENSION(:) :: tmp_product_list nimages1 = list1%nimages nimages2 = list2%nimages nproducts = 0 Zeta1 = list1%zetapzetb Eta = list2%zetapzetb EtaInv = list2%ZetaInv Zeta_C = list2%zeta Zeta_D = list2%zetb DO i = 1, nimages1 P = list1%image_list(i)%P R1 = list1%image_list(i)%R S1234a = list1%image_list(i)%S1234 pgf_max_1 = list1%image_list(i)%pgf_max ra = list1%image_list(i)%ra rb = list1%image_list(i)%rb DO j = 1, nimages2 pgf_max_2 = list2%image_list(j)%pgf_max IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE Q = list2%image_list(j)%P R2 = list2%image_list(j)%R S1234b = list2%image_list(j)%S1234 rc = list2%image_list(j)%ra rd = list2%image_list(j)%rb ZetapEtaInv = Zeta1 + Eta ZetapEtaInv = 1.0_dp/ZetapEtaInv Rho = Zeta1*Eta*ZetapEtaInv RhoInv = 1.0_dp/Rho S1234 = EXP(S1234a + S1234b) IF (do_periodic) THEN temp = P - Q PQ = pbc(temp, cell) shift = -PQ + temp temp_CC = rc + shift temp_DD = rd + shift END IF DO k = 1, SIZE(neighbor_cells) IF (do_periodic) THEN C11 = temp_CC + neighbor_cells(k)%cell_r(:) tmp_D = temp_DD + neighbor_cells(k)%cell_r(:) ELSE C11 = rc tmp_D = rd END IF Q = (Zeta_C*C11 + Zeta_D*tmp_D)*EtaInv rpq2 = (P(1) - Q(1))**2 + (P(2) - Q(2))**2 + (P(3) - Q(3))**2 IF (potential_parameter%potential_type == do_potential_truncated .OR. & potential_parameter%potential_type == do_potential_short .OR. & potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) CYCLE END IF IF (potential_parameter%potential_type == do_potential_TShPSC) THEN IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius*2.0_dp)**2) CYCLE END IF nproducts = nproducts + 1 ! allocate size as needed, ! updating the global size estimate to make this a rare event in longer simulations IF (nproducts > SIZE(product_list)) THEN !$OMP ATOMIC READ tmp_i4 = pgf_product_list_size tmp_i4 = MAX(pgf_product_list_size, (3*nproducts + 1)/2) !$OMP ATOMIC WRITE pgf_product_list_size = tmp_i4 ALLOCATE (tmp_product_list(SIZE(product_list))) tmp_product_list(:) = product_list DEALLOCATE (product_list) ALLOCATE (product_list(tmp_i4)) product_list(1:SIZE(tmp_product_list)) = tmp_product_list DEALLOCATE (tmp_product_list) ENDIF T = Rho*rpq2 SELECT CASE (potential_parameter%potential_type) CASE (do_potential_truncated) R = potential_parameter%cutoff_radius*SQRT(Rho) CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max) IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1)) factor = 2.0_dp*Pi*RhoInv CASE (do_potential_TShPSC) R = potential_parameter%cutoff_radius*SQRT(Rho) product_list(nproducts)%Fm = 0.0_dp CALL trunc_CS_poly_n20(product_list(nproducts)%Fm(1), R, T, m_max) factor = 2.0_dp*Pi*RhoInv CASE (do_potential_coulomb) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1)) factor = 2.0_dp*Pi*RhoInv CASE (do_potential_short) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1)) omega2 = potential_parameter%omega**2 omega_corr2 = omega2/(omega2 + Rho) omega_corr = SQRT(omega_corr2) T = T*omega_corr2 CALL fgamma(m_max, T, Fm) tmp = -omega_corr DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp tmp = tmp*omega_corr2 END DO factor = 2.0_dp*Pi*RhoInv CASE (do_potential_long) omega2 = potential_parameter%omega**2 omega_corr2 = omega2/(omega2 + Rho) omega_corr = SQRT(omega_corr2) T = T*omega_corr2 CALL fgamma(m_max, T, product_list(nproducts)%Fm(1)) tmp = omega_corr DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*tmp tmp = tmp*omega_corr2 END DO factor = 2.0_dp*Pi*RhoInv CASE (do_potential_mix_cl) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1)) omega2 = potential_parameter%omega**2 omega_corr2 = omega2/(omega2 + Rho) omega_corr = SQRT(omega_corr2) T = T*omega_corr2 CALL fgamma(m_max, T, Fm) tmp = omega_corr DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = & product_list(nproducts)%Fm(l)*potential_parameter%scale_coulomb & + Fm(l)*tmp*potential_parameter%scale_longrange tmp = tmp*omega_corr2 END DO factor = 2.0_dp*Pi*RhoInv CASE (do_potential_mix_cl_trunc) ! truncated R = potential_parameter%cutoff_radius*SQRT(rho) CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max) IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1)) ! Coulomb CALL fgamma(m_max, T, Fm) DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)* & (potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - & Fm(l)*potential_parameter%scale_longrange ENDDO ! longrange omega2 = potential_parameter%omega**2 omega_corr2 = omega2/(omega2 + Rho) omega_corr = SQRT(omega_corr2) T = T*omega_corr2 CALL fgamma(m_max, T, Fm) tmp = omega_corr DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange tmp = tmp*omega_corr2 END DO factor = 2.0_dp*Pi*RhoInv CASE (do_potential_gaussian) omega2 = potential_parameter%omega**2 T = -omega2*T/(Rho + omega2) tmp = 1.0_dp DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = EXP(T)*tmp tmp = tmp*omega2/(Rho + omega2) END DO factor = (Pi/(Rho + omega2))**(1.5_dp) CASE (do_potential_mix_lg) omega2 = potential_parameter%omega**2 omega_corr2 = omega2/(omega2 + Rho) omega_corr = SQRT(omega_corr2) T = T*omega_corr2 CALL fgamma(m_max, T, Fm) tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange DO l = 1, m_max + 1 Fm(l) = Fm(l)*tmp tmp = tmp*omega_corr2 END DO T = Rho*rpq2 T = -omega2*T/(Rho + omega2) tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = EXP(T)*tmp + Fm(l) tmp = tmp*omega2/(Rho + omega2) END DO factor = 1.0_dp CASE (do_potential_id) num = list1%zeta*list1%zetb den = list1%zeta + list1%zetb ssss = -num/den*SUM((ra - rb)**2) num = den*Zeta_C den = den + Zeta_C ssss = ssss - num/den*SUM((P - rc)**2) G(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + Zeta_C*rc(:))/den num = den*Zeta_D den = den + Zeta_D ssss = ssss - num/den*SUM((G - rd)**2) product_list(nproducts)%Fm(:) = EXP(ssss) factor = 1.0_dp IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234 END SELECT tmp = (Pi*ZetapEtaInv)**3 factor = factor*S1234*SQRT(tmp) DO l = 1, m_max + 1 product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*factor END DO W = (Zeta1*P + Eta*Q)*ZetapEtaInv product_list(nproducts)%ra = ra product_list(nproducts)%rb = rb product_list(nproducts)%rc = C11 product_list(nproducts)%rd = tmp_D product_list(nproducts)%ZetapEtaInv = ZetapEtaInv product_list(nproducts)%Rho = Rho product_list(nproducts)%RhoInv = RhoInv product_list(nproducts)%P = P product_list(nproducts)%Q = Q product_list(nproducts)%W = W product_list(nproducts)%AB = ra - rb product_list(nproducts)%CD = C11 - tmp_D END DO END DO END DO END SUBROUTINE build_pgf_product_list ! ************************************************************************************************** !> \brief ... !> \param npgfa ... !> \param npgfb ... !> \param list ... !> \param zeta ... !> \param zetb ... !> \param screen1 ... !> \param screen2 ... !> \param pgf ... !> \param R_pgf ... !> \param log10_pmax ... !> \param log10_eps_schwarz ... !> \param ra ... !> \param rb ... !> \param nelements ... !> \param neighbor_cells ... !> \param nimages ... !> \param do_periodic ... ! ************************************************************************************************** SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, & log10_pmax, log10_eps_schwarz, ra, rb, nelements, & neighbor_cells, nimages, do_periodic) INTEGER, INTENT(IN) :: npgfa, npgfb TYPE(hfx_pgf_list), DIMENSION(npgfa*npgfb) :: list REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb REAL(dp), INTENT(IN) :: screen1(2), screen2(2) TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & POINTER :: pgf, R_pgf REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz, ra(3), & rb(3) INTEGER, INTENT(OUT) :: nelements TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells INTEGER :: nimages(npgfa*npgfb) LOGICAL, INTENT(IN) :: do_periodic INTEGER :: element_counter, i, ipgf, j, jpgf REAL(dp) :: AB(3), im_B(3), pgf_max, rab2, Zeta1, & Zeta_A, Zeta_B, ZetaInv nimages = 0 ! ** inner loop may never be reached nelements = npgfa*npgfb DO i = 1, SIZE(neighbor_cells) IF (do_periodic) THEN im_B = rb + neighbor_cells(i)%cell_r(:) ELSE im_B = rb END IF AB = ra - im_B rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2 IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE element_counter = 0 DO ipgf = 1, npgfa DO jpgf = 1, npgfb element_counter = element_counter + 1 pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2) IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN CYCLE END IF nimages(element_counter) = nimages(element_counter) + 1 list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max list(element_counter)%image_list(nimages(element_counter))%ra = ra list(element_counter)%image_list(nimages(element_counter))%rb = im_B list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2 Zeta_A = zeta(ipgf) Zeta_B = zetb(jpgf) Zeta1 = Zeta_A + Zeta_B ZetaInv = 1.0_dp/Zeta1 IF (nimages(element_counter) == 1) THEN list(element_counter)%ipgf = ipgf list(element_counter)%jpgf = jpgf list(element_counter)%zetaInv = ZetaInv list(element_counter)%zetapzetb = Zeta1 list(element_counter)%zeta = Zeta_A list(element_counter)%zetb = Zeta_B END IF list(element_counter)%image_list(nimages(element_counter))%S1234 = (-Zeta_A*Zeta_B*ZetaInv*rab2) list(element_counter)%image_list(nimages(element_counter))%P = (Zeta_A*ra + Zeta_B*im_B)*ZetaInv list(element_counter)%image_list(nimages(element_counter))%R = & MAX(0.0_dp, R_pgf(jpgf, ipgf)%x(1)*rab2 + R_pgf(jpgf, ipgf)%x(2)) list(element_counter)%image_list(nimages(element_counter))%ra = ra list(element_counter)%image_list(nimages(element_counter))%rb = im_B list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2 list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell END DO END DO nelements = MAX(nelements, element_counter) END DO DO j = 1, nelements list(j)%nimages = nimages(j) END DO ! ** Remove unused elements element_counter = 0 DO j = 1, nelements IF (list(j)%nimages == 0) CYCLE element_counter = element_counter + 1 list(element_counter)%nimages = list(j)%nimages list(element_counter)%zetapzetb = list(j)%zetapzetb list(element_counter)%ZetaInv = list(j)%ZetaInv list(element_counter)%zeta = list(j)%zeta list(element_counter)%zetb = list(j)%zetb list(element_counter)%ipgf = list(j)%ipgf list(element_counter)%jpgf = list(j)%jpgf DO i = 1, list(j)%nimages list(element_counter)%image_list(i) = list(j)%image_list(i) END DO END DO nelements = element_counter END SUBROUTINE build_pair_list_pgf ! ************************************************************************************************** !> \brief ... !> \param natom ... !> \param list ... !> \param set_list ... !> \param i_start ... !> \param i_end ... !> \param j_start ... !> \param j_end ... !> \param kind_of ... !> \param basis_parameter ... !> \param particle_set ... !> \param do_periodic ... !> \param coeffs_set ... !> \param coeffs_kind ... !> \param coeffs_kind_max0 ... !> \param log10_eps_schwarz ... !> \param cell ... !> \param pmax_blocks ... !> \param atomic_pair_list ... ! ************************************************************************************************** SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, & do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, & pmax_blocks, atomic_pair_list) INTEGER, INTENT(IN) :: natom TYPE(pair_list_type), INTENT(OUT) :: list TYPE(pair_set_list_type), DIMENSION(*), & INTENT(OUT) :: set_list INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end INTEGER :: kind_of(*) TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter TYPE(particle_type), DIMENSION(:), POINTER :: particle_set LOGICAL, INTENT(IN) :: do_periodic TYPE(hfx_screen_coeff_type), & DIMENSION(:, :, :, :), POINTER :: coeffs_set TYPE(hfx_screen_coeff_type), DIMENSION(:, :) :: coeffs_kind REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz TYPE(cell_type), POINTER :: cell REAL(dp) :: pmax_blocks LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list INTEGER :: iatom, ikind, iset, jatom, jkind, jset, & n_element, nset_ij, nseta, nsetb REAL(KIND=dp) :: rab2 REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp n_element = 0 nset_ij = 0 DO iatom = i_start, i_end DO jatom = j_start, j_end IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE ikind = kind_of(iatom) nseta = basis_parameter(ikind)%nset ra = particle_set(iatom)%r(:) IF (jatom < iatom) CYCLE jkind = kind_of(jatom) nsetb = basis_parameter(jkind)%nset rb = particle_set(jatom)%r(:) IF (do_periodic) THEN temp = rb - ra pbc_B = pbc(temp, cell) B11 = ra + pbc_B rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2 ELSE rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2 B11 = rb ! ra - rb END IF IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + & coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE n_element = n_element + 1 list%elements(n_element)%pair(1) = iatom list%elements(n_element)%pair(2) = jatom list%elements(n_element)%kind_pair(1) = ikind list%elements(n_element)%kind_pair(2) = jkind list%elements(n_element)%r1 = ra list%elements(n_element)%r2 = B11 list%elements(n_element)%dist2 = rab2 ! build a list of guaranteed overlapping sets list%elements(n_element)%set_bounds(1) = nset_ij + 1 DO iset = 1, nseta DO jset = 1, nsetb IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + & coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE nset_ij = nset_ij + 1 set_list(nset_ij)%pair(1) = iset set_list(nset_ij)%pair(2) = jset END DO END DO list%elements(n_element)%set_bounds(2) = nset_ij END DO END DO list%n_element = n_element END SUBROUTINE build_pair_list ! ************************************************************************************************** !> \brief ... !> \param natom ... !> \param atomic_pair_list ... !> \param kind_of ... !> \param basis_parameter ... !> \param particle_set ... !> \param do_periodic ... !> \param coeffs_kind ... !> \param coeffs_kind_max0 ... !> \param log10_eps_schwarz ... !> \param cell ... !> \param blocks ... ! ************************************************************************************************** SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, & do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, & blocks) INTEGER, INTENT(IN) :: natom LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list INTEGER :: kind_of(*) TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter TYPE(particle_type), DIMENSION(:), POINTER :: particle_set LOGICAL, INTENT(IN) :: do_periodic TYPE(hfx_screen_coeff_type), DIMENSION(:, :) :: coeffs_kind REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz TYPE(cell_type), POINTER :: cell TYPE(hfx_block_range_type), DIMENSION(:), POINTER :: blocks INTEGER :: iatom, iatom_end, iatom_start, iblock, & ikind, jatom, jatom_end, jatom_start, & jblock, jkind, nseta, nsetb REAL(KIND=dp) :: rab2 REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp atomic_pair_list = .FALSE. DO iblock = 1, SIZE(blocks) iatom_start = blocks(iblock)%istart iatom_end = blocks(iblock)%iend DO jblock = 1, SIZE(blocks) jatom_start = blocks(jblock)%istart jatom_end = blocks(jblock)%iend DO iatom = iatom_start, iatom_end ikind = kind_of(iatom) nseta = basis_parameter(ikind)%nset ra = particle_set(iatom)%r(:) DO jatom = jatom_start, jatom_end IF (jatom < iatom) CYCLE jkind = kind_of(jatom) nsetb = basis_parameter(jkind)%nset rb = particle_set(jatom)%r(:) IF (do_periodic) THEN temp = rb - ra pbc_B = pbc(temp, cell) B11 = ra + pbc_B rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2 ELSE rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2 B11 = rb ! ra - rb END IF IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + & coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) CYCLE atomic_pair_list(jatom, iatom) = .TRUE. atomic_pair_list(iatom, jatom) = .TRUE. END DO END DO END DO END DO END SUBROUTINE build_atomic_pair_list ! ************************************************************************************************** !> \brief ... !> \param natom ... !> \param list ... !> \param set_list ... !> \param i_start ... !> \param i_end ... !> \param j_start ... !> \param j_end ... !> \param kind_of ... !> \param basis_parameter ... !> \param particle_set ... !> \param do_periodic ... !> \param coeffs_set ... !> \param coeffs_kind ... !> \param coeffs_kind_max0 ... !> \param log10_eps_schwarz ... !> \param cell ... !> \param pmax_blocks ... !> \param atomic_pair_list ... !> \param skip_atom_symmetry ... ! ************************************************************************************************** SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, & do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, & pmax_blocks, atomic_pair_list, skip_atom_symmetry) INTEGER, INTENT(IN) :: natom TYPE(pair_list_type_mp2) :: list TYPE(pair_set_list_type), DIMENSION(*), & INTENT(OUT) :: set_list INTEGER, INTENT(IN) :: i_start, i_end, j_start, j_end INTEGER :: kind_of(*) TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter TYPE(particle_type), DIMENSION(:), POINTER :: particle_set LOGICAL, INTENT(IN) :: do_periodic TYPE(hfx_screen_coeff_type), & DIMENSION(:, :, :, :), POINTER :: coeffs_set TYPE(hfx_screen_coeff_type), DIMENSION(:, :) :: coeffs_kind REAL(KIND=dp), INTENT(IN) :: coeffs_kind_max0, log10_eps_schwarz TYPE(cell_type), POINTER :: cell REAL(dp) :: pmax_blocks LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list LOGICAL, OPTIONAL :: skip_atom_symmetry INTEGER :: iatom, ikind, iset, jatom, jkind, jset, & n_element, nset_ij, nseta, nsetb LOGICAL :: my_skip_atom_symmetry REAL(KIND=dp) :: rab2 REAL(KIND=dp), DIMENSION(3) :: B11, pbc_B, ra, rb, temp n_element = 0 nset_ij = 0 my_skip_atom_symmetry = .FALSE. IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry DO iatom = i_start, i_end DO jatom = j_start, j_end IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE ikind = kind_of(iatom) nseta = basis_parameter(ikind)%nset ra = particle_set(iatom)%r(:) IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) CYCLE jkind = kind_of(jatom) nsetb = basis_parameter(jkind)%nset rb = particle_set(jatom)%r(:) IF (do_periodic) THEN temp = rb - ra pbc_B = pbc(temp, cell) B11 = ra + pbc_B rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2 ELSE rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2 B11 = rb ! ra - rb END IF IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + & coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE n_element = n_element + 1 list%elements(n_element)%pair(1) = iatom list%elements(n_element)%pair(2) = jatom list%elements(n_element)%kind_pair(1) = ikind list%elements(n_element)%kind_pair(2) = jkind list%elements(n_element)%r1 = ra list%elements(n_element)%r2 = B11 list%elements(n_element)%dist2 = rab2 ! build a list of guaranteed overlapping sets list%elements(n_element)%set_bounds(1) = nset_ij + 1 DO iset = 1, nseta DO jset = 1, nsetb IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + & coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE nset_ij = nset_ij + 1 set_list(nset_ij)%pair(1) = iset set_list(nset_ij)%pair(2) = jset END DO END DO list%elements(n_element)%set_bounds(2) = nset_ij END DO END DO list%n_element = n_element END SUBROUTINE build_pair_list_mp2 END MODULE hfx_pair_list_methods