1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of integrals over Cartesian Gaussian-type functions for different r12 8!> operators: 1/r12, erf(omega*r12/r12), erfc(omega*r12/r12), exp(-omega*r12^2)/r12 and 9!> exp(-omega*r12^2) 10!> \par Literature 11!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986) 12!> R. Ahlrichs, PCCP, 8, 3072 (2006) 13!> \par History 14!> 05.2019 Added the truncated Coulomb operator (A. Bussy) 15!> \par Parameters 16!> - ax,ay,az : Angular momentum index numbers of orbital a. 17!> - cx,cy,cz : Angular momentum index numbers of orbital c. 18!> - coset : Cartesian orbital set pointer. 19!> - dac : Distance between the atomic centers a and c. 20!> - l{a,c} : Angular momentum quantum number of shell a or c. 21!> - l{a,c}_max : Maximum angular momentum quantum number of shell a or c. 22!> - l{a,c}_min : Minimum angular momentum quantum number of shell a or c. 23!> - ncoset : Number of orbitals in a Cartesian orbital set. 24!> - npgf{a,c} : Degree of contraction of shell a or c. 25!> - rac : Distance vector between the atomic centers a and c. 26!> - rac2 : Square of the distance between the atomic centers a and c. 27!> - zet{a,c} : Exponents of the Gaussian-type functions a or c. 28!> - zetp : Reciprocal of the sum of the exponents of orbital a and b. 29!> - zetw : Reciprocal of the sum of the exponents of orbital a and c. 30!> - omega : Parameter in the operator 31!> - r_cutoff : The cutoff radius for the truncated Coulomb operator 32!> \author Dorothea Golze (05.2016) 33! ************************************************************************************************** 34MODULE ai_operators_r12 35 36 USE gamma, ONLY: fgamma => fgamma_0 37 USE kinds, ONLY: dp 38 USE mathconstants, ONLY: fac,& 39 pi 40 USE orbital_pointers, ONLY: coset,& 41 ncoset 42 USE t_c_g0, ONLY: get_lmax_init,& 43 t_c_g0_n 44#include "../base/base_uses.f90" 45 46 IMPLICIT NONE 47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_operators_r12' 48 PRIVATE 49 50 ! *** Public subroutines *** 51 52 PUBLIC :: operator2, cps_coulomb2, cps_verf2, cps_verfc2, cps_vgauss2, cps_gauss2, ab_sint_os, & 53 cps_truncated2 54 55 ABSTRACT INTERFACE 56! ************************************************************************************************** 57!> \brief Interface for the calculation of integrals over s-functions and the s-type auxiliary 58!> integrals using the Obara-Saika (OS) scheme 59!> \param v ... 60!> \param nmax ... 61!> \param zetp ... 62!> \param zetq ... 63!> \param zetw ... 64!> \param rho ... 65!> \param rac2 ... 66!> \param omega ... 67!> \param r_cutoff ... 68! ************************************************************************************************** 69 SUBROUTINE ab_sint_os(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 70 USE kinds, ONLY: dp 71 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 72 INTEGER, INTENT(IN) :: nmax 73 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, & 74 r_cutoff 75 76 END SUBROUTINE ab_sint_os 77 END INTERFACE 78 79CONTAINS 80 81! ************************************************************************************************** 82!> \brief Calculation of the primitive two-center integrals over Cartesian Gaussian-type 83!> functions for different r12 operators. 84!> \param cps_operator2 procedure pointer for the respective operator. The integrals evaluation 85!> differs only in the evaluation of the cartesian primitive s (cps) integrals [s|O(r12)|s] 86!> and auxiliary integrals [s|O(r12)|s]^n. This pointer selects the correct routine. 87!> \param la_max ... 88!> \param npgfa ... 89!> \param zeta ... 90!> \param la_min ... 91!> \param lc_max ... 92!> \param npgfc ... 93!> \param zetc ... 94!> \param lc_min ... 95!> \param omega ... 96!> \param r_cutoff ... 97!> \param rac ... 98!> \param rac2 ... 99!> \param vac matrix storing the integrals 100!> \param v temporary work array 101!> \param maxder maximal derivative 102!> \param vac_plus matrix storing the integrals for highler l-quantum numbers; used to 103!> construct the derivatives 104! ************************************************************************************************** 105 106 SUBROUTINE operator2(cps_operator2, la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, & 107 omega, r_cutoff, rac, rac2, vac, v, maxder, vac_plus) 108 PROCEDURE(ab_sint_os), POINTER :: cps_operator2 109 INTEGER, INTENT(IN) :: la_max, npgfa 110 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 111 INTEGER, INTENT(IN) :: la_min, lc_max, npgfc 112 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc 113 INTEGER, INTENT(IN) :: lc_min 114 REAL(KIND=dp), INTENT(IN) :: omega, r_cutoff 115 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac 116 REAL(KIND=dp), INTENT(IN) :: rac2 117 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vac 118 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 119 INTEGER, INTENT(IN), OPTIONAL :: maxder 120 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: vac_plus 121 122 CHARACTER(len=*), PARAMETER :: routineN = 'operator2', routineP = moduleN//':'//routineN 123 124 INTEGER :: ax, ay, az, coc, cocx, cocy, cocz, cx, & 125 cy, cz, i, ipgf, j, jpgf, la, lc, & 126 maxder_local, n, na, nap, nc, ncp, & 127 nmax, handle 128 REAL(KIND=dp) :: dac, f1, f2, f3, f4, f5, f6, fcx, & 129 fcy, fcz, rho, zetp, zetq, zetw 130 REAL(KIND=dp), DIMENSION(3) :: raw, rcw 131 132 CALL timeset(routineN, handle) 133 134 v = 0.0_dp 135 136 maxder_local = 0 137 IF (PRESENT(maxder)) THEN 138 maxder_local = maxder 139 vac_plus = 0.0_dp 140 END IF 141 142 nmax = la_max + lc_max + 1 143 144 ! *** Calculate the distance of the centers a and c *** 145 146 dac = SQRT(rac2) 147 148 ! *** Loop over all pairs of primitive Gaussian-type functions *** 149 150 na = 0 151 nap = 0 152 153 DO ipgf = 1, npgfa 154 155 nc = 0 156 ncp = 0 157 158 DO jpgf = 1, npgfc 159 160 ! *** Calculate some prefactors *** 161 162 zetp = 1.0_dp/zeta(ipgf) 163 zetq = 1.0_dp/zetc(jpgf) 164 zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf)) 165 166 rho = zeta(ipgf)*zetc(jpgf)*zetw 167 168 ! *** Calculate the basic two-center integrals [s||s]{n} *** 169 170 CALL cps_operator2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 171 172 ! *** Vertical recurrence steps: [s||s] -> [s||c] *** 173 174 IF (lc_max > 0) THEN 175 176 f1 = 0.5_dp*zetq 177 f2 = -rho*zetq 178 179 rcw(:) = -zeta(ipgf)*zetw*rac(:) 180 181 ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) *** 182 183 DO n = 1, nmax - 1 184 v(1, 2, n) = rcw(1)*v(1, 1, n + 1) 185 v(1, 3, n) = rcw(2)*v(1, 1, n + 1) 186 v(1, 4, n) = rcw(3)*v(1, 1, n + 1) 187 END DO 188 189 ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + *** 190 ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + *** 191 ! ** f2*[s||c-2i]{n+1} *** 192 193 DO lc = 2, lc_max 194 195 DO n = 1, nmax - lc 196 197 ! **** Increase the angular momentum component z of c *** 198 199 v(1, coset(0, 0, lc), n) = & 200 rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + & 201 f1*REAL(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + & 202 f2*v(1, coset(0, 0, lc - 2), n + 1)) 203 204 ! *** Increase the angular momentum component y of c *** 205 206 cz = lc - 1 207 v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1) 208 209 DO cy = 2, lc 210 cz = lc - cy 211 v(1, coset(0, cy, cz), n) = & 212 rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + & 213 f1*REAL(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + & 214 f2*v(1, coset(0, cy - 2, cz), n + 1)) 215 END DO 216 217 ! *** Increase the angular momentum component x of c *** 218 219 DO cy = 0, lc - 1 220 cz = lc - 1 - cy 221 v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1) 222 END DO 223 224 DO cx = 2, lc 225 f6 = f1*REAL(cx - 1, dp) 226 DO cy = 0, lc - cx 227 cz = lc - cx - cy 228 v(1, coset(cx, cy, cz), n) = & 229 rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + & 230 f6*(v(1, coset(cx - 2, cy, cz), n) + & 231 f2*v(1, coset(cx - 2, cy, cz), n + 1)) 232 END DO 233 END DO 234 235 END DO 236 237 END DO 238 239 END IF 240 241 ! *** Vertical recurrence steps: [s||c] -> [a||c] *** 242 243 IF (la_max > 0) THEN 244 245 f3 = 0.5_dp*zetp 246 f4 = -rho*zetp 247 f5 = 0.5_dp*zetw 248 249 raw(:) = zetc(jpgf)*zetw*rac(:) 250 251 ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) *** 252 253 DO n = 1, nmax - 1 254 v(2, 1, n) = raw(1)*v(1, 1, n + 1) 255 v(3, 1, n) = raw(2)*v(1, 1, n + 1) 256 v(4, 1, n) = raw(3)*v(1, 1, n + 1) 257 END DO 258 259 ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + *** 260 ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + *** 261 ! *** f4*[a-2i||s]{n+1}) *** 262 263 DO la = 2, la_max 264 265 DO n = 1, nmax - la 266 267 ! *** Increase the angular momentum component z of a *** 268 269 v(coset(0, 0, la), 1, n) = & 270 raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + & 271 f3*REAL(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + & 272 f4*v(coset(0, 0, la - 2), 1, n + 1)) 273 274 ! *** Increase the angular momentum component y of a *** 275 276 az = la - 1 277 v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1) 278 279 DO ay = 2, la 280 az = la - ay 281 v(coset(0, ay, az), 1, n) = & 282 raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + & 283 f3*REAL(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + & 284 f4*v(coset(0, ay - 2, az), 1, n + 1)) 285 END DO 286 287 ! *** Increase the angular momentum component x of a *** 288 289 DO ay = 0, la - 1 290 az = la - 1 - ay 291 v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1) 292 END DO 293 294 DO ax = 2, la 295 f6 = f3*REAL(ax - 1, dp) 296 DO ay = 0, la - ax 297 az = la - ax - ay 298 v(coset(ax, ay, az), 1, n) = & 299 raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + & 300 f6*(v(coset(ax - 2, ay, az), 1, n) + & 301 f4*v(coset(ax - 2, ay, az), 1, n + 1)) 302 END DO 303 END DO 304 305 END DO 306 307 END DO 308 309 DO lc = 1, lc_max 310 311 DO cx = 0, lc 312 DO cy = 0, lc - cx 313 cz = lc - cx - cy 314 315 coc = coset(cx, cy, cz) 316 cocx = coset(MAX(0, cx - 1), cy, cz) 317 cocy = coset(cx, MAX(0, cy - 1), cz) 318 cocz = coset(cx, cy, MAX(0, cz - 1)) 319 320 fcx = f5*REAL(cx, dp) 321 fcy = f5*REAL(cy, dp) 322 fcz = f5*REAL(cz, dp) 323 324 ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + *** 325 ! *** f5*Ni(c)*[s||c-1i]{n+1} *** 326 327 DO n = 1, nmax - 1 - lc 328 v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1) 329 v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1) 330 v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1) 331 END DO 332 333 ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + *** 334 ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + *** 335 ! *** f4*[a-2i||c]{n+1}) + *** 336 ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} *** 337 338 DO la = 2, la_max 339 340 DO n = 1, nmax - la - lc 341 342 ! *** Increase the angular momentum component z of a *** 343 344 v(coset(0, 0, la), coc, n) = & 345 raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + & 346 f3*REAL(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + & 347 f4*v(coset(0, 0, la - 2), coc, n + 1)) + & 348 fcz*v(coset(0, 0, la - 1), cocz, n + 1) 349 350 ! *** Increase the angular momentum component y of a *** 351 352 az = la - 1 353 v(coset(0, 1, az), coc, n) = & 354 raw(2)*v(coset(0, 0, az), coc, n + 1) + & 355 fcy*v(coset(0, 0, az), cocy, n + 1) 356 357 DO ay = 2, la 358 az = la - ay 359 v(coset(0, ay, az), coc, n) = & 360 raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + & 361 f3*REAL(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + & 362 f4*v(coset(0, ay - 2, az), coc, n + 1)) + & 363 fcy*v(coset(0, ay - 1, az), cocy, n + 1) 364 END DO 365 366 ! *** Increase the angular momentum component x of a *** 367 368 DO ay = 0, la - 1 369 az = la - 1 - ay 370 v(coset(1, ay, az), coc, n) = & 371 raw(1)*v(coset(0, ay, az), coc, n + 1) + & 372 fcx*v(coset(0, ay, az), cocx, n + 1) 373 END DO 374 375 DO ax = 2, la 376 f6 = f3*REAL(ax - 1, dp) 377 DO ay = 0, la - ax 378 az = la - ax - ay 379 v(coset(ax, ay, az), coc, n) = & 380 raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + & 381 f6*(v(coset(ax - 2, ay, az), coc, n) + & 382 f4*v(coset(ax - 2, ay, az), coc, n + 1)) + & 383 fcx*v(coset(ax - 1, ay, az), cocx, n + 1) 384 END DO 385 END DO 386 387 END DO 388 389 END DO 390 391 END DO 392 END DO 393 394 END DO 395 396 END IF 397 398 DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local) 399 DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local) 400 vac(na + i, nc + j) = v(i, j, 1) 401 END DO 402 END DO 403 404 IF (PRESENT(maxder)) THEN 405 DO j = 1, ncoset(lc_max) 406 DO i = 1, ncoset(la_max) 407 vac_plus(nap + i, ncp + j) = v(i, j, 1) 408 END DO 409 END DO 410 END IF 411 412 nc = nc + ncoset(lc_max - maxder_local) 413 ncp = ncp + ncoset(lc_max) 414 415 END DO 416 417 na = na + ncoset(la_max - maxder_local) 418 nap = nap + ncoset(la_max) 419 420 END DO 421 422 CALL timestop(handle) 423 424 END SUBROUTINE operator2 425 426! ************************************************************************************************** 427!> \brief Calculation of Coulomb integrals for s-function, i.e, [s|1/r12|s], and the auxiliary 428!> integrals [s|1/r12|s]^n 429!> \param v matrix storing the integrals 430!> \param nmax maximal n in the auxiliary integrals [s|1/r12|s]^n 431!> \param zetp = 1/zeta 432!> \param zetq = 1/zetc 433!> \param zetw = 1/(zeta+zetc) 434!> \param rho = zeta*zetc*zetw 435!> \param rac2 square distance between center A and C, |Ra-Rc|^2 436!> \param omega this parameter is actually not used, but included for the sake of the abstract 437!> interface 438!> \param r_cutoff same as above 439! ************************************************************************************************** 440 SUBROUTINE cps_coulomb2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 441 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 442 INTEGER, INTENT(IN) :: nmax 443 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, & 444 r_cutoff 445 446 CHARACTER(len=*), PARAMETER :: routineN = 'cps_coulomb2', routineP = moduleN//':'//routineN 447 448 INTEGER :: n 449 REAL(KIND=dp) :: f0, t 450 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 451 452 MARK_USED(omega) 453 MARK_USED(r_cutoff) 454 455 ALLOCATE (f(0:nmax)) 456 f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq 457 458 ! *** Calculate the incomplete Gamma/Boys function *** 459 460 t = rho*rac2 461 CALL fgamma(nmax - 1, t, f) 462 463 ! *** Calculate the basic two-center integrals [s||s]{n} *** 464 465 DO n = 1, nmax 466 v(1, 1, n) = f0*f(n - 1) 467 END DO 468 469 DEALLOCATE (f) 470 END SUBROUTINE cps_coulomb2 471 472! ************************************************************************************************** 473!> \brief Calculation of verf integrals for s-function, i.e, [s|erf(omega*r12)/r12|s], and the 474!> auxiliary integrals [s|erf(omega*r12)/r12|s]^n 475!> \param v matrix storing the integrals 476!> \param nmax maximal n in the auxiliary integrals [s|erf(omega*r12)/r12|s]^n 477!> \param zetp = 1/zeta 478!> \param zetq = 1/zetc 479!> \param zetw = 1/(zeta+zetc) 480!> \param rho = zeta*zetc*zetw 481!> \param rac2 square distance between center A and C, |Ra-Rc|^2 482!> \param omega parameter in the operator 483!> \param r_cutoff dummy argument for the sake of generality 484! ************************************************************************************************** 485 SUBROUTINE cps_verf2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 486 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 487 INTEGER, INTENT(IN) :: nmax 488 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, & 489 r_cutoff 490 491 CHARACTER(len=*), PARAMETER :: routineN = 'cps_verf2', routineP = moduleN//':'//routineN 492 493 INTEGER :: n 494 REAL(KIND=dp) :: arg, comega, f0, t 495 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 496 497 MARK_USED(r_cutoff) 498 499 ALLOCATE (f(0:nmax)) 500 comega = omega**2/(omega**2 + rho) 501 f0 = 2.0_dp*SQRT(pi**5*zetw*comega)*zetp*zetq 502 503 ! *** Calculate the incomplete Gamma/Boys function *** 504 505 t = rho*rac2 506 arg = comega*t 507 CALL fgamma(nmax - 1, arg, f) 508 509 ! *** Calculate the basic two-center integrals [s||s]{n} *** 510 511 DO n = 1, nmax 512 v(1, 1, n) = f0*f(n - 1)*comega**(n - 1) 513 END DO 514 515 DEALLOCATE (f) 516 517 END SUBROUTINE cps_verf2 518 519! ************************************************************************************************** 520!> \brief Calculation of verfc integrals for s-function, i.e, [s|erfc(omega*r12)/r12|s], and 521!> the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n 522!> \param v matrix storing the integrals 523!> \param nmax maximal n in the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n 524!> \param zetp = 1/zeta 525!> \param zetq = 1/zetc 526!> \param zetw = 1/(zeta+zetc) 527!> \param rho = zeta*zetc*zetw 528!> \param rac2 square distance between center A and C, |Ra-Rc|^2 529!> \param omega parameter in the operator 530!> \param r_cutoff dummy argument for the sake of generality 531! ************************************************************************************************** 532 SUBROUTINE cps_verfc2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 533 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 534 INTEGER, INTENT(IN) :: nmax 535 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, & 536 r_cutoff 537 538 CHARACTER(len=*), PARAMETER :: routineN = 'cps_verfc2', routineP = moduleN//':'//routineN 539 540 INTEGER :: n 541 REAL(KIND=dp) :: argerf, comega, f0, t 542 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf 543 544 MARK_USED(r_cutoff) 545 546 ALLOCATE (fv(0:nmax), fverf(0:nmax)) 547 comega = omega**2/(omega**2 + rho) 548 f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq 549 550 ! *** Calculate the incomplete Gamma/Boys function *** 551 552 t = rho*rac2 553 argerf = comega*t 554 555 CALL fgamma(nmax - 1, t, fv) 556 CALL fgamma(nmax - 1, argerf, fverf) 557 558 ! *** Calculate the basic two-center integrals [s||s]{n} *** 559 560 DO n = 1, nmax 561 v(1, 1, n) = f0*(fv(n - 1) - SQRT(comega)*comega**(n - 1)*fverf(n - 1)) 562 END DO 563 564 DEALLOCATE (fv, fverf) 565 566 END SUBROUTINE cps_verfc2 567 568! ************************************************************************************************** 569!> \brief Calculation of vgauss integrals for s-function, i.e, [s|exp(-omega*r12^2)/r12|s], and 570!> the auxiliary integrals [s|exp(-omega*r12^2)/r12|s] 571!> \param v matrix storing the integrals 572!> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)/r12|s] 573!> \param zetp = 1/zeta 574!> \param zetq = 1/zetc 575!> \param zetw = 1/(zeta+zetc) 576!> \param rho = zeta*zetc*zetw 577!> \param rac2 square distance between center A and C, |Ra-Rc|^2 578!> \param omega parameter in the operator 579!> \param r_cutoff dummy argument for the sake of generality 580! ************************************************************************************************** 581 SUBROUTINE cps_vgauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 582 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 583 INTEGER, INTENT(IN) :: nmax 584 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, & 585 r_cutoff 586 587 CHARACTER(len=*), PARAMETER :: routineN = 'cps_vgauss2', routineP = moduleN//':'//routineN 588 589 INTEGER :: j, n 590 REAL(KIND=dp) :: arg, dummy, eta, expT, f0, fsign, t, tau 591 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 592 593 MARK_USED(r_cutoff) 594 595 ALLOCATE (f(0:nmax)) 596 597 dummy = zetp 598 dummy = zetq 599 eta = rho/(rho + omega) 600 tau = omega/(rho + omega) 601 602 ! *** Calculate the incomplete Gamma/Boys function *** 603 604 t = rho*rac2 605 arg = eta*t 606 607 CALL fgamma(nmax - 1, arg, f) 608 609 expT = EXP(-omega/(omega + rho)*t) 610 f0 = 2.0_dp*SQRT(pi**5*zetw**3)/(rho + omega)*expT 611 612 ! *** Calculate the basic two-center integrals [s||s]{n} *** 613 v(1, 1, 1:nmax) = 0.0_dp 614 DO n = 1, nmax 615 fsign = (-1.0_dp)**(n - 1) 616 DO j = 0, n - 1 617 v(1, 1, n) = v(1, 1, n) + f0*fsign* & 618 fac(n - 1)/fac(n - j - 1)/fac(j)*(-tau)**(n - j - 1)*(-eta)**j*f(j) 619 ENDDO 620 ENDDO 621 622 DEALLOCATE (f) 623 624 END SUBROUTINE cps_vgauss2 625 626! ************************************************************************************************** 627!> \brief Calculation of gauss integrals for s-function, i.e, [s|exp(-omega*r12^2)|s], and 628!> the auxiliary integrals [s|exp(-omega*r12^2)|s] 629!> \param v matrix storing the integrals 630!> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)|s] 631!> \param zetp = 1/zeta 632!> \param zetq = 1/zetc 633!> \param zetw = 1/(zeta+zetc) 634!> \param rho = zeta*zetc*zetw 635!> \param rac2 square distance between center A and C, |Ra-Rc|^2 636!> \param omega parameter in the operator 637!> \param r_cutoff dummy argument for the sake of generality 638! ************************************************************************************************** 639 SUBROUTINE cps_gauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 640 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 641 INTEGER, INTENT(IN) :: nmax 642 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, & 643 r_cutoff 644 645 CHARACTER(len=*), PARAMETER :: routineN = 'cps_gauss2', routineP = moduleN//':'//routineN 646 647 INTEGER :: n 648 REAL(KIND=dp) :: dummy, expT, f0, t, tau 649 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 650 651 MARK_USED(r_cutoff) 652 653 ALLOCATE (f(0:nmax)) 654 655 dummy = zetp 656 dummy = zetq 657 tau = omega/(rho + omega) 658 t = rho*rac2 659 expT = EXP(-tau*t) 660 f0 = pi**3*SQRT(zetw**3/(rho + omega)**3)*expT 661 662 ! *** Calculate the basic two-center integrals [s||s]{n} *** 663 664 DO n = 1, nmax 665 v(1, 1, n) = f0*tau**(n - 1) 666 END DO 667 668 DEALLOCATE (f) 669 670 END SUBROUTINE cps_gauss2 671 672! ************************************************************************************************** 673!> \brief Calculation of truncated Coulomb integrals for s-function, i.e, [s|TC|s] where TC = 1/r12 674!> if r12 <= r_cutoff and 0 otherwise 675!> \param v matrix storing the integrals 676!> \param nmax maximal n in the auxiliary integrals [s|TC|s] 677!> \param zetp = 1/zeta 678!> \param zetq = 1/zetc 679!> \param zetw = 1/(zeta+zetc) 680!> \param rho = zeta*zetc*zetw 681!> \param rac2 square distance between center A and C, |Ra-Rc|^2 682!> \param omega dummy argument for the sake of generality 683!> \param r_cutoff the radius at which the operator is cut 684!> \note The truncated operator must have been initialized from the data file prior to this call 685! ************************************************************************************************** 686 SUBROUTINE cps_truncated2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff) 687 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v 688 INTEGER, INTENT(IN) :: nmax 689 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, & 690 r_cutoff 691 692 CHARACTER(len=*), PARAMETER :: routineN = 'cps_truncated2', routineP = moduleN//':'//routineN 693 694 INTEGER :: n 695 LOGICAL :: use_gamma 696 REAL(KIND=dp) :: f0, r, t 697 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f 698 699 MARK_USED(omega) 700 701 ALLOCATE (f(nmax + 1)) !t_c_g0 needs to start at index 1 702 703 r = r_cutoff*SQRT(rho) 704 t = rho*rac2 705 f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq 706 707 !check that the operator has been init from file 708 CPASSERT(get_lmax_init() .GE. nmax) 709 710 CALL t_c_g0_n(f, use_gamma, r, t, nmax) 711 IF (use_gamma) CALL fgamma(nmax, t, f) 712 713 DO n = 1, nmax 714 v(1, 1, n) = f0*f(n) 715 END DO 716 717 DEALLOCATE (f) 718 719 END SUBROUTINE cps_truncated2 720 721END MODULE ai_operators_r12 722