1c 2c $Id$ 3c 4 5* ****************************************************** 6* * * 7* * nwpw_WGaussian * 8* * * 9* ****************************************************** 10* 11* Calculates the two electron two center Gaussian integral 12* 13* // 14* WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 15* || ------------------------------------ dr dr' 16* // |r-r'| 17* 18* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat) 19* 20* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 21* 22* The normalization constant C_l is defined such at 23* / 24* | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1 25* / 26* 27 real*8 function nwpw_WGaussian(la,ma,sa,lb,mb,sb,Rab) 28 implicit none 29 integer la,ma,lb,mb 30 real*8 sa,sb,Rab(3) 31 32* *** local variables *** 33 integer l,m,fac,fac2 34 real*8 c,x,y,pi,tmp,mtmp,alpha 35 real*8 cos_theta,phi,R 36 37* **** external functions **** 38 integer nwpw_doublefactorial 39 external nwpw_doublefactorial 40 real*8 nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm 41 external nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm 42 43 call nwpw_timing_start(49) 44 pi = 4.0d0*datan(1.0d0) 45 x = dble(nwpw_doublefactorial(2*la+1)) 46 y = dble(nwpw_doublefactorial(2*lb+1)) 47 alpha = dsqrt(0.25d0*(sa*sa + sb*sb)) 48 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 49 cos_theta = Rab(3)/R 50 51 if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then 52 phi = 0.0d0 53 else 54 phi = datan2(Rab(2),Rab(1)) 55 end if 56 57 if (mod(2*la+lb,2).eq.1) then 58 c = -32.0d0*pi/(x*y) 59 else 60 c = 32.0d0*pi/(x*y) 61 end if 62 63 if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 64 fac = -1 65 else 66 fac = 1 67 end if 68 69 tmp = 0.0d0 70 do l = abs(la-lb), (la+lb), 2 71 mtmp = 0.0d0 72 do m=-l,l 73 mtmp = mtmp + nwpw_gaunt(.false.,l,m,la,ma,lb,mb) 74 > *Tesseral_lm(l,m,cos_theta,phi) 75 end do 76 tmp = tmp + fac * mtmp * nwpw_GaussBessel(la+lb,l,alpha,R) 77 fac = -fac 78 end do 79 call nwpw_timing_end(49) 80 81 nwpw_WGaussian = c * tmp 82 return 83 end 84 85 86* ****************************************************** 87* * * 88* * nwpw_dWGaussian * 89* * * 90* ****************************************************** 91* 92* Calculates the two electron two center Gaussian integral and it's derivative wrt to Rab 93* 94* // 95* dWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 96* || ------------------------------------ dr dr' 97* // |r-r'| 98* 99* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat) 100* 101* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 102* 103* The normalization constant C_l is defined such at 104* / 105* | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1 106* / 107* 108 subroutine nwpw_dWGaussian(la,ma,sa,lb,mb,sb,Rab,W,dW) 109 implicit none 110 integer la,ma,lb,mb 111 real*8 sa,sb,Rab(3) 112 real*8 W,dW(3) 113 114* *** local variables *** 115 integer l,m,fac 116 real*8 c,x,y,pi,tmp,mtmp,mtmpx,mtmpy,mtmpz,alpha 117 real*8 cos_theta,phi,R,gg1,gg2,gg3,Tx,Ty,Tz 118 119* **** external functions **** 120 integer nwpw_doublefactorial 121 external nwpw_doublefactorial 122 real*8 nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm 123 external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm 124 125 call nwpw_timing_start(55) 126 pi = 4.0d0*datan(1.0d0) 127 x = dble(nwpw_doublefactorial(2*la+1)) 128 y = dble(nwpw_doublefactorial(2*lb+1)) 129 alpha = dsqrt(0.25d0*(sa*sa + sb*sb)) 130 131 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 132 cos_theta = Rab(3)/R 133 phi = datan2(Rab(2),Rab(1)) 134 135 if (mod(2*la+lb,2).eq.1) then 136 c = -32.0d0*pi/(x*y) 137 else 138 c = 32.0d0*pi/(x*y) 139 end if 140 141 if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 142 fac = -1 143 else 144 fac = 1 145 end if 146 147 W = 0.0d0 148 dW(1) = 0.0d0 149 dW(2) = 0.0d0 150 dW(3) = 0.0d0 151 do l = abs(la-lb), (la+lb), 2 152 mtmp = 0.0d0 153 mtmpx = 0.0d0 154 mtmpy = 0.0d0 155 mtmpz = 0.0d0 156 do m=-l,l 157 gg1 = nwpw_gaunt(.false.,l,m,la,ma,lb,mb) 158 call dTesseral_lm(l,m,cos_theta,phi,Tx,Ty,Tz) 159 mtmp = mtmp + gg1*Tesseral_lm(l,m,cos_theta,phi) 160 mtmpx = mtmpx + gg1*Tx 161 mtmpy = mtmpy + gg1*Ty 162 mtmpz = mtmpz + gg1*Tz 163 end do 164 gg2 = nwpw_GaussBessel(la+lb,l,alpha,R) 165 gg3 = nwpw_dGaussBessel(la+lb,l,alpha,R) 166 W = W + fac*(mtmp*gg2) 167 dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3) 168 dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3) 169 dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3) 170 fac = -fac 171 end do 172 W = W*c 173 dW(1) = dW(1)*c 174 dW(2) = dW(2)*c 175 dW(3) = dW(3)*c 176 call nwpw_timing_end(55) 177 178 return 179 end 180 181* ****************************************************** 182* * * 183* * nwpw_UGaussian * 184* * * 185* ****************************************************** 186* 187* Calculates the two electron one center Gaussian integral 188* 189* // 190* WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r) * g(lb,mb,sb;r') 191* || ------------------------------------ dr dr' 192* // |r-r'| 193* 194* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat) 195* 196* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 197* 198* The normalization constant C_l is defined such at 199* / 200* | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1 201* / 202* 203* 204* Note - this routine is equivalent to find_self_energy_coeff in the paw code 205* 206 real*8 function nwpw_UGaussian(la,ma,sa,lb,mb,sb) 207 implicit none 208 integer la,ma,lb,mb 209 real*8 sa,sb 210 211* *** local variables *** 212 real*8 x,y,twopi,s,U 213 214* **** external functions **** 215 integer nwpw_doublefactorial 216 external nwpw_doublefactorial 217 218 call nwpw_timing_start(48) 219 U = 0.0d0 220 if ((la.eq.lb).and.(ma.eq.mb)) then 221 twopi = 8.0d0*datan(1.0d0) 222 x = dble((2*la+1)*nwpw_doublefactorial(2*la+1)) 223 y = (dsqrt(0.5d0*(sa*sa+sb*sb)))**(2*la+1) 224 U = 4.0d0*dsqrt(twopi)/(x*y) 225 end if 226 call nwpw_timing_end(48) 227 228 nwpw_UGaussian = U 229 return 230 end 231 232 233*************************** real combo versions *********************************** 234 235* ****************************************************** 236* * * 237* * nwpw_WGaussian3 * 238* * * 239* ****************************************************** 240* 241* Calculates the 4 terms sum 242* 243* WGaussian(la,ma,sa,lb,mb,sb,Rab) 244* - WGaussian(la,ma,sa,lb,mb,sm,Rab) 245* - WGaussian(la,ma,sm,lb,mb,sb,Rab) 246* + WGaussian(la,ma,sm,lb,mb,sm,Rab) 247* 248* of the two electron two center Gaussian integral 249* 250* // 251* WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 252* || ------------------------------------ dr dr' 253* // |r-r'| 254* 255* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat) 256* 257* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 258* 259* The normalization constant C_l is defined such at 260* / 261* | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1 262* / 263* 264 real*8 function nwpw_WGaussian3(la,ma,sa,lb,mb,sb,sm,Rab) 265 implicit none 266 integer la,ma,lb,mb 267 real*8 sa,sb,sm,Rab(3) 268 269* *** local variables *** 270 integer l,m,fac,fac2 271 real*8 c,x,y,pi,mtmp 272 real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm 273 real*8 tmp 274 real*8 cos_theta,phi,R,gg2 275 276* **** external functions **** 277 integer nwpw_doublefactorial 278 external nwpw_doublefactorial 279 real*8 nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm 280 external nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm 281 282 call nwpw_timing_start(49) 283 pi = 4.0d0*datan(1.0d0) 284 x = dble(nwpw_doublefactorial(2*la+1)) 285 y = dble(nwpw_doublefactorial(2*lb+1)) 286 alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb)) 287 alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm)) 288 alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb)) 289 alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm)) 290 291 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 292 cos_theta = Rab(3)/R 293 294 if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then 295 phi = 0.0d0 296 else 297 phi = datan2(Rab(2),Rab(1)) 298 end if 299 300 if (mod(2*la+lb,2).eq.1) then 301 c = -32.0d0*pi/(x*y) 302 else 303 c = 32.0d0*pi/(x*y) 304 end if 305 306 if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 307 fac = -1 308 else 309 fac = 1 310 end if 311 312 tmp = 0.0d0 313 do l = abs(la-lb), (la+lb), 2 314 mtmp = 0.0d0 315 do m=-l,l 316 mtmp = mtmp + nwpw_gaunt(.false.,l,m,la,ma,lb,mb) 317 > *Tesseral_lm(l,m,cos_theta,phi) 318 end do 319 gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R) 320 > - nwpw_GaussBessel(la+lb,l,alpha_am,R) 321 > - nwpw_GaussBessel(la+lb,l,alpha_mb,R) 322 > + nwpw_GaussBessel(la+lb,l,alpha_mm,R) 323 324 tmp = tmp + fac*mtmp*gg2 325 fac = -fac 326 end do 327 call nwpw_timing_end(49) 328 329 nwpw_WGaussian3 = c*tmp 330 return 331 end 332 333* ****************************************************** 334* * * 335* * nwpw_WGaussian2 * 336* * * 337* ****************************************************** 338* 339* Calculates the 4 term sum 340* 341* WGaussian(la,ma,sa,lb,mb,sa,Rab) 342* - WGaussian(la,ma,sa,lb,mb,sm,Rab) 343* - WGaussian(la,ma,sm,lb,mb,sa,Rab) 344* + WGaussian(la,ma,sm,lb,mb,sm,Rab) 345* 346* of the two electron two center Gaussian integral 347* 348* // 349* WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 350* || ------------------------------------ dr dr' 351* // |r-r'| 352* 353* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat) 354* 355* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 356* 357* The normalization constant C_l is defined such at 358* / 359* | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1 360* / 361* 362 real*8 function nwpw_WGaussian2(la,ma,sa,lb,mb,sm,Rab) 363 implicit none 364 integer la,ma,lb,mb 365 real*8 sa,sb,sm,Rab(3) 366 367* *** local variables *** 368 integer l,m,fac,fac2 369 real*8 c,x,y,pi,mtmp 370 real*8 alpha_aa,alpha_am,alpha_mm 371 real*8 tmp 372 real*8 cos_theta,phi,R,gg2 373 374* **** external functions **** 375 integer nwpw_doublefactorial 376 external nwpw_doublefactorial 377 real*8 nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm 378 external nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm 379 380 call nwpw_timing_start(49) 381 pi = 4.0d0*datan(1.0d0) 382 x = dble(nwpw_doublefactorial(2*la+1)) 383 y = dble(nwpw_doublefactorial(2*lb+1)) 384 alpha_aa = dsqrt(0.25d0*(sa*sa + sa*sa)) 385 alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm)) 386 alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm)) 387 388 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 389 cos_theta = Rab(3)/R 390 391 if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then 392 phi = 0.0d0 393 else 394 phi = datan2(Rab(2),Rab(1)) 395 end if 396 397 if (mod(2*la+lb,2).eq.1) then 398 c = -32.0d0*pi/(x*y) 399 else 400 c = 32.0d0*pi/(x*y) 401 end if 402 403 if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 404 fac = -1 405 else 406 fac = 1 407 end if 408 409 tmp = 0.0d0 410 do l = abs(la-lb), (la+lb), 2 411 mtmp = 0.0d0 412 do m=-l,l 413 mtmp = mtmp + nwpw_gaunt(.false.,l,m,la,ma,lb,mb) 414 > *Tesseral_lm(l,m,cos_theta,phi) 415 end do 416 417 gg2 = nwpw_GaussBessel(la+lb,l,alpha_aa,R) 418 > + nwpw_GaussBessel(la+lb,l,alpha_mm,R) 419 > - 2.0d0*nwpw_GaussBessel(la+lb,l,alpha_am,R) 420 421 tmp = tmp + fac*mtmp*gg2 422 423 fac = -fac 424 end do 425 call nwpw_timing_end(49) 426 427 nwpw_WGaussian2 = c*tmp 428 return 429 end 430 431* ****************************************************** 432* * * 433* * nwpw_dWGaussian3 * 434* * * 435* ****************************************************** 436* 437* Calculates the 4 terms sum 438* 439* WGaussian(la,ma,sa,lb,mb,sb,Rab) 440* - WGaussian(la,ma,sa,lb,mb,sm,Rab) 441* - WGaussian(la,ma,sm,lb,mb,sb,Rab) 442* + WGaussian(la,ma,sm,lb,mb,sm,Rab) 443* 444* of the two electron two center Gaussian integral and it's derivative wrt to Rab 445* 446* // 447* dWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 448* || ------------------------------------ dr dr' 449* // |r-r'| 450* 451* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat) 452* 453* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 454* 455* The normalization constant C_l is defined such at 456* / 457* | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1 458* / 459* 460 subroutine nwpw_dWGaussian3(la,ma,sa,lb,mb,sb,sm,Rab,W,dW) 461 implicit none 462 integer la,ma,lb,mb 463 real*8 sa,sb,sm,Rab(3) 464 real*8 W,dW(3) 465 466* *** local variables *** 467 integer l,m,fac 468 real*8 c,x,y,pi,tmp,mtmp,mtmpx,mtmpy,mtmpz,alpha 469 real*8 cos_theta,phi,R,gg1,gg2,gg3,Tx,Ty,Tz 470 real*8 W_ab,W_am,W_mb,W_mm 471 real*8 dW_ab(3),dW_am(3),dW_mb(3),dW_mm(3) 472 real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm 473 474* **** external functions **** 475 integer nwpw_doublefactorial 476 external nwpw_doublefactorial 477 real*8 nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm 478 external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm 479 480 call nwpw_timing_start(55) 481 pi = 4.0d0*datan(1.0d0) 482 x = dble(nwpw_doublefactorial(2*la+1)) 483 y = dble(nwpw_doublefactorial(2*lb+1)) 484 alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb)) 485 alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm)) 486 alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb)) 487 alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm)) 488 489 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 490 cos_theta = Rab(3)/R 491 phi = datan2(Rab(2),Rab(1)) 492 493 if (mod(2*la+lb,2).eq.1) then 494 c = -32.0d0*pi/(x*y) 495 else 496 c = 32.0d0*pi/(x*y) 497 end if 498 499 if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 500 fac = -1 501 else 502 fac = 1 503 end if 504 505 W = 0.0d0 506 dW(1) = 0.0d0 507 dW(2) = 0.0d0 508 dW(3) = 0.0d0 509 510 do l = abs(la-lb), (la+lb), 2 511 mtmp = 0.0d0 512 mtmpx = 0.0d0 513 mtmpy = 0.0d0 514 mtmpz = 0.0d0 515 do m=-l,l 516 gg1 = nwpw_gaunt(.false.,l,m,la,ma,lb,mb) 517 call dTesseral_lm(l,m,cos_theta,phi,Tx,Ty,Tz) 518 mtmp = mtmp + gg1*Tesseral_lm(l,m,cos_theta,phi) 519 mtmpx = mtmpx + gg1*Tx 520 mtmpy = mtmpy + gg1*Ty 521 mtmpz = mtmpz + gg1*Tz 522 end do 523 524 gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R) 525 > - nwpw_GaussBessel(la+lb,l,alpha_am,R) 526 > - nwpw_GaussBessel(la+lb,l,alpha_mb,R) 527 > + nwpw_GaussBessel(la+lb,l,alpha_mm,R) 528 529 gg3 = nwpw_dGaussBessel(la+lb,l,alpha_ab,R) 530 > - nwpw_dGaussBessel(la+lb,l,alpha_am,R) 531 > - nwpw_dGaussBessel(la+lb,l,alpha_mb,R) 532 > + nwpw_dGaussBessel(la+lb,l,alpha_mm,R) 533 534 W = W + fac*(mtmp*gg2) 535 dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3) 536 dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3) 537 dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3) 538 539 fac = -fac 540 end do 541 W = W*c 542 dW(1) = dW(1)*c 543 dW(2) = dW(2)*c 544 dW(3) = dW(3)*c 545 call nwpw_timing_end(55) 546 547 return 548 end 549 550 551*************************** real combo versions *********************************** 552 553 554 555*************************** complex versions *********************************** 556 557 558* ****************************************************** 559* * * 560* * nwpw_CWGaussian * 561* * * 562* ****************************************************** 563* 564* Calculates the two electron two center Gaussian integral 565* 566* // 567* CWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 568* || ------------------------------------ dr dr' 569* // |r-r'| 570* 571* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat) 572* 573* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 574* 575* The normalization constant C_l is defined such at 576* / 577* | g(l,m,s;r) * |r|**l *conjg(Ylm(rhat)) dr = 1 578* / 579* 580 complex*16 function nwpw_CWGaussian(la,ma,sa,lb,mb,sb,Rab) 581 implicit none 582 integer la,ma,lb,mb 583 real*8 sa,sb,Rab(3) 584 585* *** local variables *** 586 integer l,m,fac,fac2 587 real*8 c,x,y,pi,alpha 588 real*8 cos_theta,phi,R 589 complex*16 mtmp,tmp 590 591* **** external functions **** 592 integer nwpw_doublefactorial 593 external nwpw_doublefactorial 594 real*8 nwpw_gaunt,nwpw_GaussBessel 595 external nwpw_gaunt,nwpw_GaussBessel 596 complex*16 Yspherical_lm 597 external Yspherical_lm 598 real*8 nwpw_gaunt_sub,gen_gaunt_coeff_sub 599 external nwpw_gaunt_sub,gen_gaunt_coeff_sub 600 601 602 call nwpw_timing_start(49) 603 pi = 4.0d0*datan(1.0d0) 604 x = dble(nwpw_doublefactorial(2*la+1)) 605 y = dble(nwpw_doublefactorial(2*lb+1)) 606 alpha = dsqrt(0.25d0*(sa*sa + sb*sb)) 607 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 608 cos_theta = Rab(3)/R 609 610 if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then 611 phi = 0.0d0 612 else 613 phi = datan2(Rab(2),Rab(1)) 614 end if 615 616 if (mod(2*la+lb,2).eq.1) then 617 c = -32.0d0*pi/(x*y) 618 else 619 c = 32.0d0*pi/(x*y) 620 end if 621 622 !phase_factor = (-1)**(m1+l1)/sigma**(l1+l2+1) 623 !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 624 !if (mod(lb+mb,2).eq.1) then 625 626 627 !if (mod(lb+mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 628 if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 629 fac = -1 630 else 631 fac = 1 632 end if 633 634 m = ma + mb 635 tmp = 0.0d0 636 do l = abs(la-lb), (la+lb), 2 637 if (abs(m).le.l) then 638c write(*,*) "l,m=",l,m,la,ma,lb,mb, 639c > nwpw_gaunt(.true.,l,m,la,ma,lb,-mb), 640c > nwpw_gaunt_sub(.true.,l,m,la,ma,lb,-mb), 641c > gen_gaunt_coeff_sub(l,m,la,ma,lb,-mb) 642 mtmp = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb) 643 > *YSpherical_lm(l,m,cos_theta,phi) 644 tmp = tmp + fac * mtmp * nwpw_GaussBessel(la+lb,l,alpha,R) 645 end if 646 fac = -fac 647 end do 648 call nwpw_timing_end(49) 649 650 nwpw_CWGaussian = c * tmp 651 return 652 end 653 654 655* ****************************************************** 656* * * 657* * nwpw_dCWGaussian * 658* * * 659* ****************************************************** 660* 661* Calculates the two electron two center Gaussian integral and it's derivative wrt to Rab 662* 663* // 664* dCWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 665* || ------------------------------------ dr dr' 666* // |r-r'| 667* 668* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat) 669* 670* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 671* 672* The normalization constant C_l is defined such at 673* / 674* | g(l,m,s;r) * |r|**l * congj(Ylm(rhat)) dr = 1 675* / 676* 677 subroutine nwpw_dCWGaussian(la,ma,sa,lb,mb,sb,Rab,W,dW) 678 implicit none 679 integer la,ma,lb,mb 680 real*8 sa,sb,Rab(3) 681 complex*16 W,dW(3) 682 683* *** local variables *** 684 integer l,m,fac 685 real*8 c,x,y,pi,alpha 686 real*8 cos_theta,phi,R,gg1,gg2,gg3 687 complex*16 mtmp,mtmpx,mtmpy,mtmpz,Tx,Ty,Tz 688 689* **** external functions **** 690 integer nwpw_doublefactorial 691 external nwpw_doublefactorial 692 real*8 nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel 693 external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel 694 complex*16 Yspherical_lm 695 external Yspherical_lm 696 697 call nwpw_timing_start(55) 698 pi = 4.0d0*datan(1.0d0) 699 x = dble(nwpw_doublefactorial(2*la+1)) 700 y = dble(nwpw_doublefactorial(2*lb+1)) 701 alpha = dsqrt(0.25d0*(sa*sa + sb*sb)) 702 703 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 704 cos_theta = Rab(3)/R 705 phi = datan2(Rab(2),Rab(1)) 706 707 if (mod(2*la+lb,2).eq.1) then 708 c = -32.0d0*pi/(x*y) 709 else 710 c = 32.0d0*pi/(x*y) 711 end if 712 713 !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 714 !if (mod(lb+mb,2).eq.1) then 715 716 if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 717 fac = -1 718 else 719 fac = 1 720 end if 721 722 W = dcmplx(0.0d0,0.0d0) 723 dW(1) = dcmplx(0.0d0,0.0d0) 724 dW(2) = dcmplx(0.0d0,0.0d0) 725 dW(3) = dcmplx(0.0d0,0.0d0) 726 m = ma + mb 727 do l = abs(la-lb), (la+lb), 2 728 if (abs(m).le.l) then 729 gg1 = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb) 730 call dYspherical_lm(l,m,cos_theta,phi,Tx,Ty,Tz) 731 mtmp = gg1*Yspherical_lm(l,m,cos_theta,phi) 732 mtmpx = gg1*Tx 733 mtmpy = gg1*Ty 734 mtmpz = gg1*Tz 735 gg2 = nwpw_GaussBessel(la+lb,l,alpha,R) 736 gg3 = nwpw_dGaussBessel(la+lb,l,alpha,R) 737 W = W + fac*(mtmp*gg2) 738 dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3) 739 dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3) 740 dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3) 741 end if 742 fac = -fac 743 end do 744 W = W*c 745 dW(1) = dW(1)*c 746 dW(2) = dW(2)*c 747 dW(3) = dW(3)*c 748 call nwpw_timing_end(55) 749 750 return 751 end 752 753 754*************************** complex versions *********************************** 755 756*************************** complex combo versions *********************************** 757 758* ****************************************************** 759* * * 760* * nwpw_CWGaussian3 * 761* * * 762* ****************************************************** 763* 764* Calculates the two electron two center Gaussian integral 765* 766* // 767* CWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 768* || ------------------------------------ dr dr' 769* // |r-r'| 770* 771* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat) 772* 773* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 774* 775* The normalization constant C_l is defined such at 776* / 777* | g(l,m,s;r) * |r|**l *conjg(Ylm(rhat)) dr = 1 778* / 779* 780 complex*16 function nwpw_CWGaussian3(la,ma,sa,lb,mb,sb,sm,Rab) 781 implicit none 782 integer la,ma,lb,mb 783 real*8 sa,sb,sm,Rab(3) 784 785* *** local variables *** 786 integer l,m,fac,fac2 787 real*8 c,x,y,pi 788 real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm 789 real*8 cos_theta,phi,R,gg2 790 complex*16 mtmp,tmp 791 792* **** external functions **** 793 integer nwpw_doublefactorial 794 external nwpw_doublefactorial 795 real*8 nwpw_gaunt,nwpw_GaussBessel 796 external nwpw_gaunt,nwpw_GaussBessel 797 complex*16 Yspherical_lm 798 external Yspherical_lm 799 real*8 nwpw_gaunt_sub,gen_gaunt_coeff_sub 800 external nwpw_gaunt_sub,gen_gaunt_coeff_sub 801 802 803 call nwpw_timing_start(49) 804 pi = 4.0d0*datan(1.0d0) 805 x = dble(nwpw_doublefactorial(2*la+1)) 806 y = dble(nwpw_doublefactorial(2*lb+1)) 807 alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb)) 808 alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm)) 809 alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb)) 810 alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm)) 811 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 812 cos_theta = Rab(3)/R 813 814 if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then 815 phi = 0.0d0 816 else 817 phi = datan2(Rab(2),Rab(1)) 818 end if 819 820 if (mod(2*la+lb,2).eq.1) then 821 c = -32.0d0*pi/(x*y) 822 else 823 c = 32.0d0*pi/(x*y) 824 end if 825 826 !phase_factor = (-1)**(m1+l1)/sigma**(l1+l2+1) 827 !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 828 !if (mod(lb+mb,2).eq.1) then 829 830 831 !if (mod(lb+mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 832 if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 833 fac = -1 834 else 835 fac = 1 836 end if 837 838 m = ma + mb 839 tmp = dcmplx(0.0d0,0.0d0) 840 do l = abs(la-lb), (la+lb), 2 841 if (abs(m).le.l) then 842c write(*,*) "l,m=",l,m,la,ma,lb,mb, 843c > nwpw_gaunt(.true.,l,m,la,ma,lb,-mb), 844c > nwpw_gaunt_sub(.true.,l,m,la,ma,lb,-mb), 845c > gen_gaunt_coeff_sub(l,m,la,ma,lb,-mb) 846 mtmp = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb) 847 > *YSpherical_lm(l,m,cos_theta,phi) 848 849 gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R) 850 > - nwpw_GaussBessel(la+lb,l,alpha_am,R) 851 > - nwpw_GaussBessel(la+lb,l,alpha_mb,R) 852 > + nwpw_GaussBessel(la+lb,l,alpha_mm,R) 853 854 tmp = tmp + fac*mtmp*gg2 855 856 end if 857 fac = -fac 858 end do 859 call nwpw_timing_end(49) 860 861 nwpw_CWGaussian3 = c*tmp 862 return 863 end 864 865 866* ****************************************************** 867* * * 868* * nwpw_CWGaussian2 * 869* * * 870* ****************************************************** 871* 872* Calculates the two electron two center Gaussian integral 873* 874* // 875* CWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 876* || ------------------------------------ dr dr' 877* // |r-r'| 878* 879* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat) 880* 881* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 882* 883* The normalization constant C_l is defined such at 884* / 885* | g(l,m,s;r) * |r|**l *conjg(Ylm(rhat)) dr = 1 886* / 887* 888 complex*16 function nwpw_CWGaussian2(la,ma,sa,lb,mb,sm,Rab) 889 implicit none 890 integer la,ma,lb,mb 891 real*8 sa,sm,Rab(3) 892 893* *** local variables *** 894 integer l,m,fac,fac2 895 real*8 c,x,y,pi 896 real*8 alpha_aa,alpha_am,alpha_mm 897 real*8 cos_theta,phi,R,gg2 898 complex*16 mtmp,tmp 899 900* **** external functions **** 901 integer nwpw_doublefactorial 902 external nwpw_doublefactorial 903 real*8 nwpw_gaunt,nwpw_GaussBessel 904 external nwpw_gaunt,nwpw_GaussBessel 905 complex*16 Yspherical_lm 906 external Yspherical_lm 907 real*8 nwpw_gaunt_sub,gen_gaunt_coeff_sub 908 external nwpw_gaunt_sub,gen_gaunt_coeff_sub 909 910 911 call nwpw_timing_start(49) 912 pi = 4.0d0*datan(1.0d0) 913 x = dble(nwpw_doublefactorial(2*la+1)) 914 y = dble(nwpw_doublefactorial(2*lb+1)) 915 alpha_aa = dsqrt(0.25d0*(sa*sa + sa*sa)) 916 alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm)) 917 alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm)) 918 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 919 cos_theta = Rab(3)/R 920 921 if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then 922 phi = 0.0d0 923 else 924 phi = datan2(Rab(2),Rab(1)) 925 end if 926 927 if (mod(2*la+lb,2).eq.1) then 928 c = -32.0d0*pi/(x*y) 929 else 930 c = 32.0d0*pi/(x*y) 931 end if 932 933 !phase_factor = (-1)**(m1+l1)/sigma**(l1+l2+1) 934 !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 935 !if (mod(lb+mb,2).eq.1) then 936 937 938 !if (mod(lb+mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 939 if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 940 fac = -1 941 else 942 fac = 1 943 end if 944 945 m = ma + mb 946 tmp = dcmplx(0.0d0,0.0d0) 947 do l = abs(la-lb), (la+lb), 2 948 if (abs(m).le.l) then 949c write(*,*) "l,m=",l,m,la,ma,lb,mb, 950c > nwpw_gaunt(.true.,l,m,la,ma,lb,-mb), 951c > nwpw_gaunt_sub(.true.,l,m,la,ma,lb,-mb), 952c > gen_gaunt_coeff_sub(l,m,la,ma,lb,-mb) 953 mtmp = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb) 954 > *YSpherical_lm(l,m,cos_theta,phi) 955 956 gg2 = nwpw_GaussBessel(la+lb,l,alpha_aa,R) 957 > + nwpw_GaussBessel(la+lb,l,alpha_mm,R) 958 > - 2.0d0*nwpw_GaussBessel(la+lb,l,alpha_am,R) 959 960 tmp = tmp + fac*mtmp*gg2 961 962 end if 963 fac = -fac 964 end do 965 call nwpw_timing_end(49) 966 967 nwpw_CWGaussian2 = c*tmp 968 return 969 end 970 971 972 973* ****************************************************** 974* * * 975* * nwpw_dCWGaussian3 * 976* * * 977* ****************************************************** 978* 979* Calculates the two electron two center Gaussian integral and it's derivative wrt to Rab 980* 981* // 982* dCWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb) = || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb) 983* || ------------------------------------ dr dr' 984* // |r-r'| 985* 986* where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat) 987* 988* and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) ) 989* 990* The normalization constant C_l is defined such at 991* / 992* | g(l,m,s;r) * |r|**l * congj(Ylm(rhat)) dr = 1 993* / 994* 995 subroutine nwpw_dCWGaussian3(la,ma,sa,lb,mb,sb,sm,Rab,W,dW) 996 implicit none 997 integer la,ma,lb,mb 998 real*8 sa,sb,sm,Rab(3) 999 complex*16 W,dW(3) 1000 1001* *** local variables *** 1002 integer l,m,fac 1003 real*8 c,x,y,pi 1004 real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm 1005 real*8 cos_theta,phi,R,gg1,gg2,gg3 1006 complex*16 mtmp,mtmpx,mtmpy,mtmpz,Tx,Ty,Tz 1007 1008* **** external functions **** 1009 integer nwpw_doublefactorial 1010 external nwpw_doublefactorial 1011 real*8 nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel 1012 external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel 1013 complex*16 Yspherical_lm 1014 external Yspherical_lm 1015 1016 call nwpw_timing_start(55) 1017 pi = 4.0d0*datan(1.0d0) 1018 x = dble(nwpw_doublefactorial(2*la+1)) 1019 y = dble(nwpw_doublefactorial(2*lb+1)) 1020 alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb)) 1021 alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm)) 1022 alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb)) 1023 alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm)) 1024 1025 R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3)) 1026 cos_theta = Rab(3)/R 1027 phi = datan2(Rab(2),Rab(1)) 1028 1029 if (mod(2*la+lb,2).eq.1) then 1030 c = -32.0d0*pi/(x*y) 1031 else 1032 c = 32.0d0*pi/(x*y) 1033 end if 1034 1035 !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then 1036 !if (mod(lb+mb,2).eq.1) then 1037 1038 if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then 1039 fac = -1 1040 else 1041 fac = 1 1042 end if 1043 1044 W = dcmplx(0.0d0,0.0d0) 1045 dW(1) = dcmplx(0.0d0,0.0d0) 1046 dW(2) = dcmplx(0.0d0,0.0d0) 1047 dW(3) = dcmplx(0.0d0,0.0d0) 1048 1049 m = ma + mb 1050 do l = abs(la-lb), (la+lb), 2 1051 if (abs(m).le.l) then 1052 gg1 = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb) 1053 call dYspherical_lm(l,m,cos_theta,phi,Tx,Ty,Tz) 1054 mtmp = gg1*Yspherical_lm(l,m,cos_theta,phi) 1055 mtmpx = gg1*Tx 1056 mtmpy = gg1*Ty 1057 mtmpz = gg1*Tz 1058 1059 gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R) 1060 > - nwpw_GaussBessel(la+lb,l,alpha_am,R) 1061 > - nwpw_GaussBessel(la+lb,l,alpha_mb,R) 1062 > + nwpw_GaussBessel(la+lb,l,alpha_mm,R) 1063 1064 gg3 = nwpw_dGaussBessel(la+lb,l,alpha_ab,R) 1065 > - nwpw_dGaussBessel(la+lb,l,alpha_am,R) 1066 > - nwpw_dGaussBessel(la+lb,l,alpha_mb,R) 1067 > + nwpw_dGaussBessel(la+lb,l,alpha_mm,R) 1068 1069 W = W + fac*(mtmp*gg2) 1070 dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3) 1071 dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3) 1072 dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3) 1073 end if 1074 fac = -fac 1075 end do 1076 W = W*c 1077 dW(1) = dW(1)*c 1078 dW(2) = dW(2)*c 1079 dW(3) = dW(3)*c 1080 call nwpw_timing_end(55) 1081 1082 return 1083 end 1084 1085 1086*************************** complex combo versions *********************************** 1087