1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for calculating the s-integrals and their scalar derivatives with respect to rab2 8!> over solid harmonic Gaussian (SHG) functions + contraction routines for these integrals 9!> i) (s|O(r12)|s) where O(r12) is the overlap, coulomb operator etc. 10!> ii) (aba) and (abb) s-overlaps 11!> \par Literature (partly) 12!> T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008) 13!> T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure 14!> Theory, Wiley 15!> R. Ahlrichs, PCCP, 8, 3072 (2006) 16!> \par History 17!> created [05.2016] 18!> \author Dorothea Golze 19! ************************************************************************************************** 20MODULE s_contract_shg 21 USE gamma, ONLY: fgamma => fgamma_0 22 USE kinds, ONLY: dp 23 USE mathconstants, ONLY: dfac,& 24 fac,& 25 pi 26#include "../base/base_uses.f90" 27 28 IMPLICIT NONE 29 30 PRIVATE 31 32 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 's_contract_shg' 33 34! *** Public subroutines *** 35 PUBLIC :: s_overlap_ab, s_overlap_abb, s_coulomb_ab, s_verf_ab, s_verfc_ab, & 36 s_vgauss_ab, s_gauss_ab, s_ra2m_ab 37 38 PUBLIC :: contract_sint_ab_clow, contract_sint_ab_chigh, contract_s_overlap_aba, & 39 contract_s_overlap_abb, contract_s_ra2m_ab 40 41CONTAINS 42 43! ************************************************************************************************** 44!> \brief calculates the uncontracted, not normalized [s|s] overlap 45!> \param la_max maximal l quantum number on a 46!> \param npgfa number of primitive Gaussian on a 47!> \param zeta set of exponents on a 48!> \param lb_max maximal l quantum number on b 49!> \param npgfb number of primitive Gaussian on a 50!> \param zetb set of exponents on a 51!> \param rab distance vector between a and b 52!> \param s uncontracted overlap of s functions 53!> \param calculate_forces ... 54! ************************************************************************************************** 55 SUBROUTINE s_overlap_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, rab, s, calculate_forces) 56 57 INTEGER, INTENT(IN) :: la_max, npgfa 58 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 59 INTEGER, INTENT(IN) :: lb_max, npgfb 60 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 61 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 62 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s 63 LOGICAL, INTENT(IN) :: calculate_forces 64 65 CHARACTER(len=*), PARAMETER :: routineN = 's_overlap_ab', routineP = moduleN//':'//routineN 66 67 INTEGER :: ids, ipgfa, jpgfb, ndev 68 REAL(KIND=dp) :: a, b, rab2, xhi, zet 69 70 ! Distance of the centers a and b 71 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 72 ndev = 0 73 IF (calculate_forces) ndev = 1 74 ! Loops over all pairs of primitive Gaussian-type functions 75 DO ipgfa = 1, npgfa 76 DO jpgfb = 1, npgfb 77 78 ! Distance Screening !maybe later 79 80 ! Calculate some prefactors 81 a = zeta(ipgfa) 82 b = zetb(jpgfb) 83 zet = a + b 84 xhi = a*b/zet 85 86 ! [s|s] integral 87 s(ipgfa, jpgfb, 1) = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2) 88 89 DO ids = 2, la_max + lb_max + ndev + 1 90 s(ipgfa, jpgfb, ids) = -xhi*s(ipgfa, jpgfb, ids - 1) 91 ENDDO 92 93 END DO 94 END DO 95 96 END SUBROUTINE s_overlap_ab 97 98! ************************************************************************************************** 99!> \brief calculates [s|ra^n|s] integrals for [aba] and the [s|rb^n|s] 100!> integrals for [abb] 101!> \param la_max maximal l quantum number on a, orbital basis 102!> \param npgfa number of primitive Gaussian on a, orbital basis 103!> \param zeta set of exponents on a, orbital basis 104!> \param lb_max maximal l quantum number on b, orbital basis 105!> \param npgfb number of primitive Gaussian on a, orbital basis 106!> \param zetb set of exponents on b, orbital basis 107!> \param lcb_max maximal l quantum number of aux basis on b 108!> \param npgfcb number of primitive Gaussian on b. aux basis 109!> \param zetcb set of exponents on b, aux basis 110!> \param rab distance vector between a and b 111!> \param s uncontracted [s|r^n|s] integrals 112!> \param calculate_forces ... 113! ************************************************************************************************** 114 SUBROUTINE s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, & 115 rab, s, calculate_forces) 116 117 INTEGER, INTENT(IN) :: la_max, npgfa 118 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 119 INTEGER, INTENT(IN) :: lb_max, npgfb 120 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 121 INTEGER, INTENT(IN) :: lcb_max, npgfcb 122 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetcb 123 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 124 REAL(KIND=dp), DIMENSION(:, :, :, :, :), & 125 INTENT(INOUT) :: s 126 LOGICAL, INTENT(IN) :: calculate_forces 127 128 CHARACTER(len=*), PARAMETER :: routineN = 's_overlap_abb', routineP = moduleN//':'//routineN 129 130 INTEGER :: i, ids, il, ipgfa, j, jpgfb, kpgfb, & 131 lbb_max, lmax, n, ndev, nds, nfac, nl 132 REAL(KIND=dp) :: a, b, dfsr_int, exp_rab2, k, pfac, & 133 prefac, rab2, sqrt_pi3, sqrt_zet, & 134 sr_int, temp, xhi, zet 135 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dsr_int, dtemp 136 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coeff_srs 137 138 ! Distance of the centers a and b 139 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 140 ndev = 0 141 IF (calculate_forces) ndev = 1 142 143 lbb_max = lb_max + lcb_max 144 nl = INT(lbb_max/2) 145 IF (lb_max == 0 .OR. lcb_max == 0) nl = 0 146 lmax = la_max + lbb_max 147 148 ALLOCATE (dtemp(nl + 1), dsr_int(nl + 1)) 149 ALLOCATE (coeff_srs(nl + 1, nl + 1, nl + 1)) 150 IF (nl > 5) CALL get_prefac_sabb(nl, coeff_srs) 151 152 sqrt_pi3 = SQRT(pi**3) 153 154 ! Loops over all pairs of primitive Gaussian-type functions 155 DO kpgfb = 1, npgfcb 156 DO jpgfb = 1, npgfb 157 DO ipgfa = 1, npgfa 158 159 !Calculate some prefactors 160 a = zeta(ipgfa) 161 b = zetb(jpgfb) + zetcb(kpgfb) 162 163 zet = a + b 164 xhi = a*b/zet 165 exp_rab2 = EXP(-xhi*rab2) 166 167 pfac = a**2/zet 168 sqrt_zet = SQRT(zet) 169 170 DO il = 0, nl 171 nds = lmax - 2*il + ndev + 1 172 SELECT CASE (il) 173 CASE (0) 174 ! [s|s] integral 175 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2 176 DO ids = 2, nds 177 n = ids - 1 178 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, kpgfb, il + 1, 1) 179 ENDDO 180 CASE (1) 181 ![s|r^2|s] integral 182 sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp 183 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int 184 k = sqrt_pi3*a**2/sqrt_zet**7 185 DO ids = 2, nds 186 n = ids - 1 187 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*exp_rab2*sr_int & 188 + n*(-xhi)**(n - 1)*k*exp_rab2 189 ENDDO 190 CASE (2) 191 ![s|r^4|s] integral 192 prefac = sqrt_pi3/4.0_dp/sqrt_zet**7 193 temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2 194 sr_int = prefac*temp 195 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int 196 !** derivatives 197 k = sqrt_pi3*a**4/sqrt_zet**11 198 dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2) 199 DO ids = 2, nds 200 n = ids - 1 201 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 202 dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 203 dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2 204 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) 205 ENDDO 206 CASE (3) 207 ![s|r^6|s] integral 208 prefac = sqrt_pi3/8.0_dp/sqrt_zet**9 209 temp = 105.0_dp + 210.0_dp*pfac*rab2 210 temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3 211 sr_int = prefac*temp 212 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int 213 !** derivatives 214 k = sqrt_pi3*a**6/sqrt_zet**15 215 dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 & 216 + 24.0_dp*pfac**3*rab2**2) 217 dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2) 218 DO ids = 2, nds 219 n = ids - 1 220 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 221 dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 222 dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) & 223 *exp_rab2*dsr_int(2) 224 dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)*(-xhi)**(n - 3)*k*exp_rab2 225 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) & 226 + dtemp(3) + dtemp(4) 227 ENDDO 228 CASE (4) 229 ![s|r^8|s] integral 230 prefac = sqrt_pi3/16.0_dp/sqrt_zet**11 231 temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2 232 temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4 233 sr_int = prefac*temp 234 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int 235 !** derivatives 236 k = sqrt_pi3*a**8/sqrt_zet**19 237 dsr_int(1) = 2520.0_dp*pfac + 3024.0_dp*pfac**2*rab2 238 dsr_int(1) = dsr_int(1) + 864.0_dp*pfac**3*rab2**2 & 239 + 64.0_dp*pfac**4*rab2**3 240 dsr_int(1) = prefac*dsr_int(1) 241 dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2 242 dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2 243 dsr_int(2) = prefac*dsr_int(2) 244 dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2 245 dsr_int(3) = prefac*dsr_int(3) 246 DO ids = 2, nds 247 n = ids - 1 248 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 249 dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 250 dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) & 251 *exp_rab2*dsr_int(2) 252 dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) & 253 *exp_rab2*dsr_int(3) 254 dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) & 255 *k*exp_rab2 256 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) & 257 + dtemp(4) + dtemp(5) 258 ENDDO 259 CASE (5) 260 ![s|r^10|s] integral 261 prefac = sqrt_pi3/32.0_dp/sqrt_zet**13 262 temp = 10395.0_dp + 34650.0_dp*pfac*rab2 263 temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3 264 temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5 265 sr_int = prefac*temp 266 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int 267 !** derivatives 268 k = sqrt_pi3*a**10/sqrt_zet**23 269 dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2 270 dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2 271 dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3 272 dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4 273 dsr_int(1) = prefac*dsr_int(1) 274 dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2 275 dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2 276 dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3 277 dsr_int(2) = prefac*dsr_int(2) 278 dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2 279 dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2 280 dsr_int(3) = prefac*dsr_int(3) 281 dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2 282 dsr_int(4) = prefac*dsr_int(4) 283 DO ids = 2, nds 284 n = ids - 1 285 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 286 dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 287 dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) & 288 *exp_rab2*dsr_int(2) 289 dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) & 290 *exp_rab2*dsr_int(3) 291 dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) & 292 *exp_rab2*dsr_int(4) 293 dtemp(6) = REAL(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) & 294 *k*exp_rab2 295 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) & 296 + dtemp(4) + dtemp(5) + dtemp(6) 297 ENDDO 298 CASE DEFAULT 299 !*** general formula; factor 1.5-2 slower than explicit expressions 300 prefac = exp_rab2/sqrt_zet**(2*il + 3) 301 sr_int = 0.0_dp 302 DO i = 0, il 303 sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i 304 ENDDO 305 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = prefac*sr_int 306 DO ids = 2, nds 307 n = ids - 1 308 nfac = 1 309 dfsr_int = (-xhi)**n*sr_int 310 DO j = 1, il 311 temp = 0.0_dp 312 DO i = j, il 313 temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j) 314 ENDDO 315 nfac = nfac*(n - j + 1) 316 dfsr_int = dfsr_int + temp*REAL(nfac, dp)/fac(j)*(-xhi)**(n - j) 317 ENDDO 318 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = prefac*dfsr_int 319 ENDDO 320 321 END SELECT 322 323 ENDDO 324 325 END DO 326 END DO 327 END DO 328 329 DEALLOCATE (dtemp, dsr_int) 330 DEALLOCATE (coeff_srs) 331 332 END SUBROUTINE s_overlap_abb 333 334! ************************************************************************************************** 335!> \brief calculates the uncontracted, not normalized [0a|ra^(2m)|0b] two-center integral, 336!> where ra = r-Ra and Ra center of a 337!> \param la_max maximal l quantum number on a 338!> \param npgfa number of primitive Gaussian on a 339!> \param zeta set of exponents on a 340!> \param lb_max maximal l quantum number on b 341!> \param npgfb number of primitive Gaussian on a 342!> \param zetb set of exponents on a 343!> \param m exponent of the ra operator 344!> \param rab distance vector between a and b 345!> \param s ... 346!> \param calculate_forces ... 347! ************************************************************************************************** 348 SUBROUTINE s_ra2m_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, m, rab, s, calculate_forces) 349 350 INTEGER, INTENT(IN) :: la_max, npgfa 351 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 352 INTEGER, INTENT(IN) :: lb_max, npgfb 353 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 354 INTEGER, INTENT(IN) :: m 355 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 356 REAL(KIND=dp), DIMENSION(:, :, :, :), & 357 INTENT(INOUT) :: s 358 LOGICAL, INTENT(IN) :: calculate_forces 359 360 CHARACTER(len=*), PARAMETER :: routineN = 's_ra2m_ab', routineP = moduleN//':'//routineN 361 362 INTEGER :: i, ids, il, ipgfa, j, jpgfb, n, ndev, & 363 nds, nfac 364 REAL(KIND=dp) :: a, b, dfsr_int, exp_rab2, k, pfac, & 365 prefac, rab2, sqrt_pi3, sqrt_zet, & 366 sr_int, temp, xhi, zet 367 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dsr_int, dtemp 368 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coeff_srs 369 370 ! Distance of the centers a and b 371 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 372 ndev = 0 373 IF (calculate_forces) ndev = 1 374 375 ALLOCATE (dtemp(m + 1), dsr_int(m + 1)) 376 ALLOCATE (coeff_srs(m + 1, m + 1, m + 1)) 377 CALL get_prefac_sabb(m, coeff_srs) 378 !IF(m > 5) CALL get_prefac_sabb(m, coeff_srs) 379 sqrt_pi3 = SQRT(pi**3) 380 381 ! Loops over all pairs of primitive Gaussian-type functions 382 DO ipgfa = 1, npgfa 383 DO jpgfb = 1, npgfb 384 385 ! Calculate some prefactors 386 a = zeta(ipgfa) 387 b = zetb(jpgfb) 388 zet = a + b 389 xhi = a*b/zet 390 exp_rab2 = EXP(-xhi*rab2) 391 392 sqrt_zet = SQRT(zet) 393 pfac = b**2/zet 394 395 nds = la_max + lb_max + ndev + 1 396 DO il = 0, m 397 SELECT CASE (il) 398 CASE (0) 399 ! [s|s] integral 400 s(ipgfa, jpgfb, m - il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2 401 DO ids = 2, nds 402 n = ids - 1 403 s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, m - il + 1, 1) 404 ENDDO 405 CASE (1) 406 ![s|r^2|s] integral 407 sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp 408 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int 409 k = sqrt_pi3*b**2/sqrt_zet**7 410 DO ids = 2, nds 411 n = ids - 1 412 s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*exp_rab2*sr_int & 413 + n*(-xhi)**(n - 1)*k*exp_rab2 414 ENDDO 415 CASE (2) 416 ![s|r^4|s] integral 417 prefac = sqrt_pi3/4.0_dp/sqrt_zet**7 418 temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2 419 sr_int = prefac*temp 420 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int 421 !** derivatives 422 k = sqrt_pi3*b**4/sqrt_zet**11 423 dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2) 424 DO ids = 2, nds 425 n = ids - 1 426 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 427 dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 428 dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2 429 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) 430 ENDDO 431 CASE (3) 432 ![s|r^6|s] integral 433 prefac = sqrt_pi3/8.0_dp/sqrt_zet**9 434 temp = 105.0_dp + 210.0_dp*pfac*rab2 435 temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3 436 sr_int = prefac*temp 437 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int 438 !** derivatives 439 k = sqrt_pi3*b**6/sqrt_zet**15 440 dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 & 441 + 24.0_dp*pfac**3*rab2**2) 442 dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2) 443 DO ids = 2, nds 444 n = ids - 1 445 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 446 dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 447 dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) & 448 *exp_rab2*dsr_int(2) 449 dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)*(-xhi)**(n - 3)*k*exp_rab2 450 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) & 451 + dtemp(3) + dtemp(4) 452 ENDDO 453 CASE (4) 454 ![s|r^8|s] integral 455 prefac = sqrt_pi3/16.0_dp/sqrt_zet**11 456 temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2 457 temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4 458 sr_int = prefac*temp 459 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int 460 !** derivatives 461 k = sqrt_pi3*b**8/sqrt_zet**19 462 dsr_int(1) = 2520.0_dp*pfac + 3024.0_dp*pfac**2*rab2 463 dsr_int(1) = dsr_int(1) + 864.0_dp*pfac**3*rab2**2 & 464 + 64.0_dp*pfac**4*rab2**3 465 dsr_int(1) = prefac*dsr_int(1) 466 dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2 467 dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2 468 dsr_int(2) = prefac*dsr_int(2) 469 dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2 470 dsr_int(3) = prefac*dsr_int(3) 471 DO ids = 2, nds 472 n = ids - 1 473 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 474 dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 475 dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) & 476 *exp_rab2*dsr_int(2) 477 dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) & 478 *exp_rab2*dsr_int(3) 479 dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) & 480 *k*exp_rab2 481 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) & 482 + dtemp(4) + dtemp(5) 483 ENDDO 484 CASE (5) 485 ![s|r^10|s] integral 486 prefac = sqrt_pi3/32.0_dp/sqrt_zet**13 487 temp = 10395.0_dp + 34650.0_dp*pfac*rab2 488 temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3 489 temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5 490 sr_int = prefac*temp 491 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int 492 !** derivatives 493 k = sqrt_pi3*b**10/sqrt_zet**23 494 dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2 495 dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2 496 dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3 497 dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4 498 dsr_int(1) = prefac*dsr_int(1) 499 dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2 500 dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2 501 dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3 502 dsr_int(2) = prefac*dsr_int(2) 503 dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2 504 dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2 505 dsr_int(3) = prefac*dsr_int(3) 506 dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2 507 dsr_int(4) = prefac*dsr_int(4) 508 DO ids = 2, nds 509 n = ids - 1 510 dtemp(1) = (-xhi)**n*exp_rab2*sr_int 511 dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1) 512 dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) & 513 *exp_rab2*dsr_int(2) 514 dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) & 515 *exp_rab2*dsr_int(3) 516 dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) & 517 *exp_rab2*dsr_int(4) 518 dtemp(6) = REAL(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) & 519 *k*exp_rab2 520 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) & 521 + dtemp(4) + dtemp(5) + dtemp(6) 522 ENDDO 523 CASE DEFAULT 524 prefac = exp_rab2/sqrt_zet**(2*il + 3) 525 sr_int = 0.0_dp 526 DO i = 0, il 527 sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i 528 ENDDO 529 s(ipgfa, jpgfb, m - il + 1, 1) = prefac*sr_int 530 DO ids = 2, nds 531 n = ids - 1 532 nfac = 1 533 dfsr_int = (-xhi)**n*sr_int 534 DO j = 1, il 535 temp = 0.0_dp 536 DO i = j, il 537 temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j) 538 ENDDO 539 nfac = nfac*(n - j + 1) 540 dfsr_int = dfsr_int + temp*REAL(nfac, dp)/fac(j)*(-xhi)**(n - j) 541 ENDDO 542 s(ipgfa, jpgfb, m - il + 1, ids) = prefac*dfsr_int 543 ENDDO 544 END SELECT 545 ENDDO 546 ENDDO 547 ENDDO 548 549 DEALLOCATE (coeff_srs) 550 DEALLOCATE (dtemp, dsr_int) 551 552 END SUBROUTINE s_ra2m_ab 553 554! ************************************************************************************************** 555!> \brief prefactor for the general formula to calculate the (0a|0b|0b) overlap integrals 556!> \param nl ... 557!> \param prefac ... 558! ************************************************************************************************** 559 SUBROUTINE get_prefac_sabb(nl, prefac) 560 INTEGER, INTENT(IN) :: nl 561 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: prefac 562 563 CHARACTER(len=*), PARAMETER :: routineN = 'get_prefac_sabb', & 564 routineP = moduleN//':'//routineN 565 566 INTEGER :: il, j, k 567 REAL(KIND=dp) :: sqrt_pi3, temp 568 569 sqrt_pi3 = SQRT(pi**3) 570 571 DO il = 0, nl 572 temp = dfac(2*il + 1)*sqrt_pi3*fac(il)/2.0_dp**il 573 DO j = 0, il 574 DO k = j, il 575 prefac(k + 1, j + 1, il + 1) = temp*2.0_dp**k/dfac(2*k + 1)/fac(il - k)/fac(k - j) 576 ENDDO 577 ENDDO 578 ENDDO 579 580 END SUBROUTINE get_prefac_sabb 581 582! ************************************************************************************************** 583!> \brief calculates the uncontracted, not normalized [s|1/r12|s] two-center coulomb integral 584!> \param la_max maximal l quantum number on a 585!> \param npgfa number of primitive Gaussian on a 586!> \param zeta set of exponents on a 587!> \param lb_max maximal l quantum number on b 588!> \param npgfb number of primitive Gaussian on a 589!> \param zetb set of exponents on a 590!> \param omega parameter not needed, but given for the sake of the abstract interface 591!> \param rab distance vector between a and b 592!> \param v uncontracted coulomb integral of s functions 593!> \param calculate_forces ... 594! ************************************************************************************************** 595 SUBROUTINE s_coulomb_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces) 596 597 INTEGER, INTENT(IN) :: la_max, npgfa 598 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 599 INTEGER, INTENT(IN) :: lb_max, npgfb 600 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 601 REAL(KIND=dp), INTENT(IN) :: omega 602 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 603 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 604 LOGICAL, INTENT(IN) :: calculate_forces 605 606 CHARACTER(len=*), PARAMETER :: routineN = 's_coulomb_ab', routineP = moduleN//':'//routineN 607 608 INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax 609 REAL(KIND=dp) :: a, b, dummy, prefac, rab2, T, xhi, zet 610 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 611 612 dummy = omega 613 614 ! Distance of the centers a and b 615 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 616 ndev = 0 617 IF (calculate_forces) ndev = 1 618 nmax = la_max + lb_max + ndev + 1 619 ALLOCATE (f(0:nmax)) 620 ! Loops over all pairs of primitive Gaussian-type functions 621 DO ipgfa = 1, npgfa 622 DO jpgfb = 1, npgfb 623 624 ! Calculate some prefactors 625 a = zeta(ipgfa) 626 b = zetb(jpgfb) 627 zet = a + b 628 xhi = a*b/zet 629 prefac = 2.0_dp*SQRT(pi**5)/(a*b)/SQRT(zet) 630 T = xhi*rab2 631 CALL fgamma(nmax - 1, T, f) 632 633 DO ids = 1, nmax 634 n = ids - 1 635 v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*f(n) 636 ENDDO 637 638 END DO 639 END DO 640 DEALLOCATE (f) 641 642 END SUBROUTINE s_coulomb_ab 643 644! ************************************************************************************************** 645!> \brief calculates the uncontracted, not normalized [s|1/erf(omega*r12)/r12|s] two-center integral 646!> \param la_max maximal l quantum number on a 647!> \param npgfa number of primitive Gaussian on a 648!> \param zeta set of exponents on a 649!> \param lb_max maximal l quantum number on b 650!> \param npgfb number of primitive Gaussian on a 651!> \param zetb set of exponents on a 652!> \param omega parameter in the operator 653!> \param rab distance vector between a and b 654!> \param v uncontracted erf(r)/r integral of s functions 655!> \param calculate_forces ... 656! ************************************************************************************************** 657 SUBROUTINE s_verf_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces) 658 659 INTEGER, INTENT(IN) :: la_max, npgfa 660 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 661 INTEGER, INTENT(IN) :: lb_max, npgfb 662 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 663 REAL(KIND=dp), INTENT(IN) :: omega 664 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 665 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 666 LOGICAL, INTENT(IN) :: calculate_forces 667 668 CHARACTER(len=*), PARAMETER :: routineN = 's_verf_ab', routineP = moduleN//':'//routineN 669 670 INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax 671 REAL(KIND=dp) :: a, Arg, b, comega, prefac, rab2, T, xhi, & 672 zet 673 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 674 675 ! Distance of the centers a and b 676 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 677 ndev = 0 678 IF (calculate_forces) ndev = 1 679 nmax = la_max + lb_max + ndev + 1 680 ALLOCATE (f(0:nmax)) 681 ! Loops over all pairs of primitive Gaussian-type functions 682 DO ipgfa = 1, npgfa 683 DO jpgfb = 1, npgfb 684 685 ! Calculate some prefactors 686 a = zeta(ipgfa) 687 b = zetb(jpgfb) 688 zet = a + b 689 xhi = a*b/zet 690 comega = omega**2/(omega**2 + xhi) 691 prefac = 2.0_dp*SQRT(pi**5)*SQRT(comega)/(a*b)/SQRT(zet) 692 T = xhi*rab2 693 Arg = comega*T 694 CALL fgamma(nmax - 1, Arg, f) 695 696 DO ids = 1, nmax 697 n = ids - 1 698 v(ipgfa, jpgfb, ids) = prefac*(-xhi*comega)**n*f(n) 699 ENDDO 700 701 END DO 702 END DO 703 DEALLOCATE (f) 704 705 END SUBROUTINE s_verf_ab 706 707! ************************************************************************************************** 708!> \brief calculates the uncontracted, not normalized [s|1/erfc(omega*r12)/r12|s] two-center 709!> integral 710!> \param la_max maximal l quantum number on a 711!> \param npgfa number of primitive Gaussian on a 712!> \param zeta set of exponents on a 713!> \param lb_max maximal l quantum number on b 714!> \param npgfb number of primitive Gaussian on a 715!> \param zetb set of exponents on a 716!> \param omega parameter in the operator 717!> \param rab distance vector between a and b 718!> \param v uncontracted erf(r)/r integral of s functions 719!> \param calculate_forces ... 720! ************************************************************************************************** 721 SUBROUTINE s_verfc_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces) 722 723 INTEGER, INTENT(IN) :: la_max, npgfa 724 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 725 INTEGER, INTENT(IN) :: lb_max, npgfb 726 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 727 REAL(KIND=dp), INTENT(IN) :: omega 728 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 729 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 730 LOGICAL, INTENT(IN) :: calculate_forces 731 732 CHARACTER(len=*), PARAMETER :: routineN = 's_verfc_ab', routineP = moduleN//':'//routineN 733 734 INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax 735 REAL(KIND=dp) :: a, b, comega, comegaT, prefac, rab2, T, & 736 xhi, zet 737 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf 738 739 ! Distance of the centers a and b 740 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 741 ndev = 0 742 IF (calculate_forces) ndev = 1 743 nmax = la_max + lb_max + ndev + 1 744 ALLOCATE (fv(0:nmax), fverf(0:nmax)) 745 ! Loops over all pairs of primitive Gaussian-type functions 746 DO ipgfa = 1, npgfa 747 DO jpgfb = 1, npgfb 748 749 ! Calculate some prefactors 750 a = zeta(ipgfa) 751 b = zetb(jpgfb) 752 zet = a + b 753 xhi = a*b/zet 754 comega = omega**2/(omega**2 + xhi) 755 prefac = 2.0_dp*SQRT(pi**5)/(a*b)/SQRT(zet) 756 T = xhi*rab2 757 comegaT = comega*T 758 CALL fgamma(nmax - 1, T, fv) 759 CALL fgamma(nmax - 1, comegaT, fverf) 760 761 DO ids = 1, nmax 762 n = ids - 1 763 v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*(fv(n) - SQRT(comega)*comega**n*fverf(n)) 764 ENDDO 765 766 END DO 767 END DO 768 DEALLOCATE (fv, fverf) 769 770 END SUBROUTINE s_verfc_ab 771 772! ************************************************************************************************** 773!> \brief calculates the uncontracted, not normalized [s|exp(-omega*r12**2)/r12|s] two-center 774!> integral 775!> \param la_max maximal l quantum number on a 776!> \param npgfa number of primitive Gaussian on a 777!> \param zeta set of exponents on a 778!> \param lb_max maximal l quantum number on b 779!> \param npgfb number of primitive Gaussian on a 780!> \param zetb set of exponents on a 781!> \param omega parameter in the operator 782!> \param rab distance vector between a and b 783!> \param v uncontracted exp(-omega*r**2)/r integral of s functions 784!> \param calculate_forces ... 785! ************************************************************************************************** 786 SUBROUTINE s_vgauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces) 787 788 INTEGER, INTENT(IN) :: la_max, npgfa 789 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 790 INTEGER, INTENT(IN) :: lb_max, npgfb 791 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 792 REAL(KIND=dp), INTENT(IN) :: omega 793 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 794 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 795 LOGICAL, INTENT(IN) :: calculate_forces 796 797 CHARACTER(len=*), PARAMETER :: routineN = 's_vgauss_ab', routineP = moduleN//':'//routineN 798 799 INTEGER :: ids, ipgfa, j, jpgfb, n, ndev, nmax 800 REAL(KIND=dp) :: a, b, eta, etaT, expT, oeta, prefac, & 801 rab2, T, xeta, xhi, zet 802 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 803 804 ! Distance of the centers a and b 805 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 806 ndev = 0 807 IF (calculate_forces) ndev = 1 808 nmax = la_max + lb_max + ndev + 1 809 ALLOCATE (f(0:nmax)) 810 ! Loops over all pairs of primitive Gaussian-type functions 811 v = 0.0_dp 812 DO ipgfa = 1, npgfa 813 DO jpgfb = 1, npgfb 814 815 ! Calculate some prefactors 816 a = zeta(ipgfa) 817 b = zetb(jpgfb) 818 zet = a + b 819 xhi = a*b/zet 820 eta = xhi/(xhi + omega) 821 oeta = omega*eta 822 xeta = xhi*eta 823 T = xhi*rab2 824 expT = EXP(-omega/(omega + xhi)*T) 825 prefac = 2.0_dp*SQRT(pi**5/zet**3)/(xhi + omega)*expT 826 etaT = eta*T 827 CALL fgamma(nmax - 1, etaT, f) 828 829 DO ids = 1, nmax 830 n = ids - 1 831 DO j = 0, n 832 v(ipgfa, jpgfb, ids) = v(ipgfa, jpgfb, ids) & 833 + prefac*fac(n)/fac(j)/fac(n - j)*(-oeta)**(n - j)*(-xeta)**j*f(j) 834 ENDDO 835 ENDDO 836 837 END DO 838 END DO 839 DEALLOCATE (f) 840 841 END SUBROUTINE s_vgauss_ab 842 843! ************************************************************************************************** 844!> \brief calculates the uncontracted, not normalized [s|exp(-omega*r12**2)|s] two-center 845!> integral 846!> \param la_max maximal l quantum number on a 847!> \param npgfa number of primitive Gaussian on a 848!> \param zeta set of exponents on a 849!> \param lb_max maximal l quantum number on b 850!> \param npgfb number of primitive Gaussian on a 851!> \param zetb set of exponents on a 852!> \param omega parameter in the operator 853!> \param rab distance vector between a and b 854!> \param v uncontracted exp(-omega*r**2) integral of s functions 855!> \param calculate_forces ... 856! ************************************************************************************************** 857 SUBROUTINE s_gauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces) 858 859 INTEGER, INTENT(IN) :: la_max, npgfa 860 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 861 INTEGER, INTENT(IN) :: lb_max, npgfb 862 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 863 REAL(KIND=dp), INTENT(IN) :: omega 864 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 865 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 866 LOGICAL, INTENT(IN) :: calculate_forces 867 868 CHARACTER(len=*), PARAMETER :: routineN = 's_gauss_ab', routineP = moduleN//':'//routineN 869 870 INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax 871 REAL(KIND=dp) :: a, b, eta, expT, oeta, prefac, rab2, T, & 872 xhi, zet 873 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 874 875 ! Distance of the centers a and b 876 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 877 ndev = 0 878 IF (calculate_forces) ndev = 1 879 nmax = la_max + lb_max + ndev + 1 880 ALLOCATE (f(0:nmax)) 881 ! Loops over all pairs of primitive Gaussian-type functions 882 DO ipgfa = 1, npgfa 883 DO jpgfb = 1, npgfb 884 885 ! Calculate some prefactors 886 a = zeta(ipgfa) 887 b = zetb(jpgfb) 888 zet = a + b 889 xhi = a*b/zet 890 eta = xhi/(xhi + omega) 891 oeta = omega*eta 892 T = xhi*rab2 893 expT = EXP(-omega/(omega + xhi)*T) 894 prefac = pi**3/SQRT(zet**3)/SQRT((xhi + omega)**3)*expT 895 896 DO ids = 1, nmax 897 n = ids - 1 898 v(ipgfa, jpgfb, ids) = prefac*(-oeta)**n 899 ENDDO 900 901 END DO 902 END DO 903 DEALLOCATE (f) 904 905 END SUBROUTINE s_gauss_ab 906 907! ************************************************************************************************** 908!> \brief Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; 909!> this routine is more efficient for uncontracted basis sets (clow), e.g. for ri basis sets 910!> \param la set of l quantum numbers for a 911!> \param npgfa number of primitive Gaussian on a 912!> \param nshella number of shells for a 913!> \param scona_shg SHG contraction/normalization matrix for a 914!> \param lb set of l quantum numbers for b 915!> \param npgfb number of primitive Gaussian on b 916!> \param nshellb number of shells for b 917!> \param sconb_shg SHG contraction/normalization matrix for b 918!> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives 919!> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives 920!> \param calculate_forces ... 921! ************************************************************************************************** 922 SUBROUTINE contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, & 923 swork, swork_cont, calculate_forces) 924 925 INTEGER, DIMENSION(:), INTENT(IN) :: la 926 INTEGER, INTENT(IN) :: npgfa, nshella 927 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg 928 INTEGER, DIMENSION(:), INTENT(IN) :: lb 929 INTEGER, INTENT(IN) :: npgfb, nshellb 930 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg 931 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork 932 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont 933 LOGICAL, INTENT(IN) :: calculate_forces 934 935 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_sint_ab_clow', & 936 routineP = moduleN//':'//routineN 937 938 INTEGER :: ids, ids_start, ipgfa, ishella, jpgfb, & 939 jshellb, lai, lbj, ndev, nds 940 941 ndev = 0 942 IF (calculate_forces) ndev = 1 943 944 swork_cont = 0.0_dp 945 DO ishella = 1, nshella 946 lai = la(ishella) 947 DO jshellb = 1, nshellb 948 lbj = lb(jshellb) 949 nds = lai + lbj + 1 950 ids_start = nds - MIN(lai, lbj) 951 DO ipgfa = 1, npgfa 952 DO jpgfb = 1, npgfb 953 DO ids = ids_start, nds + ndev 954 swork_cont(ids, ishella, jshellb) = swork_cont(ids, ishella, jshellb) & 955 + scona_shg(ipgfa, ishella) & 956 *sconb_shg(jpgfb, jshellb) & 957 *swork(ipgfa, jpgfb, ids) 958 ENDDO 959 ENDDO 960 ENDDO 961 ENDDO 962 ENDDO 963 964 END SUBROUTINE contract_sint_ab_clow 965 966! ************************************************************************************************** 967!> \brief Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; 968!> this routine is more efficient for highly contracted basis sets (chigh) 969!> \param npgfa number of primitive Gaussian on a 970!> \param nshella number of shells for a 971!> \param scona SHG contraction/normalization matrix for a 972!> \param npgfb number of primitive Gaussian on b 973!> \param nshellb number of shells for b 974!> \param sconb SHG contraction/normalization matrix for b 975!> \param nds maximal derivative of [s|O(r12)|s] with respect to rab2 976!> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives 977!> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives 978! ************************************************************************************************** 979 SUBROUTINE contract_sint_ab_chigh(npgfa, nshella, scona, & 980 npgfb, nshellb, sconb, & 981 nds, swork, swork_cont) 982 983 INTEGER, INTENT(IN) :: npgfa, nshella 984 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona 985 INTEGER, INTENT(IN) :: npgfb, nshellb 986 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb 987 INTEGER, INTENT(IN) :: nds 988 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork 989 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont 990 991 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_sint_ab_chigh', & 992 routineP = moduleN//':'//routineN 993 994 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc 995 996 swork_cont = 0.0_dp 997 ALLOCATE (work_pc(npgfb, nds, nshella)) 998 999 CALL dgemm("T", "N", npgfb*nds, nshella, npgfa, 1._dp, swork, npgfa, scona, npgfa, & 1000 0.0_dp, work_pc, npgfb*nds) 1001 CALL dgemm("T", "N", nds*nshella, nshellb, npgfb, 1._dp, work_pc, npgfb, sconb, npgfb, & 1002 0.0_dp, swork_cont, nds*nshella) 1003 1004 DEALLOCATE (work_pc) 1005 1006 END SUBROUTINE contract_sint_ab_chigh 1007 1008! ************************************************************************************************** 1009!> \brief Contraction, normalization and combinatiorial combination of the [0a|(r-Ra)^(2m)|0b] 1010!> integrals and their scalar derivatives 1011!> \param npgfa number of primitive Gaussian on a 1012!> \param nshella number of shells for a 1013!> \param scon_ra2m contraction matrix on a containg the combinatorial factors 1014!> \param npgfb number of primitive Gaussian on b 1015!> \param nshellb number of shells for b 1016!> \param sconb SHG contraction/normalization matrix for b 1017!> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives 1018!> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives 1019!> \param m exponent in operator (r-Ra)^(2m) 1020!> \param nds_max maximal derivative with respect to rab2 1021! ************************************************************************************************** 1022 SUBROUTINE contract_s_ra2m_ab(npgfa, nshella, scon_ra2m, npgfb, nshellb, sconb, swork, swork_cont, & 1023 m, nds_max) 1024 1025 INTEGER, INTENT(IN) :: npgfa, nshella 1026 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_ra2m 1027 INTEGER, INTENT(IN) :: npgfb, nshellb 1028 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb 1029 REAL(KIND=dp), DIMENSION(:, :, :, :), & 1030 INTENT(INOUT) :: swork 1031 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont 1032 INTEGER, INTENT(IN) :: m, nds_max 1033 1034 CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_s_ra2m_ab', & 1035 routineP = moduleN//':'//routineN 1036 1037 INTEGER :: i, my_m 1038 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc 1039 1040 my_m = m + 1 1041 ALLOCATE (work_pc(npgfb, nds_max, nshella)) 1042 work_pc = 0.0_dp 1043 swork_cont = 0.0_dp 1044 DO i = 1, my_m 1045 CALL dgemm("T", "N", npgfb*nds_max, nshella, npgfa, 1.0_dp, swork(1:npgfa, 1:npgfb, i, 1:nds_max), npgfa, & 1046 scon_ra2m(1:npgfa, i, 1:nshella), npgfa, 1.0_dp, work_pc, npgfb*nds_max) 1047 ENDDO 1048 CALL dgemm("T", "N", nds_max*nshella, nshellb, npgfb, 1.0_dp, work_pc, npgfb, sconb, npgfb, 0.0_dp, & 1049 swork_cont, nds_max*nshella) 1050 1051 DEALLOCATE (work_pc) 1052 1053 END SUBROUTINE contract_s_ra2m_ab 1054 1055! ************************************************************************************************** 1056!> \brief Contraction and normalization of the (abb) overlap 1057!> \param la set of l quantum numbers on a 1058!> \param npgfa number of primitive Gaussians functions on a; orbital basis 1059!> \param nshella number of shells for a; orbital basis 1060!> \param scona_shg SHG contraction/normalization matrix for a; orbital basis 1061!> \param lb set of l quantum numbers on b; orbital basis 1062!> \param npgfb number of primitive Gaussians functions on b; orbital basis 1063!> \param nshellb number of shells for b; orbital basis 1064!> \param lcb set of l quantum numbers on b; ri basis 1065!> \param npgfcb number of primitive Gaussians functions on b; ri basis 1066!> \param nshellcb number of shells for b; ri basis 1067!> \param orbb_index index for orbital basis at B for sconb_mix 1068!> \param rib_index index for ri basis at B for sconb_mix 1069!> \param sconb_mix mixed contraction matrix for orbital and ri basis at B 1070!> \param nl_max related to the parameter m in (a|rb^(2m)|b) 1071!> \param nds_max derivative with respect to rab**2 1072!> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives 1073!> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives 1074!> \param calculate_forces ... 1075! ************************************************************************************************** 1076 SUBROUTINE contract_s_overlap_abb(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, & 1077 lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, & 1078 nl_max, nds_max, swork, swork_cont, calculate_forces) 1079 1080 INTEGER, DIMENSION(:), INTENT(IN) :: la 1081 INTEGER, INTENT(IN) :: npgfa, nshella 1082 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg 1083 INTEGER, DIMENSION(:), INTENT(IN) :: lb 1084 INTEGER, INTENT(IN) :: npgfb, nshellb 1085 INTEGER, DIMENSION(:), INTENT(IN) :: lcb 1086 INTEGER, INTENT(IN) :: npgfcb, nshellcb 1087 INTEGER, DIMENSION(:, :), INTENT(IN) :: orbb_index, rib_index 1088 REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix 1089 INTEGER, INTENT(IN) :: nl_max, nds_max 1090 REAL(KIND=dp), DIMENSION(:, :, :, :, :), & 1091 INTENT(IN) :: swork 1092 REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), & 1093 INTENT(INOUT) :: swork_cont 1094 LOGICAL, INTENT(IN) :: calculate_forces 1095 1096 CHARACTER(len=*), PARAMETER :: routineN = 'contract_s_overlap_abb', & 1097 routineP = moduleN//':'//routineN 1098 1099 INTEGER :: forb, fri, ids, ids_start, iil, il, & 1100 ishella, jpgfb, jshellb, kpgfb, & 1101 kshellb, lai, lbb, lbj, lbk, ndev, & 1102 nds, nl 1103 REAL(KIND=dp), ALLOCATABLE, & 1104 DIMENSION(:, :, :, :, :) :: work_ppc 1105 1106 ndev = 0 1107 IF (calculate_forces) ndev = 1 1108 1109 ALLOCATE (work_ppc(npgfb, npgfcb, nl_max, nds_max, nshella)) 1110 work_ppc = 0.0_dp 1111 CALL dgemm("T", "N", npgfb*npgfcb*nl_max*nds_max, nshella, npgfa, 1._dp, swork, npgfa, & 1112 scona_shg, npgfa, 0.0_dp, work_ppc, npgfb*npgfcb*nl_max*nds_max) 1113 swork_cont = 0.0_dp 1114 DO kshellb = 1, nshellcb 1115 lbk = lcb(kshellb) 1116 DO jshellb = 1, nshellb 1117 lbj = lb(jshellb) 1118 nl = INT((lbj + lbk)/2) 1119 IF (lbj == 0 .OR. lbk == 0) nl = 0 1120 DO ishella = 1, nshella 1121 lai = la(ishella) 1122 DO il = 0, nl 1123 lbb = lbj + lbk - 2*il 1124 nds = lai + lbb + 1 1125 ids_start = nds - MIN(lai, lbb) 1126 DO jpgfb = 1, npgfb 1127 forb = orbb_index(jpgfb, jshellb) 1128 DO kpgfb = 1, npgfcb 1129 fri = rib_index(kpgfb, kshellb) 1130 DO ids = ids_start, nds + ndev 1131 DO iil = 0, il 1132 swork_cont(ids, il, ishella, jshellb, kshellb) = & 1133 swork_cont(ids, il, ishella, jshellb, kshellb) & 1134 + sconb_mix(iil + 1, fri, forb, il + 1)*work_ppc(jpgfb, kpgfb, il - iil + 1, ids, ishella) 1135 ENDDO 1136 ENDDO 1137 ENDDO 1138 ENDDO 1139 ENDDO 1140 ENDDO 1141 ENDDO 1142 ENDDO 1143 DEALLOCATE (work_ppc) 1144 1145 END SUBROUTINE contract_s_overlap_abb 1146 1147! ************************************************************************************************** 1148!> \brief Contraction and normalization of the (aba) overlap 1149!> \param la set of l quantum numbers on a; orbital basis 1150!> \param npgfa number of primitive Gaussians functions on a; orbital basis 1151!> \param nshella number of shells for a; orbital basis 1152!> \param lb set of l quantum numbers on b; orbital basis 1153!> \param npgfb number of primitive Gaussians functions on b; orbital basis 1154!> \param nshellb number of shells for b; orbital basis 1155!> \param sconb_shg SHG contraction/normalization matrix for b; orbital basis 1156!> \param lca set of l quantum numbers on a; ri basis 1157!> \param npgfca number of primitive Gaussians functions on a; ri basis 1158!> \param nshellca number of shells for a; ri basis 1159!> \param orba_index index for orbital basis at A for scona_mix 1160!> \param ria_index index for ri basis at A for scona_mix 1161!> \param scona_mix mixed contraction matrix for orbital and ri basis at A 1162!> \param nl_max related to the parameter m in (a|ra^(2m)|b) 1163!> \param nds_max maximal derivative with respect to rab**2 1164!> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives 1165!> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives 1166!> \param calculate_forces ... 1167! ************************************************************************************************** 1168 SUBROUTINE contract_s_overlap_aba(la, npgfa, nshella, lb, npgfb, nshellb, sconb_shg, & 1169 lca, npgfca, nshellca, orba_index, ria_index, scona_mix, & 1170 nl_max, nds_max, swork, swork_cont, calculate_forces) 1171 1172 INTEGER, DIMENSION(:), INTENT(IN) :: la 1173 INTEGER, INTENT(IN) :: npgfa, nshella 1174 INTEGER, DIMENSION(:), INTENT(IN) :: lb 1175 INTEGER, INTENT(IN) :: npgfb, nshellb 1176 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg 1177 INTEGER, DIMENSION(:), INTENT(IN) :: lca 1178 INTEGER, INTENT(IN) :: npgfca, nshellca 1179 INTEGER, DIMENSION(:, :), INTENT(IN) :: orba_index, ria_index 1180 REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix 1181 INTEGER, INTENT(IN) :: nl_max, nds_max 1182 REAL(KIND=dp), DIMENSION(:, :, :, :, :), & 1183 INTENT(IN) :: swork 1184 REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), & 1185 INTENT(INOUT) :: swork_cont 1186 LOGICAL, INTENT(IN) :: calculate_forces 1187 1188 CHARACTER(len=*), PARAMETER :: routineN = 'contract_s_overlap_aba', & 1189 routineP = moduleN//':'//routineN 1190 1191 INTEGER :: forb, fri, ids, ids_start, iil, il, & 1192 ipgfa, ishella, jshellb, kpgfa, & 1193 kshella, laa, lai, lak, lbj, ndev, & 1194 nds, nl 1195 REAL(KIND=dp), ALLOCATABLE, & 1196 DIMENSION(:, :, :, :, :) :: work_ppc 1197 1198 ndev = 0 1199 IF (calculate_forces) ndev = 1 1200 1201 ALLOCATE (work_ppc(npgfa, npgfca, nl_max, nds_max, nshellb)) 1202 work_ppc = 0.0_dp 1203 CALL dgemm("T", "N", npgfa*npgfca*nl_max*nds_max, nshellb, npgfb, 1._dp, swork, npgfb, & 1204 sconb_shg, npgfb, 0.0_dp, work_ppc, npgfa*npgfca*nl_max*nds_max) 1205 swork_cont = 0.0_dp 1206 DO kshella = 1, nshellca 1207 lak = lca(kshella) 1208 DO jshellb = 1, nshellb 1209 lbj = lb(jshellb) 1210 DO ishella = 1, nshella 1211 lai = la(ishella) 1212 nl = INT((lai + lak)/2) 1213 IF (lai == 0 .OR. lak == 0) nl = 0 1214 DO il = 0, nl 1215 laa = lai + lak - 2*il 1216 nds = laa + lbj + 1 1217 ids_start = nds - MIN(laa, lbj) 1218 DO kpgfa = 1, npgfca 1219 fri = ria_index(kpgfa, kshella) 1220 DO ipgfa = 1, npgfa 1221 forb = orba_index(ipgfa, ishella) 1222 DO ids = ids_start, nds + ndev 1223 DO iil = 0, il 1224 swork_cont(ids, il, ishella, jshellb, kshella) = & 1225 swork_cont(ids, il, ishella, jshellb, kshella) & 1226 + scona_mix(iil + 1, fri, forb, il + 1)*work_ppc(ipgfa, kpgfa, il - iil + 1, ids, jshellb) 1227 ENDDO 1228 ENDDO 1229 ENDDO 1230 ENDDO 1231 ENDDO 1232 ENDDO 1233 ENDDO 1234 ENDDO 1235 1236 DEALLOCATE (work_ppc) 1237 END SUBROUTINE contract_s_overlap_aba 1238 1239END MODULE s_contract_shg 1240