1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Interface to the Libint-Library 8!> \par History 9!> 11.2006 created [Manuel Guidon] 10!> 11.2019 Fixed potential_id initial values (A. Bussy) 11!> \author Manuel Guidon 12! ************************************************************************************************** 13MODULE hfx_libint_interface 14 15 USE cell_types, ONLY: cell_type,& 16 real_to_scaled 17 USE gamma, ONLY: fgamma => fgamma_0 18 USE hfx_contraction_methods, ONLY: contract 19 USE hfx_types, ONLY: hfx_pgf_product_list,& 20 hfx_potential_type 21 USE input_constants, ONLY: & 22 do_potential_coulomb, do_potential_gaussian, do_potential_id, do_potential_long, & 23 do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, & 24 do_potential_truncated 25 USE kinds, ONLY: dp,& 26 int_8 27 USE libint_wrapper, ONLY: cp_libint_get_derivs,& 28 cp_libint_get_eris,& 29 cp_libint_set_params_eri,& 30 cp_libint_set_params_eri_deriv,& 31 cp_libint_set_params_eri_screen,& 32 cp_libint_t,& 33 get_ssss_f_val,& 34 prim_data_f_size 35 USE mathconstants, ONLY: pi 36 USE orbital_pointers, ONLY: nco 37 USE t_c_g0, ONLY: t_c_g0_n 38#include "./base/base_uses.f90" 39 40 IMPLICIT NONE 41 PRIVATE 42 PUBLIC :: evaluate_eri, & 43 evaluate_deriv_eri, & 44 evaluate_eri_screen 45 46 INTEGER, DIMENSION(12), PARAMETER :: full_perm1 = (/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12/) 47 INTEGER, DIMENSION(12), PARAMETER :: full_perm2 = (/4, 5, 6, 1, 2, 3, 7, 8, 9, 10, 11, 12/) 48 INTEGER, DIMENSION(12), PARAMETER :: full_perm3 = (/1, 2, 3, 4, 5, 6, 10, 11, 12, 7, 8, 9/) 49 INTEGER, DIMENSION(12), PARAMETER :: full_perm4 = (/4, 5, 6, 1, 2, 3, 10, 11, 12, 7, 8, 9/) 50 INTEGER, DIMENSION(12), PARAMETER :: full_perm5 = (/7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6/) 51 INTEGER, DIMENSION(12), PARAMETER :: full_perm6 = (/7, 8, 9, 10, 11, 12, 4, 5, 6, 1, 2, 3/) 52 INTEGER, DIMENSION(12), PARAMETER :: full_perm7 = (/10, 11, 12, 7, 8, 9, 1, 2, 3, 4, 5, 6/) 53 INTEGER, DIMENSION(12), PARAMETER :: full_perm8 = (/10, 11, 12, 7, 8, 9, 4, 5, 6, 1, 2, 3/) 54 55!*** 56 57CONTAINS 58 59! ************************************************************************************************** 60!> \brief Fill data structure used in libint 61!> \param A ... 62!> \param B ... 63!> \param C ... 64!> \param D ... 65!> \param Zeta_A ... 66!> \param Zeta_B ... 67!> \param Zeta_C ... 68!> \param Zeta_D ... 69!> \param m_max ... 70!> \param potential_parameter ... 71!> \param libint ... 72!> \param R11 ... 73!> \param R22 ... 74!> \par History 75!> 03.2007 created [Manuel Guidon] 76!> \author Manuel Guidon 77! ************************************************************************************************** 78 SUBROUTINE build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, & 79 potential_parameter, libint, R11, R22) 80 81 REAL(KIND=dp) :: A(3), B(3), C(3), D(3) 82 REAL(KIND=dp), INTENT(IN) :: Zeta_A, Zeta_B, Zeta_C, Zeta_D 83 INTEGER, INTENT(IN) :: m_max 84 TYPE(hfx_potential_type) :: potential_parameter 85 TYPE(cp_libint_t) :: libint 86 REAL(dp) :: R11, R22 87 88 INTEGER :: i 89 LOGICAL :: use_gamma 90 REAL(KIND=dp) :: AB(3), AB2, CD(3), CD2, den, Eta, EtaInv, factor, G(3), num, omega2, & 91 omega_corr, omega_corr2, P(3), PQ(3), PQ2, Q(3), R, R1, R2, Rho, RhoInv, S1234, ssss, T, & 92 tmp, W(3), Zeta, ZetaInv, ZetapEtaInv 93 REAL(KIND=dp), DIMENSION(prim_data_f_size) :: F, Fm 94 95 Zeta = Zeta_A + Zeta_B 96 ZetaInv = 1.0_dp/Zeta 97 Eta = Zeta_C + Zeta_D 98 EtaInv = 1.0_dp/Eta 99 ZetapEtaInv = Zeta + Eta 100 ZetapEtaInv = 1.0_dp/ZetapEtaInv 101 Rho = Zeta*Eta*ZetapEtaInv 102 RhoInv = 1.0_dp/Rho 103 104 DO i = 1, 3 105 P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv 106 Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv 107 AB(i) = A(i) - B(i) 108 CD(i) = C(i) - D(i) 109 PQ(i) = P(i) - Q(i) 110 W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv 111 END DO 112 113 AB2 = DOT_PRODUCT(AB, AB) 114 CD2 = DOT_PRODUCT(CD, CD) 115 PQ2 = DOT_PRODUCT(PQ, PQ) 116 117 S1234 = EXP((-Zeta_A*Zeta_B*ZetaInv*AB2) + (-Zeta_C*Zeta_D*EtaInv*CD2)) 118 T = Rho*PQ2 119 120 SELECT CASE (potential_parameter%potential_type) 121 CASE (do_potential_truncated) 122 R = potential_parameter%cutoff_radius*SQRT(rho) 123 R1 = R11 124 R2 = R22 125 IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN 126 RETURN 127 END IF 128 CALL t_c_g0_n(F(1), use_gamma, R, T, m_max) 129 IF (use_gamma) CALL fgamma(m_max, T, F(1)) 130 factor = 2.0_dp*Pi*RhoInv 131 CASE (do_potential_coulomb) 132 CALL fgamma(m_max, T, F(1)) 133 factor = 2.0_dp*Pi*RhoInv 134 CASE (do_potential_short) 135 R = potential_parameter%cutoff_radius*SQRT(rho) 136 R1 = R11 137 R2 = R22 138 IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN 139 RETURN 140 END IF 141 CALL fgamma(m_max, T, F(1)) 142 omega2 = potential_parameter%omega**2 143 omega_corr2 = omega2/(omega2 + Rho) 144 omega_corr = SQRT(omega_corr2) 145 T = T*omega_corr2 146 CALL fgamma(m_max, T, Fm) 147 tmp = -omega_corr 148 DO i = 1, m_max + 1 149 F(i) = F(i) + Fm(i)*tmp 150 tmp = tmp*omega_corr2 151 END DO 152 factor = 2.0_dp*Pi*RhoInv 153 CASE (do_potential_long) 154 omega2 = potential_parameter%omega**2 155 omega_corr2 = omega2/(omega2 + Rho) 156 omega_corr = SQRT(omega_corr2) 157 T = T*omega_corr2 158 CALL fgamma(m_max, T, F(1)) 159 tmp = omega_corr 160 DO i = 1, m_max + 1 161 F(i) = F(i)*tmp 162 tmp = tmp*omega_corr2 163 END DO 164 factor = 2.0_dp*Pi*RhoInv 165 CASE (do_potential_mix_cl) 166 CALL fgamma(m_max, T, F(1)) 167 omega2 = potential_parameter%omega**2 168 omega_corr2 = omega2/(omega2 + Rho) 169 omega_corr = SQRT(omega_corr2) 170 T = T*omega_corr2 171 CALL fgamma(m_max, T, Fm) 172 tmp = omega_corr 173 DO i = 1, m_max + 1 174 F(i) = F(i)*potential_parameter%scale_coulomb + Fm(i)*tmp*potential_parameter%scale_longrange 175 tmp = tmp*omega_corr2 176 END DO 177 factor = 2.0_dp*Pi*RhoInv 178 CASE (do_potential_mix_cl_trunc) 179 180 ! truncated 181 R = potential_parameter%cutoff_radius*SQRT(rho) 182 R1 = R11 183 R2 = R22 184 IF (PQ2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) THEN 185 RETURN 186 END IF 187 CALL t_c_g0_n(F(1), use_gamma, R, T, m_max) 188 IF (use_gamma) CALL fgamma(m_max, T, F(1)) 189 190 ! Coulomb 191 CALL fgamma(m_max, T, Fm) 192 193 DO i = 1, m_max + 1 194 F(i) = F(i)*(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - & 195 Fm(i)*potential_parameter%scale_longrange 196 ENDDO 197 198 ! longrange 199 omega2 = potential_parameter%omega**2 200 omega_corr2 = omega2/(omega2 + Rho) 201 omega_corr = SQRT(omega_corr2) 202 T = T*omega_corr2 203 CALL fgamma(m_max, T, Fm) 204 tmp = omega_corr 205 DO i = 1, m_max + 1 206 F(i) = F(i) + Fm(i)*tmp*potential_parameter%scale_longrange 207 tmp = tmp*omega_corr2 208 END DO 209 factor = 2.0_dp*Pi*RhoInv 210 211 CASE (do_potential_gaussian) 212 omega2 = potential_parameter%omega**2 213 T = -omega2*T/(Rho + omega2) 214 tmp = 1.0_dp 215 DO i = 1, m_max + 1 216 F(i) = EXP(T)*tmp 217 tmp = tmp*omega2/(Rho + omega2) 218 END DO 219 factor = (Pi/(Rho + omega2))**(1.5_dp) 220 CASE (do_potential_mix_lg) 221 omega2 = potential_parameter%omega**2 222 omega_corr2 = omega2/(omega2 + Rho) 223 omega_corr = SQRT(omega_corr2) 224 T = T*omega_corr2 225 CALL fgamma(m_max, T, Fm) 226 tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange 227 DO i = 1, m_max + 1 228 Fm(i) = Fm(i)*tmp 229 tmp = tmp*omega_corr2 230 END DO 231 T = Rho*PQ2 232 T = -omega2*T/(Rho + omega2) 233 tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian 234 DO i = 1, m_max + 1 235 F(i) = EXP(T)*tmp + Fm(i) 236 tmp = tmp*omega2/(Rho + omega2) 237 END DO 238 factor = 1.0_dp 239 CASE (do_potential_id) 240 ssss = -Zeta_A*Zeta_B*ZetaInv*AB2 241 242 num = (Zeta_A + Zeta_B)*Zeta_C 243 den = Zeta_A + Zeta_B + Zeta_C 244 ssss = ssss - num/den*SUM((P - C)**2) 245 246 G(:) = (Zeta_A*A(:) + Zeta_B*B(:) + Zeta_C*C(:))/den 247 num = den*Zeta_D 248 den = den + Zeta_D 249 ssss = ssss - num/den*SUM((G - D)**2) 250 251 F(:) = EXP(ssss) 252 factor = 1.0_dp 253 IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234 254 END SELECT 255 256 tmp = (Pi*ZetapEtaInv)**3 257 factor = factor*S1234*SQRT(tmp) 258 259 DO i = 1, m_max + 1 260 F(i) = F(i)*factor 261 ENDDO 262 263 CALL cp_libint_set_params_eri_screen(libint, A, B, C, D, P, Q, W, ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F) 264 265 END SUBROUTINE build_quartet_data_screen 266 267! ************************************************************************************************** 268!> \brief Fill data structure used in libderiv 269!> \param libint ... 270!> \param A ... 271!> \param B ... 272!> \param C ... 273!> \param D ... 274!> \param Zeta_A ... 275!> \param Zeta_B ... 276!> \param Zeta_C ... 277!> \param Zeta_D ... 278!> \param ZetaInv ... 279!> \param EtaInv ... 280!> \param ZetapEtaInv ... 281!> \param Rho ... 282!> \param RhoInv ... 283!> \param P ... 284!> \param Q ... 285!> \param W ... 286!> \param m_max ... 287!> \param F ... 288!> \par History 289!> 03.2007 created [Manuel Guidon] 290!> \author Manuel Guidon 291! ************************************************************************************************** 292 SUBROUTINE build_deriv_data(libint, A, B, C, D, & 293 Zeta_A, Zeta_B, Zeta_C, Zeta_D, & 294 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, & 295 P, Q, W, m_max, F) 296 297 TYPE(cp_libint_t) :: libint 298 REAL(KIND=dp), INTENT(IN) :: A(3), B(3), C(3), D(3) 299 REAL(dp), INTENT(IN) :: Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, & 300 EtaInv, ZetapEtaInv, Rho, RhoInv, & 301 P(3), Q(3), W(3) 302 INTEGER, INTENT(IN) :: m_max 303 REAL(KIND=dp), DIMENSION(:) :: F 304 305 MARK_USED(D) ! todo: fix 306 MARK_USED(Zeta_C) 307 MARK_USED(RhoInv) 308 309 CALL cp_libint_set_params_eri_deriv(libint, A, B, C, D, P, Q, W, zeta_A, zeta_B, zeta_C, zeta_D, & 310 ZetaInv, EtaInv, ZetapEtaInv, Rho, m_max, F) 311 312 END SUBROUTINE build_deriv_data 313 314! ************************************************************************************************** 315!> \brief Evaluates derivatives of electron repulsion integrals for a primitive quartet 316!> \param deriv ... 317!> \param nproducts ... 318!> \param pgf_product_list ... 319!> \param n_a ... 320!> \param n_b ... 321!> \param n_c ... 322!> \param n_d ... 323!> \param ncoa ... 324!> \param ncob ... 325!> \param ncoc ... 326!> \param ncod ... 327!> \param nsgfa ... 328!> \param nsgfb ... 329!> \param nsgfc ... 330!> \param nsgfd ... 331!> \param primitives ... 332!> \param max_contraction ... 333!> \param tmp_max_all ... 334!> \param eps_schwarz ... 335!> \param neris ... 336!> \param Zeta_A ... 337!> \param Zeta_B ... 338!> \param Zeta_C ... 339!> \param Zeta_D ... 340!> \param ZetaInv ... 341!> \param EtaInv ... 342!> \param s_offset_a ... 343!> \param s_offset_b ... 344!> \param s_offset_c ... 345!> \param s_offset_d ... 346!> \param nl_a ... 347!> \param nl_b ... 348!> \param nl_c ... 349!> \param nl_d ... 350!> \param nsoa ... 351!> \param nsob ... 352!> \param nsoc ... 353!> \param nsod ... 354!> \param sphi_a ... 355!> \param sphi_b ... 356!> \param sphi_c ... 357!> \param sphi_d ... 358!> \param work ... 359!> \param work2 ... 360!> \param work_forces ... 361!> \param buffer1 ... 362!> \param buffer2 ... 363!> \param primitives_tmp ... 364!> \param use_virial ... 365!> \param work_virial ... 366!> \param work2_virial ... 367!> \param primitives_tmp_virial ... 368!> \param primitives_virial ... 369!> \param cell ... 370!> \param tmp_max_all_virial ... 371!> \par History 372!> 03.2007 created [Manuel Guidon] 373!> 08.2007 refactured permutation part [Manuel Guidon] 374!> \author Manuel Guidon 375! ************************************************************************************************** 376 SUBROUTINE evaluate_deriv_eri(deriv, nproducts, pgf_product_list, & 377 n_a, n_b, n_c, n_d, & 378 ncoa, ncob, ncoc, ncod, & 379 nsgfa, nsgfb, nsgfc, nsgfd, & 380 primitives, max_contraction, tmp_max_all, & 381 eps_schwarz, neris, & 382 Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, EtaInv, & 383 s_offset_a, s_offset_b, s_offset_c, s_offset_d, & 384 nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, & 385 sphi_a, sphi_b, sphi_c, sphi_d, & 386 work, work2, work_forces, buffer1, buffer2, primitives_tmp, & 387 use_virial, work_virial, work2_virial, primitives_tmp_virial, & 388 primitives_virial, cell, tmp_max_all_virial) 389 390 TYPE(cp_libint_t) :: deriv 391 INTEGER, INTENT(IN) :: nproducts 392 TYPE(hfx_pgf_product_list), DIMENSION(nproducts) :: pgf_product_list 393 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, & 394 ncod, nsgfa, nsgfb, nsgfc, nsgfd 395 REAL(dp), & 396 DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12) :: primitives 397 REAL(dp) :: max_contraction, tmp_max_all, eps_schwarz 398 INTEGER(int_8) :: neris 399 REAL(dp), INTENT(IN) :: Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, & 400 EtaInv 401 INTEGER :: s_offset_a, s_offset_b, s_offset_c, & 402 s_offset_d, nl_a, nl_b, nl_c, nl_d, & 403 nsoa, nsob, nsoc, nsod 404 REAL(dp), DIMENSION(ncoa, nsoa*nl_a) :: sphi_a 405 REAL(dp), DIMENSION(ncob, nsob*nl_b) :: sphi_b 406 REAL(dp), DIMENSION(ncoc, nsoc*nl_c) :: sphi_c 407 REAL(dp), DIMENSION(ncod, nsod*nl_d) :: sphi_d 408 REAL(dp) :: work(ncoa*ncob*ncoc*ncod, 12), work2(ncoa, ncob, ncoc, ncod, 12), & 409 work_forces(ncoa*ncob*ncoc*ncod, 12) 410 REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2 411 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*& 412 nl_c, nsod*nl_d) :: primitives_tmp 413 LOGICAL, INTENT(IN) :: use_virial 414 REAL(dp) :: work_virial(ncoa*ncob*ncoc*ncod, 12, 3), & 415 work2_virial(ncoa, ncob, ncoc, ncod, 12, 3) 416 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*& 417 nl_c, nsod*nl_d) :: primitives_tmp_virial 418 REAL(dp), & 419 DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12, 3) :: primitives_virial 420 TYPE(cell_type), POINTER :: cell 421 REAL(dp) :: tmp_max_all_virial 422 423 INTEGER :: a_mysize(1), i, j, k, l, m, m_max, & 424 mysize, n, p1, p2, p3, p4, perm_case 425 REAL(dp) :: A(3), AB(3), B(3), C(3), CD(3), D(3), & 426 P(3), Q(3), Rho, RhoInv, scoord(12), & 427 tmp_max, tmp_max_virial, W(3), & 428 ZetapEtaInv 429 430 m_max = n_a + n_b + n_c + n_d 431 m_max = m_max + 1 432 mysize = ncoa*ncob*ncoc*ncod 433 a_mysize = mysize 434 435 work = 0.0_dp 436 work2 = 0.0_dp 437 438 IF (use_virial) THEN 439 work_virial = 0.0_dp 440 work2_virial = 0.0_dp 441 END IF 442 443 perm_case = 1 444 IF (n_a < n_b) THEN 445 perm_case = perm_case + 1 446 END IF 447 IF (n_c < n_d) THEN 448 perm_case = perm_case + 2 449 END IF 450 IF (n_a + n_b > n_c + n_d) THEN 451 perm_case = perm_case + 4 452 END IF 453 454 SELECT CASE (perm_case) 455 CASE (1) 456 DO i = 1, nproducts 457 A = pgf_product_list(i)%ra 458 B = pgf_product_list(i)%rb 459 C = pgf_product_list(i)%rc 460 D = pgf_product_list(i)%rd 461 Rho = pgf_product_list(i)%Rho 462 RhoInv = pgf_product_list(i)%RhoInv 463 P = pgf_product_list(i)%P 464 Q = pgf_product_list(i)%Q 465 W = pgf_product_list(i)%W 466 AB = pgf_product_list(i)%AB 467 CD = pgf_product_list(i)%CD 468 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 469 470 CALL build_deriv_data(deriv, A, B, C, D, & 471 Zeta_A, Zeta_B, Zeta_C, Zeta_D, & 472 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, & 473 P, Q, W, m_max, pgf_product_list(i)%Fm) 474 CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize) 475 DO k = 4, 6 476 DO j = 1, mysize 477 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 478 work_forces(j, k + 3) + & 479 work_forces(j, k + 6)) 480 END DO 481 END DO 482 DO k = 1, 12 483 DO j = 1, mysize 484 work(j, k) = work(j, k) + work_forces(j, k) 485 END DO 486 END DO 487 neris = neris + 12*mysize 488 IF (use_virial) THEN 489 CALL real_to_scaled(scoord(1:3), A, cell) 490 CALL real_to_scaled(scoord(4:6), B, cell) 491 CALL real_to_scaled(scoord(7:9), C, cell) 492 CALL real_to_scaled(scoord(10:12), D, cell) 493 DO k = 1, 12 494 DO j = 1, mysize 495 DO m = 1, 3 496 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 497 END DO 498 END DO 499 END DO 500 END IF 501 END DO 502 503 DO n = 1, 12 504 tmp_max = 0.0_dp 505 DO i = 1, mysize 506 tmp_max = MAX(tmp_max, ABS(work(i, n))) 507 END DO 508 tmp_max = tmp_max*max_contraction 509 tmp_max_all = MAX(tmp_max_all, tmp_max) 510 511 DO i = 1, ncoa 512 p1 = (i - 1)*ncob 513 DO j = 1, ncob 514 p2 = (p1 + j - 1)*ncoc 515 DO k = 1, ncoc 516 p3 = (p2 + k - 1)*ncod 517 DO l = 1, ncod 518 p4 = p3 + l 519 work2(i, j, k, l, full_perm1(n)) = work(p4, n) 520 END DO 521 END DO 522 END DO 523 END DO 524 END DO 525 IF (use_virial) THEN 526 DO n = 1, 12 527 tmp_max_virial = 0.0_dp 528 DO i = 1, mysize 529 tmp_max_virial = MAX(tmp_max_virial, & 530 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 531 END DO 532 tmp_max_virial = tmp_max_virial*max_contraction 533 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 534 535 DO i = 1, ncoa 536 p1 = (i - 1)*ncob 537 DO j = 1, ncob 538 p2 = (p1 + j - 1)*ncoc 539 DO k = 1, ncoc 540 p3 = (p2 + k - 1)*ncod 541 DO l = 1, ncod 542 p4 = p3 + l 543 work2_virial(i, j, k, l, full_perm1(n), 1:3) = work_virial(p4, n, 1:3) 544 END DO 545 END DO 546 END DO 547 END DO 548 END DO 549 END IF 550 CASE (2) 551 DO i = 1, nproducts 552 A = pgf_product_list(i)%ra 553 B = pgf_product_list(i)%rb 554 C = pgf_product_list(i)%rc 555 D = pgf_product_list(i)%rd 556 Rho = pgf_product_list(i)%Rho 557 RhoInv = pgf_product_list(i)%RhoInv 558 P = pgf_product_list(i)%P 559 Q = pgf_product_list(i)%Q 560 W = pgf_product_list(i)%W 561 AB = pgf_product_list(i)%AB 562 CD = pgf_product_list(i)%CD 563 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 564 565 CALL build_deriv_data(deriv, B, A, C, D, & 566 Zeta_B, Zeta_A, Zeta_C, Zeta_D, & 567 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, & 568 P, Q, W, m_max, pgf_product_list(i)%Fm) 569 CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize) 570 DO k = 4, 6 571 DO j = 1, mysize 572 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 573 work_forces(j, k + 3) + & 574 work_forces(j, k + 6)) 575 ENDDO 576 END DO 577 DO k = 1, 12 578 DO j = 1, mysize 579 work(j, k) = work(j, k) + work_forces(j, k) 580 END DO 581 END DO 582 neris = neris + 12*mysize 583 IF (use_virial) THEN 584 CALL real_to_scaled(scoord(1:3), B, cell) 585 CALL real_to_scaled(scoord(4:6), A, cell) 586 CALL real_to_scaled(scoord(7:9), C, cell) 587 CALL real_to_scaled(scoord(10:12), D, cell) 588 DO k = 1, 12 589 DO j = 1, mysize 590 DO m = 1, 3 591 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 592 END DO 593 END DO 594 END DO 595 END IF 596 597 END DO 598 DO n = 1, 12 599 tmp_max = 0.0_dp 600 DO i = 1, mysize 601 tmp_max = MAX(tmp_max, ABS(work(i, n))) 602 END DO 603 tmp_max = tmp_max*max_contraction 604 tmp_max_all = MAX(tmp_max_all, tmp_max) 605 606 DO j = 1, ncob 607 p1 = (j - 1)*ncoa 608 DO i = 1, ncoa 609 p2 = (p1 + i - 1)*ncoc 610 DO k = 1, ncoc 611 p3 = (p2 + k - 1)*ncod 612 DO l = 1, ncod 613 p4 = p3 + l 614 work2(i, j, k, l, full_perm2(n)) = work(p4, n) 615 END DO 616 END DO 617 END DO 618 END DO 619 END DO 620 IF (use_virial) THEN 621 DO n = 1, 12 622 tmp_max_virial = 0.0_dp 623 DO i = 1, mysize 624 tmp_max_virial = MAX(tmp_max_virial, & 625 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 626 END DO 627 tmp_max_virial = tmp_max_virial*max_contraction 628 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 629 630 DO j = 1, ncob 631 p1 = (j - 1)*ncoa 632 DO i = 1, ncoa 633 p2 = (p1 + i - 1)*ncoc 634 DO k = 1, ncoc 635 p3 = (p2 + k - 1)*ncod 636 DO l = 1, ncod 637 p4 = p3 + l 638 work2_virial(i, j, k, l, full_perm2(n), 1:3) = work_virial(p4, n, 1:3) 639 END DO 640 END DO 641 END DO 642 END DO 643 END DO 644 END IF 645 CASE (3) 646 DO i = 1, nproducts 647 A = pgf_product_list(i)%ra 648 B = pgf_product_list(i)%rb 649 C = pgf_product_list(i)%rc 650 D = pgf_product_list(i)%rd 651 Rho = pgf_product_list(i)%Rho 652 RhoInv = pgf_product_list(i)%RhoInv 653 P = pgf_product_list(i)%P 654 Q = pgf_product_list(i)%Q 655 W = pgf_product_list(i)%W 656 AB = pgf_product_list(i)%AB 657 CD = pgf_product_list(i)%CD 658 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 659 660 CALL build_deriv_data(deriv, A, B, D, C, & 661 Zeta_A, Zeta_B, Zeta_D, Zeta_C, & 662 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, & 663 P, Q, W, m_max, pgf_product_list(i)%Fm) 664 CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize) 665 DO k = 4, 6 666 DO j = 1, mysize 667 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 668 work_forces(j, k + 3) + & 669 work_forces(j, k + 6)) 670 END DO 671 END DO 672 DO k = 1, 12 673 DO j = 1, mysize 674 work(j, k) = work(j, k) + work_forces(j, k) 675 END DO 676 END DO 677 neris = neris + 12*mysize 678 IF (use_virial) THEN 679 CALL real_to_scaled(scoord(1:3), A, cell) 680 CALL real_to_scaled(scoord(4:6), B, cell) 681 CALL real_to_scaled(scoord(7:9), D, cell) 682 CALL real_to_scaled(scoord(10:12), C, cell) 683 DO k = 1, 12 684 DO j = 1, mysize 685 DO m = 1, 3 686 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 687 END DO 688 END DO 689 END DO 690 END IF 691 692 END DO 693 DO n = 1, 12 694 tmp_max = 0.0_dp 695 DO i = 1, mysize 696 tmp_max = MAX(tmp_max, ABS(work(i, n))) 697 END DO 698 tmp_max = tmp_max*max_contraction 699 tmp_max_all = MAX(tmp_max_all, tmp_max) 700 701 DO i = 1, ncoa 702 p1 = (i - 1)*ncob 703 DO j = 1, ncob 704 p2 = (p1 + j - 1)*ncod 705 DO l = 1, ncod 706 p3 = (p2 + l - 1)*ncoc 707 DO k = 1, ncoc 708 p4 = p3 + k 709 work2(i, j, k, l, full_perm3(n)) = work(p4, n) 710 END DO 711 END DO 712 END DO 713 END DO 714 END DO 715 IF (use_virial) THEN 716 DO n = 1, 12 717 tmp_max_virial = 0.0_dp 718 DO i = 1, mysize 719 tmp_max_virial = MAX(tmp_max_virial, & 720 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 721 END DO 722 tmp_max_virial = tmp_max_virial*max_contraction 723 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 724 725 DO i = 1, ncoa 726 p1 = (i - 1)*ncob 727 DO j = 1, ncob 728 p2 = (p1 + j - 1)*ncod 729 DO l = 1, ncod 730 p3 = (p2 + l - 1)*ncoc 731 DO k = 1, ncoc 732 p4 = p3 + k 733 work2_virial(i, j, k, l, full_perm3(n), 1:3) = work_virial(p4, n, 1:3) 734 END DO 735 END DO 736 END DO 737 END DO 738 END DO 739 END IF 740 CASE (4) 741 DO i = 1, nproducts 742 A = pgf_product_list(i)%ra 743 B = pgf_product_list(i)%rb 744 C = pgf_product_list(i)%rc 745 D = pgf_product_list(i)%rd 746 Rho = pgf_product_list(i)%Rho 747 RhoInv = pgf_product_list(i)%RhoInv 748 P = pgf_product_list(i)%P 749 Q = pgf_product_list(i)%Q 750 W = pgf_product_list(i)%W 751 AB = pgf_product_list(i)%AB 752 CD = pgf_product_list(i)%CD 753 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 754 CALL build_deriv_data(deriv, B, A, D, C, & 755 Zeta_B, Zeta_A, Zeta_D, Zeta_C, & 756 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, & 757 P, Q, W, m_max, pgf_product_list(i)%Fm) 758 CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize) 759 DO k = 4, 6 760 DO j = 1, mysize 761 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 762 work_forces(j, k + 3) + & 763 work_forces(j, k + 6)) 764 END DO 765 END DO 766 DO k = 1, 12 767 DO j = 1, mysize 768 work(j, k) = work(j, k) + work_forces(j, k) 769 END DO 770 END DO 771 neris = neris + 12*mysize 772 IF (use_virial) THEN 773 CALL real_to_scaled(scoord(1:3), B, cell) 774 CALL real_to_scaled(scoord(4:6), A, cell) 775 CALL real_to_scaled(scoord(7:9), D, cell) 776 CALL real_to_scaled(scoord(10:12), C, cell) 777 DO k = 1, 12 778 DO j = 1, mysize 779 DO m = 1, 3 780 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 781 END DO 782 END DO 783 END DO 784 END IF 785 786 END DO 787 DO n = 1, 12 788 tmp_max = 0.0_dp 789 DO i = 1, mysize 790 tmp_max = MAX(tmp_max, ABS(work(i, n))) 791 END DO 792 tmp_max = tmp_max*max_contraction 793 tmp_max_all = MAX(tmp_max_all, tmp_max) 794 795 DO j = 1, ncob 796 p1 = (j - 1)*ncoa 797 DO i = 1, ncoa 798 p2 = (p1 + i - 1)*ncod 799 DO l = 1, ncod 800 p3 = (p2 + l - 1)*ncoc 801 DO k = 1, ncoc 802 p4 = p3 + k 803 work2(i, j, k, l, full_perm4(n)) = work(p4, n) 804 END DO 805 END DO 806 END DO 807 END DO 808 END DO 809 IF (use_virial) THEN 810 DO n = 1, 12 811 tmp_max_virial = 0.0_dp 812 DO i = 1, mysize 813 tmp_max_virial = MAX(tmp_max_virial, & 814 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 815 END DO 816 tmp_max_virial = tmp_max_virial*max_contraction 817 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 818 819 DO j = 1, ncob 820 p1 = (j - 1)*ncoa 821 DO i = 1, ncoa 822 p2 = (p1 + i - 1)*ncod 823 DO l = 1, ncod 824 p3 = (p2 + l - 1)*ncoc 825 DO k = 1, ncoc 826 p4 = p3 + k 827 work2_virial(i, j, k, l, full_perm4(n), 1:3) = work_virial(p4, n, 1:3) 828 END DO 829 END DO 830 END DO 831 END DO 832 END DO 833 END IF 834 CASE (5) 835 DO i = 1, nproducts 836 A = pgf_product_list(i)%ra 837 B = pgf_product_list(i)%rb 838 C = pgf_product_list(i)%rc 839 D = pgf_product_list(i)%rd 840 Rho = pgf_product_list(i)%Rho 841 RhoInv = pgf_product_list(i)%RhoInv 842 P = pgf_product_list(i)%P 843 Q = pgf_product_list(i)%Q 844 W = pgf_product_list(i)%W 845 AB = pgf_product_list(i)%AB 846 CD = pgf_product_list(i)%CD 847 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 848 CALL build_deriv_data(deriv, C, D, A, B, & 849 Zeta_C, Zeta_D, Zeta_A, Zeta_B, & 850 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, & 851 Q, P, W, m_max, pgf_product_list(i)%Fm) 852 CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize) 853 DO k = 4, 6 854 DO j = 1, mysize 855 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 856 work_forces(j, k + 3) + & 857 work_forces(j, k + 6)) 858 END DO 859 END DO 860 DO k = 1, 12 861 DO j = 1, mysize 862 work(j, k) = work(j, k) + work_forces(j, k) 863 END DO 864 END DO 865 neris = neris + 12*mysize 866 IF (use_virial) THEN 867 CALL real_to_scaled(scoord(1:3), C, cell) 868 CALL real_to_scaled(scoord(4:6), D, cell) 869 CALL real_to_scaled(scoord(7:9), A, cell) 870 CALL real_to_scaled(scoord(10:12), B, cell) 871 DO k = 1, 12 872 DO j = 1, mysize 873 DO m = 1, 3 874 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 875 END DO 876 END DO 877 END DO 878 END IF 879 880 END DO 881 DO n = 1, 12 882 tmp_max = 0.0_dp 883 DO i = 1, mysize 884 tmp_max = MAX(tmp_max, ABS(work(i, n))) 885 END DO 886 tmp_max = tmp_max*max_contraction 887 tmp_max_all = MAX(tmp_max_all, tmp_max) 888 889 DO k = 1, ncoc 890 p1 = (k - 1)*ncod 891 DO l = 1, ncod 892 p2 = (p1 + l - 1)*ncoa 893 DO i = 1, ncoa 894 p3 = (p2 + i - 1)*ncob 895 DO j = 1, ncob 896 p4 = p3 + j 897 work2(i, j, k, l, full_perm5(n)) = work(p4, n) 898 END DO 899 END DO 900 END DO 901 END DO 902 END DO 903 IF (use_virial) THEN 904 DO n = 1, 12 905 tmp_max_virial = 0.0_dp 906 DO i = 1, mysize 907 tmp_max_virial = MAX(tmp_max_virial, & 908 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 909 END DO 910 tmp_max_virial = tmp_max_virial*max_contraction 911 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 912 913 DO k = 1, ncoc 914 p1 = (k - 1)*ncod 915 DO l = 1, ncod 916 p2 = (p1 + l - 1)*ncoa 917 DO i = 1, ncoa 918 p3 = (p2 + i - 1)*ncob 919 DO j = 1, ncob 920 p4 = p3 + j 921 work2_virial(i, j, k, l, full_perm5(n), 1:3) = work_virial(p4, n, 1:3) 922 END DO 923 END DO 924 END DO 925 END DO 926 END DO 927 END IF 928 CASE (6) 929 DO i = 1, nproducts 930 A = pgf_product_list(i)%ra 931 B = pgf_product_list(i)%rb 932 C = pgf_product_list(i)%rc 933 D = pgf_product_list(i)%rd 934 Rho = pgf_product_list(i)%Rho 935 RhoInv = pgf_product_list(i)%RhoInv 936 P = pgf_product_list(i)%P 937 Q = pgf_product_list(i)%Q 938 W = pgf_product_list(i)%W 939 AB = pgf_product_list(i)%AB 940 CD = pgf_product_list(i)%CD 941 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 942 943 CALL build_deriv_data(deriv, C, D, B, A, & 944 Zeta_C, Zeta_D, Zeta_B, Zeta_A, & 945 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, & 946 Q, P, W, m_max, pgf_product_list(i)%Fm) 947 CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize) 948 DO k = 4, 6 949 DO j = 1, mysize 950 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 951 work_forces(j, k + 3) + & 952 work_forces(j, k + 6)) 953 END DO 954 END DO 955 DO k = 1, 12 956 DO j = 1, mysize 957 work(j, k) = work(j, k) + work_forces(j, k) 958 END DO 959 END DO 960 neris = neris + 12*mysize 961 IF (use_virial) THEN 962 CALL real_to_scaled(scoord(1:3), C, cell) 963 CALL real_to_scaled(scoord(4:6), D, cell) 964 CALL real_to_scaled(scoord(7:9), B, cell) 965 CALL real_to_scaled(scoord(10:12), A, cell) 966 DO k = 1, 12 967 DO j = 1, mysize 968 DO m = 1, 3 969 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 970 END DO 971 END DO 972 END DO 973 END IF 974 975 END DO 976 DO n = 1, 12 977 tmp_max = 0.0_dp 978 DO i = 1, mysize 979 tmp_max = MAX(tmp_max, ABS(work(i, n))) 980 END DO 981 tmp_max = tmp_max*max_contraction 982 tmp_max_all = MAX(tmp_max_all, tmp_max) 983 984 DO k = 1, ncoc 985 p1 = (k - 1)*ncod 986 DO l = 1, ncod 987 p2 = (p1 + l - 1)*ncob 988 DO j = 1, ncob 989 p3 = (p2 + j - 1)*ncoa 990 DO i = 1, ncoa 991 p4 = p3 + i 992 work2(i, j, k, l, full_perm6(n)) = work(p4, n) 993 END DO 994 END DO 995 END DO 996 END DO 997 END DO 998 IF (use_virial) THEN 999 DO n = 1, 12 1000 tmp_max_virial = 0.0_dp 1001 DO i = 1, mysize 1002 tmp_max_virial = MAX(tmp_max_virial, & 1003 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 1004 END DO 1005 tmp_max_virial = tmp_max_virial*max_contraction 1006 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 1007 1008 DO k = 1, ncoc 1009 p1 = (k - 1)*ncod 1010 DO l = 1, ncod 1011 p2 = (p1 + l - 1)*ncob 1012 DO j = 1, ncob 1013 p3 = (p2 + j - 1)*ncoa 1014 DO i = 1, ncoa 1015 p4 = p3 + i 1016 work2_virial(i, j, k, l, full_perm6(n), 1:3) = work_virial(p4, n, 1:3) 1017 END DO 1018 END DO 1019 END DO 1020 END DO 1021 END DO 1022 END IF 1023 CASE (7) 1024 DO i = 1, nproducts 1025 A = pgf_product_list(i)%ra 1026 B = pgf_product_list(i)%rb 1027 C = pgf_product_list(i)%rc 1028 D = pgf_product_list(i)%rd 1029 Rho = pgf_product_list(i)%Rho 1030 RhoInv = pgf_product_list(i)%RhoInv 1031 P = pgf_product_list(i)%P 1032 Q = pgf_product_list(i)%Q 1033 W = pgf_product_list(i)%W 1034 AB = pgf_product_list(i)%AB 1035 CD = pgf_product_list(i)%CD 1036 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1037 1038 CALL build_deriv_data(deriv, D, C, A, B, & 1039 Zeta_D, Zeta_C, Zeta_A, Zeta_B, & 1040 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, & 1041 Q, P, W, m_max, pgf_product_list(i)%Fm) 1042 CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize) 1043 DO k = 4, 6 1044 DO j = 1, mysize 1045 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 1046 work_forces(j, k + 3) + & 1047 work_forces(j, k + 6)) 1048 END DO 1049 END DO 1050 DO k = 1, 12 1051 DO j = 1, mysize 1052 work(j, k) = work(j, k) + work_forces(j, k) 1053 END DO 1054 END DO 1055 neris = neris + 12*mysize 1056 IF (use_virial) THEN 1057 CALL real_to_scaled(scoord(1:3), D, cell) 1058 CALL real_to_scaled(scoord(4:6), C, cell) 1059 CALL real_to_scaled(scoord(7:9), A, cell) 1060 CALL real_to_scaled(scoord(10:12), B, cell) 1061 DO k = 1, 12 1062 DO j = 1, mysize 1063 DO m = 1, 3 1064 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 1065 END DO 1066 END DO 1067 END DO 1068 END IF 1069 1070 END DO 1071 DO n = 1, 12 1072 tmp_max = 0.0_dp 1073 DO i = 1, mysize 1074 tmp_max = MAX(tmp_max, ABS(work(i, n))) 1075 END DO 1076 tmp_max = tmp_max*max_contraction 1077 tmp_max_all = MAX(tmp_max_all, tmp_max) 1078 1079 DO l = 1, ncod 1080 p1 = (l - 1)*ncoc 1081 DO k = 1, ncoc 1082 p2 = (p1 + k - 1)*ncoa 1083 DO i = 1, ncoa 1084 p3 = (p2 + i - 1)*ncob 1085 DO j = 1, ncob 1086 p4 = p3 + j 1087 work2(i, j, k, l, full_perm7(n)) = work(p4, n) 1088 END DO 1089 END DO 1090 END DO 1091 END DO 1092 END DO 1093 IF (use_virial) THEN 1094 DO n = 1, 12 1095 tmp_max_virial = 0.0_dp 1096 DO i = 1, mysize 1097 tmp_max_virial = MAX(tmp_max_virial, & 1098 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 1099 END DO 1100 tmp_max_virial = tmp_max_virial*max_contraction 1101 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 1102 1103 DO l = 1, ncod 1104 p1 = (l - 1)*ncoc 1105 DO k = 1, ncoc 1106 p2 = (p1 + k - 1)*ncoa 1107 DO i = 1, ncoa 1108 p3 = (p2 + i - 1)*ncob 1109 DO j = 1, ncob 1110 p4 = p3 + j 1111 work2_virial(i, j, k, l, full_perm7(n), 1:3) = work_virial(p4, n, 1:3) 1112 END DO 1113 END DO 1114 END DO 1115 END DO 1116 END DO 1117 END IF 1118 CASE (8) 1119 DO i = 1, nproducts 1120 A = pgf_product_list(i)%ra 1121 B = pgf_product_list(i)%rb 1122 C = pgf_product_list(i)%rc 1123 D = pgf_product_list(i)%rd 1124 Rho = pgf_product_list(i)%Rho 1125 RhoInv = pgf_product_list(i)%RhoInv 1126 P = pgf_product_list(i)%P 1127 Q = pgf_product_list(i)%Q 1128 W = pgf_product_list(i)%W 1129 AB = pgf_product_list(i)%AB 1130 CD = pgf_product_list(i)%CD 1131 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1132 1133 CALL build_deriv_data(deriv, D, C, B, A, & 1134 Zeta_D, Zeta_C, Zeta_B, Zeta_A, & 1135 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv, & 1136 Q, P, W, m_max, pgf_product_list(i)%Fm) 1137 CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize) 1138 DO k = 4, 6 1139 DO j = 1, mysize 1140 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + & 1141 work_forces(j, k + 3) + & 1142 work_forces(j, k + 6)) 1143 END DO 1144 END DO 1145 DO k = 1, 12 1146 DO j = 1, mysize 1147 work(j, k) = work(j, k) + work_forces(j, k) 1148 END DO 1149 END DO 1150 neris = neris + 12*mysize 1151 IF (use_virial) THEN 1152 CALL real_to_scaled(scoord(1:3), D, cell) 1153 CALL real_to_scaled(scoord(4:6), C, cell) 1154 CALL real_to_scaled(scoord(7:9), B, cell) 1155 CALL real_to_scaled(scoord(10:12), A, cell) 1156 DO k = 1, 12 1157 DO j = 1, mysize 1158 DO m = 1, 3 1159 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(INT((k - 1)/3)*3 + m) 1160 END DO 1161 END DO 1162 END DO 1163 END IF 1164 1165 END DO 1166 DO n = 1, 12 1167 tmp_max = 0.0_dp 1168 DO i = 1, mysize 1169 tmp_max = MAX(tmp_max, ABS(work(i, n))) 1170 END DO 1171 tmp_max = tmp_max*max_contraction 1172 tmp_max_all = MAX(tmp_max_all, tmp_max) 1173 1174 DO l = 1, ncod 1175 p1 = (l - 1)*ncoc 1176 DO k = 1, ncoc 1177 p2 = (p1 + k - 1)*ncob 1178 DO j = 1, ncob 1179 p3 = (p2 + j - 1)*ncoa 1180 DO i = 1, ncoa 1181 p4 = p3 + i 1182 work2(i, j, k, l, full_perm8(n)) = work(p4, n) 1183 END DO 1184 END DO 1185 END DO 1186 END DO 1187 END DO 1188 IF (use_virial) THEN 1189 DO n = 1, 12 1190 tmp_max_virial = 0.0_dp 1191 DO i = 1, mysize 1192 tmp_max_virial = MAX(tmp_max_virial, & 1193 ABS(work_virial(i, n, 1)), ABS(work_virial(i, n, 2)), ABS(work_virial(i, n, 3))) 1194 END DO 1195 tmp_max_virial = tmp_max_virial*max_contraction 1196 tmp_max_all_virial = MAX(tmp_max_all_virial, tmp_max_virial) 1197 1198 DO l = 1, ncod 1199 p1 = (l - 1)*ncoc 1200 DO k = 1, ncoc 1201 p2 = (p1 + k - 1)*ncob 1202 DO j = 1, ncob 1203 p3 = (p2 + j - 1)*ncoa 1204 DO i = 1, ncoa 1205 p4 = p3 + i 1206 work2_virial(i, j, k, l, full_perm8(n), 1:3) = work_virial(p4, n, 1:3) 1207 END DO 1208 END DO 1209 END DO 1210 END DO 1211 END DO 1212 END IF 1213 END SELECT 1214 1215 IF (.NOT. use_virial) THEN 1216 IF (tmp_max_all < eps_schwarz) RETURN 1217 END IF 1218 1219 IF (tmp_max_all >= eps_schwarz) THEN 1220 DO i = 1, 12 1221 primitives_tmp(:, :, :, :) = 0.0_dp 1222 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, & 1223 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2(1, 1, 1, 1, i), & 1224 sphi_a, & 1225 sphi_b, & 1226 sphi_c, & 1227 sphi_d, & 1228 primitives_tmp(1, 1, 1, 1), & 1229 buffer1, buffer2) 1230 1231 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, & 1232 s_offset_b + 1:s_offset_b + nsob*nl_b, & 1233 s_offset_c + 1:s_offset_c + nsoc*nl_c, & 1234 s_offset_d + 1:s_offset_d + nsod*nl_d, i) = & 1235 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, & 1236 s_offset_b + 1:s_offset_b + nsob*nl_b, & 1237 s_offset_c + 1:s_offset_c + nsoc*nl_c, & 1238 s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :) 1239 END DO 1240 END IF 1241 1242 IF (use_virial .AND. tmp_max_all_virial >= eps_schwarz) THEN 1243 DO i = 1, 12 1244 DO m = 1, 3 1245 primitives_tmp_virial(:, :, :, :) = 0.0_dp 1246 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, & 1247 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2_virial(1, 1, 1, 1, i, m), & 1248 sphi_a, & 1249 sphi_b, & 1250 sphi_c, & 1251 sphi_d, & 1252 primitives_tmp_virial(1, 1, 1, 1), & 1253 buffer1, buffer2) 1254 1255 primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, & 1256 s_offset_b + 1:s_offset_b + nsob*nl_b, & 1257 s_offset_c + 1:s_offset_c + nsoc*nl_c, & 1258 s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) = & 1259 primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, & 1260 s_offset_b + 1:s_offset_b + nsob*nl_b, & 1261 s_offset_c + 1:s_offset_c + nsoc*nl_c, & 1262 s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) + primitives_tmp_virial(:, :, :, :) 1263 END DO 1264 END DO 1265 END IF 1266 1267 END SUBROUTINE evaluate_deriv_eri 1268 1269! ************************************************************************************************** 1270!> \brief Evaluates max-abs values of electron repulsion integrals for a primitive quartet 1271!> \param libint ... 1272!> \param A ... 1273!> \param B ... 1274!> \param C ... 1275!> \param D ... 1276!> \param Zeta_A ... 1277!> \param Zeta_B ... 1278!> \param Zeta_C ... 1279!> \param Zeta_D ... 1280!> \param n_a ... 1281!> \param n_b ... 1282!> \param n_c ... 1283!> \param n_d ... 1284!> \param max_val ... 1285!> \param potential_parameter ... 1286!> \param R1 ... 1287!> \param R2 ... 1288!> \param p_work ... 1289!> \par History 1290!> 03.2007 created [Manuel Guidon] 1291!> 08.2007 refactured permutation part [Manuel Guidon] 1292!> \author Manuel Guidon 1293! ************************************************************************************************** 1294 SUBROUTINE evaluate_eri_screen(libint, A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, & 1295 n_a, n_b, n_c, n_d, & 1296 max_val, potential_parameter, R1, R2, & 1297 p_work) 1298 1299 TYPE(cp_libint_t) :: libint 1300 REAL(dp), INTENT(IN) :: A(3), B(3), C(3), D(3), Zeta_A, Zeta_B, & 1301 Zeta_C, Zeta_D 1302 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d 1303 REAL(dp), INTENT(INOUT) :: max_val 1304 TYPE(hfx_potential_type) :: potential_parameter 1305 REAL(dp) :: R1, R2 1306 REAL(dp), DIMENSION(:), POINTER :: p_work 1307 1308 INTEGER :: a_mysize(1), i, m_max, mysize, perm_case 1309 1310 m_max = n_a + n_b + n_c + n_d 1311 mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d) 1312 a_mysize = mysize 1313 1314 IF (m_max /= 0) THEN 1315 perm_case = 1 1316 IF (n_a < n_b) THEN 1317 perm_case = perm_case + 1 1318 END IF 1319 IF (n_c < n_d) THEN 1320 perm_case = perm_case + 2 1321 END IF 1322 IF (n_a + n_b > n_c + n_d) THEN 1323 perm_case = perm_case + 4 1324 END IF 1325 1326 SELECT CASE (perm_case) 1327 CASE (1) 1328 CALL build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, & 1329 potential_parameter, libint, R1, R2) 1330 1331 CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize) 1332 DO i = 1, mysize 1333 max_val = MAX(max_val, ABS(p_work(i))) 1334 END DO 1335 CASE (2) 1336 CALL build_quartet_data_screen(B, A, C, D, Zeta_B, Zeta_A, Zeta_C, Zeta_D, m_max, & 1337 potential_parameter, libint, R1, R2) 1338 1339 CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize) 1340 DO i = 1, mysize 1341 max_val = MAX(max_val, ABS(p_work(i))) 1342 END DO 1343 CASE (3) 1344 CALL build_quartet_data_screen(A, B, D, C, Zeta_A, Zeta_B, Zeta_D, Zeta_C, m_max, & 1345 potential_parameter, libint, R1, R2) 1346 1347 CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize) 1348 DO i = 1, mysize 1349 max_val = MAX(max_val, ABS(p_work(i))) 1350 END DO 1351 CASE (4) 1352 CALL build_quartet_data_screen(B, A, D, C, Zeta_B, Zeta_A, Zeta_D, Zeta_C, m_max, & 1353 potential_parameter, libint, R1, R2) 1354 1355 CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize) 1356 DO i = 1, mysize 1357 max_val = MAX(max_val, ABS(p_work(i))) 1358 END DO 1359 CASE (5) 1360 CALL build_quartet_data_screen(C, D, A, B, Zeta_C, Zeta_D, Zeta_A, Zeta_B, m_max, & 1361 potential_parameter, libint, R1, R2) 1362 1363 CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize) 1364 DO i = 1, mysize 1365 max_val = MAX(max_val, ABS(p_work(i))) 1366 END DO 1367 CASE (6) 1368 CALL build_quartet_data_screen(C, D, B, A, Zeta_C, Zeta_D, Zeta_B, Zeta_A, m_max, & 1369 potential_parameter, libint, R1, R2) 1370 1371 CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize) 1372 DO i = 1, mysize 1373 max_val = MAX(max_val, ABS(p_work(i))) 1374 END DO 1375 CASE (7) 1376 CALL build_quartet_data_screen(D, C, A, B, Zeta_D, Zeta_C, Zeta_A, Zeta_B, m_max, & 1377 potential_parameter, libint, R1, R2) 1378 1379 CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize) 1380 DO i = 1, mysize 1381 max_val = MAX(max_val, ABS(p_work(i))) 1382 END DO 1383 CASE (8) 1384 CALL build_quartet_data_screen(D, C, B, A, Zeta_D, Zeta_C, Zeta_B, Zeta_A, m_max, & 1385 potential_parameter, libint, R1, R2) 1386 1387 CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize) 1388 DO i = 1, mysize 1389 max_val = MAX(max_val, ABS(p_work(i))) 1390 END DO 1391 END SELECT 1392 ELSE 1393 CALL build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, & 1394 potential_parameter, libint, R1, R2) 1395 max_val = ABS(get_ssss_f_val(libint)) 1396 END IF 1397 1398 END SUBROUTINE evaluate_eri_screen 1399 1400! ************************************************************************************************** 1401!> \brief Evaluate electron repulsion integrals for a primitive quartet 1402!> \param libint ... 1403!> \param nproducts ... 1404!> \param pgf_product_list ... 1405!> \param n_a ... 1406!> \param n_b ... 1407!> \param n_c ... 1408!> \param n_d ... 1409!> \param ncoa ... 1410!> \param ncob ... 1411!> \param ncoc ... 1412!> \param ncod ... 1413!> \param nsgfa ... 1414!> \param nsgfb ... 1415!> \param nsgfc ... 1416!> \param nsgfd ... 1417!> \param primitives ... 1418!> \param max_contraction ... 1419!> \param tmp_max ... 1420!> \param eps_schwarz ... 1421!> \param neris ... 1422!> \param ZetaInv ... 1423!> \param EtaInv ... 1424!> \param s_offset_a ... 1425!> \param s_offset_b ... 1426!> \param s_offset_c ... 1427!> \param s_offset_d ... 1428!> \param nl_a ... 1429!> \param nl_b ... 1430!> \param nl_c ... 1431!> \param nl_d ... 1432!> \param nsoa ... 1433!> \param nsob ... 1434!> \param nsoc ... 1435!> \param nsod ... 1436!> \param sphi_a ... 1437!> \param sphi_b ... 1438!> \param sphi_c ... 1439!> \param sphi_d ... 1440!> \param work ... 1441!> \param work2 ... 1442!> \param buffer1 ... 1443!> \param buffer2 ... 1444!> \param primitives_tmp ... 1445!> \param p_work ... 1446!> \par History 1447!> 11.2006 created [Manuel Guidon] 1448!> 08.2007 refactured permutation part [Manuel Guidon] 1449!> \author Manuel Guidon 1450! ************************************************************************************************** 1451 SUBROUTINE evaluate_eri(libint, nproducts, pgf_product_list, & 1452 n_a, n_b, n_c, n_d, & 1453 ncoa, ncob, ncoc, ncod, & 1454 nsgfa, nsgfb, nsgfc, nsgfd, & 1455 primitives, max_contraction, tmp_max, & 1456 eps_schwarz, neris, & 1457 ZetaInv, EtaInv, & 1458 s_offset_a, s_offset_b, s_offset_c, s_offset_d, & 1459 nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, & 1460 sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, & 1461 primitives_tmp, p_work) 1462 1463 TYPE(cp_libint_t) :: libint 1464 INTEGER, INTENT(IN) :: nproducts 1465 TYPE(hfx_pgf_product_list), DIMENSION(nproducts) :: pgf_product_list 1466 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, & 1467 ncod, nsgfa, nsgfb, nsgfc, nsgfd 1468 REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitives 1469 REAL(dp) :: max_contraction, tmp_max, eps_schwarz 1470 INTEGER(int_8) :: neris 1471 REAL(dp), INTENT(IN) :: ZetaInv, EtaInv 1472 INTEGER :: s_offset_a, s_offset_b, s_offset_c, & 1473 s_offset_d, nl_a, nl_b, nl_c, nl_d, & 1474 nsoa, nsob, nsoc, nsod 1475 REAL(dp), DIMENSION(ncoa, nsoa*nl_a), INTENT(IN) :: sphi_a 1476 REAL(dp), DIMENSION(ncob, nsob*nl_b), INTENT(IN) :: sphi_b 1477 REAL(dp), DIMENSION(ncoc, nsoc*nl_c), INTENT(IN) :: sphi_c 1478 REAL(dp), DIMENSION(ncod, nsod*nl_d), INTENT(IN) :: sphi_d 1479 REAL(dp) :: work(ncoa*ncob*ncoc*ncod), & 1480 work2(ncoa, ncob, ncoc, ncod) 1481 REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2 1482 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*& 1483 nl_c, nsod*nl_d) :: primitives_tmp 1484 REAL(dp), DIMENSION(:), POINTER :: p_work 1485 1486 INTEGER :: a_mysize(1), i, j, k, l, m_max, mysize, & 1487 p1, p2, p3, p4, perm_case 1488 REAL(dp) :: A(3), AB(3), B(3), C(3), CD(3), D(3), & 1489 P(3), Q(3), Rho, RhoInv, W(3), & 1490 ZetapEtaInv 1491 REAL(KIND=dp), DIMENSION(prim_data_f_size) :: F 1492 1493 m_max = n_a + n_b + n_c + n_d 1494 mysize = ncoa*ncob*ncoc*ncod 1495 a_mysize = mysize 1496 1497 work = 0.0_dp 1498 IF (m_max /= 0) THEN 1499 perm_case = 1 1500 IF (n_a < n_b) THEN 1501 perm_case = perm_case + 1 1502 END IF 1503 IF (n_c < n_d) THEN 1504 perm_case = perm_case + 2 1505 END IF 1506 IF (n_a + n_b > n_c + n_d) THEN 1507 perm_case = perm_case + 4 1508 END IF 1509 SELECT CASE (perm_case) 1510 CASE (1) 1511 DO i = 1, nproducts 1512 A = pgf_product_list(i)%ra 1513 B = pgf_product_list(i)%rb 1514 C = pgf_product_list(i)%rc 1515 D = pgf_product_list(i)%rd 1516 Rho = pgf_product_list(i)%Rho 1517 RhoInv = pgf_product_list(i)%RhoInv 1518 P = pgf_product_list(i)%P 1519 Q = pgf_product_list(i)%Q 1520 W = pgf_product_list(i)%W 1521 AB = pgf_product_list(i)%AB 1522 CD = pgf_product_list(i)%CD 1523 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1524 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1525 1526 CALL build_quartet_data(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, & 1527 P, Q, W, m_max, F) 1528 CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize) 1529 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1530 neris = neris + mysize 1531 END DO 1532 DO i = 1, mysize 1533 tmp_max = MAX(tmp_max, ABS(work(i))) 1534 END DO 1535 tmp_max = tmp_max*max_contraction 1536 IF (tmp_max < eps_schwarz) THEN 1537 RETURN 1538 END IF 1539 1540 DO i = 1, ncoa 1541 p1 = (i - 1)*ncob 1542 DO j = 1, ncob 1543 p2 = (p1 + j - 1)*ncoc 1544 DO k = 1, ncoc 1545 p3 = (p2 + k - 1)*ncod 1546 DO l = 1, ncod 1547 p4 = p3 + l 1548 work2(i, j, k, l) = work(p4) 1549 END DO 1550 END DO 1551 END DO 1552 END DO 1553 CASE (2) 1554 DO i = 1, nproducts 1555 A = pgf_product_list(i)%ra 1556 B = pgf_product_list(i)%rb 1557 C = pgf_product_list(i)%rc 1558 D = pgf_product_list(i)%rd 1559 Rho = pgf_product_list(i)%Rho 1560 RhoInv = pgf_product_list(i)%RhoInv 1561 P = pgf_product_list(i)%P 1562 Q = pgf_product_list(i)%Q 1563 W = pgf_product_list(i)%W 1564 AB = pgf_product_list(i)%AB 1565 CD = pgf_product_list(i)%CD 1566 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1567 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1568 1569 CALL build_quartet_data(libint, B, A, C, D, & 1570 ZetaInv, EtaInv, ZetapEtaInv, Rho, & 1571 P, Q, W, m_max, F) 1572 1573 CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize) 1574 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1575 neris = neris + mysize 1576 END DO 1577 DO i = 1, mysize 1578 tmp_max = MAX(tmp_max, ABS(work(i))) 1579 END DO 1580 tmp_max = tmp_max*max_contraction 1581 IF (tmp_max < eps_schwarz) THEN 1582 RETURN 1583 END IF 1584 1585 DO j = 1, ncob 1586 p1 = (j - 1)*ncoa 1587 DO i = 1, ncoa 1588 p2 = (p1 + i - 1)*ncoc 1589 DO k = 1, ncoc 1590 p3 = (p2 + k - 1)*ncod 1591 DO l = 1, ncod 1592 p4 = p3 + l 1593 work2(i, j, k, l) = work(p4) 1594 END DO 1595 END DO 1596 END DO 1597 END DO 1598 CASE (3) 1599 DO i = 1, nproducts 1600 A = pgf_product_list(i)%ra 1601 B = pgf_product_list(i)%rb 1602 C = pgf_product_list(i)%rc 1603 D = pgf_product_list(i)%rd 1604 Rho = pgf_product_list(i)%Rho 1605 RhoInv = pgf_product_list(i)%RhoInv 1606 P = pgf_product_list(i)%P 1607 Q = pgf_product_list(i)%Q 1608 W = pgf_product_list(i)%W 1609 AB = pgf_product_list(i)%AB 1610 CD = pgf_product_list(i)%CD 1611 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1612 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1613 1614 CALL build_quartet_data(libint, A, B, D, C, & 1615 ZetaInv, EtaInv, ZetapEtaInv, Rho, & 1616 P, Q, W, m_max, F) 1617 1618 CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize) 1619 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1620 neris = neris + mysize 1621 END DO 1622 DO i = 1, mysize 1623 tmp_max = MAX(tmp_max, ABS(work(i))) 1624 END DO 1625 tmp_max = tmp_max*max_contraction 1626 IF (tmp_max < eps_schwarz) THEN 1627 RETURN 1628 END IF 1629 1630 DO i = 1, ncoa 1631 p1 = (i - 1)*ncob 1632 DO j = 1, ncob 1633 p2 = (p1 + j - 1)*ncod 1634 DO l = 1, ncod 1635 p3 = (p2 + l - 1)*ncoc 1636 DO k = 1, ncoc 1637 p4 = p3 + k 1638 work2(i, j, k, l) = work(p4) 1639 END DO 1640 END DO 1641 END DO 1642 END DO 1643 CASE (4) 1644 DO i = 1, nproducts 1645 A = pgf_product_list(i)%ra 1646 B = pgf_product_list(i)%rb 1647 C = pgf_product_list(i)%rc 1648 D = pgf_product_list(i)%rd 1649 Rho = pgf_product_list(i)%Rho 1650 RhoInv = pgf_product_list(i)%RhoInv 1651 P = pgf_product_list(i)%P 1652 Q = pgf_product_list(i)%Q 1653 W = pgf_product_list(i)%W 1654 AB = pgf_product_list(i)%AB 1655 CD = pgf_product_list(i)%CD 1656 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1657 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1658 1659 CALL build_quartet_data(libint, B, A, D, C, & 1660 ZetaInv, EtaInv, ZetapEtaInv, Rho, & 1661 P, Q, W, m_max, F) 1662 1663 CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize) 1664 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1665 neris = neris + mysize 1666 END DO 1667 DO i = 1, mysize 1668 tmp_max = MAX(tmp_max, ABS(work(i))) 1669 END DO 1670 tmp_max = tmp_max*max_contraction 1671 IF (tmp_max < eps_schwarz) THEN 1672 RETURN 1673 END IF 1674 1675 DO j = 1, ncob 1676 p1 = (j - 1)*ncoa 1677 DO i = 1, ncoa 1678 p2 = (p1 + i - 1)*ncod 1679 DO l = 1, ncod 1680 p3 = (p2 + l - 1)*ncoc 1681 DO k = 1, ncoc 1682 p4 = p3 + k 1683 work2(i, j, k, l) = work(p4) 1684 END DO 1685 END DO 1686 END DO 1687 END DO 1688 CASE (5) 1689 DO i = 1, nproducts 1690 A = pgf_product_list(i)%ra 1691 B = pgf_product_list(i)%rb 1692 C = pgf_product_list(i)%rc 1693 D = pgf_product_list(i)%rd 1694 Rho = pgf_product_list(i)%Rho 1695 RhoInv = pgf_product_list(i)%RhoInv 1696 P = pgf_product_list(i)%P 1697 Q = pgf_product_list(i)%Q 1698 W = pgf_product_list(i)%W 1699 AB = pgf_product_list(i)%AB 1700 CD = pgf_product_list(i)%CD 1701 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1702 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1703 1704 CALL build_quartet_data(libint, C, D, A, B, & 1705 EtaInv, ZetaInv, ZetapEtaInv, Rho, & 1706 Q, P, W, m_max, F) 1707 1708 CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize) 1709 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1710 neris = neris + mysize 1711 END DO 1712 DO i = 1, mysize 1713 tmp_max = MAX(tmp_max, ABS(work(i))) 1714 END DO 1715 tmp_max = tmp_max*max_contraction 1716 IF (tmp_max < eps_schwarz) THEN 1717 RETURN 1718 END IF 1719 1720 DO k = 1, ncoc 1721 p1 = (k - 1)*ncod 1722 DO l = 1, ncod 1723 p2 = (p1 + l - 1)*ncoa 1724 DO i = 1, ncoa 1725 p3 = (p2 + i - 1)*ncob 1726 DO j = 1, ncob 1727 p4 = p3 + j 1728 work2(i, j, k, l) = work(p4) 1729 END DO 1730 END DO 1731 END DO 1732 END DO 1733 CASE (6) 1734 DO i = 1, nproducts 1735 A = pgf_product_list(i)%ra 1736 B = pgf_product_list(i)%rb 1737 C = pgf_product_list(i)%rc 1738 D = pgf_product_list(i)%rd 1739 Rho = pgf_product_list(i)%Rho 1740 RhoInv = pgf_product_list(i)%RhoInv 1741 P = pgf_product_list(i)%P 1742 Q = pgf_product_list(i)%Q 1743 W = pgf_product_list(i)%W 1744 AB = pgf_product_list(i)%AB 1745 CD = pgf_product_list(i)%CD 1746 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1747 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1748 1749 CALL build_quartet_data(libint, C, D, B, A, & 1750 EtaInv, ZetaInv, ZetapEtaInv, Rho, & 1751 Q, P, W, m_max, F) 1752 1753 CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize) 1754 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1755 neris = neris + mysize 1756 END DO 1757 DO i = 1, mysize 1758 tmp_max = MAX(tmp_max, ABS(work(i))) 1759 END DO 1760 tmp_max = tmp_max*max_contraction 1761 IF (tmp_max < eps_schwarz) THEN 1762 RETURN 1763 END IF 1764 1765 DO k = 1, ncoc 1766 p1 = (k - 1)*ncod 1767 DO l = 1, ncod 1768 p2 = (p1 + l - 1)*ncob 1769 DO j = 1, ncob 1770 p3 = (p2 + j - 1)*ncoa 1771 DO i = 1, ncoa 1772 p4 = p3 + i 1773 work2(i, j, k, l) = work(p4) 1774 END DO 1775 END DO 1776 END DO 1777 END DO 1778 CASE (7) 1779 DO i = 1, nproducts 1780 A = pgf_product_list(i)%ra 1781 B = pgf_product_list(i)%rb 1782 C = pgf_product_list(i)%rc 1783 D = pgf_product_list(i)%rd 1784 Rho = pgf_product_list(i)%Rho 1785 RhoInv = pgf_product_list(i)%RhoInv 1786 P = pgf_product_list(i)%P 1787 Q = pgf_product_list(i)%Q 1788 W = pgf_product_list(i)%W 1789 AB = pgf_product_list(i)%AB 1790 CD = pgf_product_list(i)%CD 1791 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1792 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1793 1794 CALL build_quartet_data(libint, D, C, A, B, & 1795 EtaInv, ZetaInv, ZetapEtaInv, Rho, & 1796 Q, P, W, m_max, F) 1797 1798 CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize) 1799 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1800 neris = neris + mysize 1801 END DO 1802 DO i = 1, mysize 1803 tmp_max = MAX(tmp_max, ABS(work(i))) 1804 END DO 1805 tmp_max = tmp_max*max_contraction 1806 IF (tmp_max < eps_schwarz) THEN 1807 RETURN 1808 END IF 1809 1810 DO l = 1, ncod 1811 p1 = (l - 1)*ncoc 1812 DO k = 1, ncoc 1813 p2 = (p1 + k - 1)*ncoa 1814 DO i = 1, ncoa 1815 p3 = (p2 + i - 1)*ncob 1816 DO j = 1, ncob 1817 p4 = p3 + j 1818 work2(i, j, k, l) = work(p4) 1819 END DO 1820 END DO 1821 END DO 1822 END DO 1823 CASE (8) 1824 DO i = 1, nproducts 1825 A = pgf_product_list(i)%ra 1826 B = pgf_product_list(i)%rb 1827 C = pgf_product_list(i)%rc 1828 D = pgf_product_list(i)%rd 1829 Rho = pgf_product_list(i)%Rho 1830 RhoInv = pgf_product_list(i)%RhoInv 1831 P = pgf_product_list(i)%P 1832 Q = pgf_product_list(i)%Q 1833 W = pgf_product_list(i)%W 1834 AB = pgf_product_list(i)%AB 1835 CD = pgf_product_list(i)%CD 1836 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1837 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1838 1839 CALL build_quartet_data(libint, D, C, B, A, & 1840 EtaInv, ZetaInv, ZetapEtaInv, Rho, & 1841 Q, P, W, m_max, F) 1842 1843 CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize) 1844 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 1845 neris = neris + mysize 1846 END DO 1847 DO i = 1, mysize 1848 tmp_max = MAX(tmp_max, ABS(work(i))) 1849 END DO 1850 tmp_max = tmp_max*max_contraction 1851 IF (tmp_max < eps_schwarz) THEN 1852 RETURN 1853 END IF 1854 1855 DO l = 1, ncod 1856 p1 = (l - 1)*ncoc 1857 DO k = 1, ncoc 1858 p2 = (p1 + k - 1)*ncob 1859 DO j = 1, ncob 1860 p3 = (p2 + j - 1)*ncoa 1861 DO i = 1, ncoa 1862 p4 = p3 + i 1863 work2(i, j, k, l) = work(p4) 1864 END DO 1865 END DO 1866 END DO 1867 END DO 1868 END SELECT 1869 ELSE 1870 DO i = 1, nproducts 1871 A = pgf_product_list(i)%ra 1872 B = pgf_product_list(i)%rb 1873 C = pgf_product_list(i)%rc 1874 D = pgf_product_list(i)%rd 1875 Rho = pgf_product_list(i)%Rho 1876 RhoInv = pgf_product_list(i)%RhoInv 1877 P = pgf_product_list(i)%P 1878 Q = pgf_product_list(i)%Q 1879 W = pgf_product_list(i)%W 1880 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 1881 F(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1) 1882 1883 CALL build_quartet_data(libint, A, B, C, D, & !todo: check 1884 ZetaInv, EtaInv, ZetapEtaInv, Rho, & 1885 P, Q, W, m_max, F) 1886 work(1) = work(1) + F(1) 1887 neris = neris + mysize 1888 END DO 1889 work2(1, 1, 1, 1) = work(1) 1890 tmp_max = max_contraction*ABS(work(1)) 1891 IF (tmp_max < eps_schwarz) RETURN 1892 END IF 1893 1894 IF (tmp_max < eps_schwarz) RETURN 1895 primitives_tmp = 0.0_dp 1896 1897 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, & 1898 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2, & 1899 sphi_a, & 1900 sphi_b, & 1901 sphi_c, & 1902 sphi_d, & 1903 primitives_tmp, & 1904 buffer1, buffer2) 1905 1906 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, & 1907 s_offset_b + 1:s_offset_b + nsob*nl_b, & 1908 s_offset_c + 1:s_offset_c + nsoc*nl_c, & 1909 s_offset_d + 1:s_offset_d + nsod*nl_d) = & 1910 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, & 1911 s_offset_b + 1:s_offset_b + nsob*nl_b, & 1912 s_offset_c + 1:s_offset_c + nsoc*nl_c, & 1913 s_offset_d + 1:s_offset_d + nsod*nl_d) + primitives_tmp(:, :, :, :) 1914 1915 END SUBROUTINE evaluate_eri 1916 1917! ************************************************************************************************** 1918!> \brief Fill data structure used in libint 1919!> \param libint ... 1920!> \param A ... 1921!> \param B ... 1922!> \param C ... 1923!> \param D ... 1924!> \param ZetaInv ... 1925!> \param EtaInv ... 1926!> \param ZetapEtaInv ... 1927!> \param Rho ... 1928!> \param P ... 1929!> \param Q ... 1930!> \param W ... 1931!> \param m_max ... 1932!> \param F ... 1933!> \par History 1934!> 03.2007 created [Manuel Guidon] 1935!> \author Manuel Guidon 1936! ************************************************************************************************** 1937 SUBROUTINE build_quartet_data(libint, A, B, C, D, & 1938 ZetaInv, EtaInv, ZetapEtaInv, Rho, & 1939 P, Q, W, m_max, F) 1940 1941 TYPE(cp_libint_t) :: libint 1942 REAL(KIND=dp), INTENT(IN) :: A(3), B(3), C(3), D(3) 1943 REAL(dp), INTENT(IN) :: ZetaInv, EtaInv, ZetapEtaInv, Rho, P(3), & 1944 Q(3), W(3) 1945 INTEGER, INTENT(IN) :: m_max 1946 REAL(KIND=dp), DIMENSION(:) :: F 1947 1948 CALL cp_libint_set_params_eri(libint, A, B, C, D, ZetaInv, EtaInv, ZetapEtaInv, Rho, P, Q, W, m_max, F) 1949 END SUBROUTINE build_quartet_data 1950 1951END MODULE hfx_libint_interface 1952