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