1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################## 9c ## ## 10c ## subroutine overlap -- p-orbital overlap for pisystem ## 11c ## ## 12c ############################################################## 13c 14c 15c "overlap" computes the overlap for two parallel p-orbitals 16c given the atomic numbers and distance of separation 17c 18c 19 subroutine overlap (atmnum1,atmnum2,rang,ovlap) 20 use units 21 implicit none 22 integer atmnum1 23 integer atmnum2 24 integer na,nb,la,lb 25 real*8 ovlap 26 real*8 rbohr,rang 27 real*8 za,zb,s(3) 28 real*8 zeta(18) 29 save zeta 30c 31c Slater orbital exponents for hydrogen through argon 32c 33 data zeta / 1.000, 1.700, 0.650, 0.975, 1.300, 1.625, 34 & 1.950, 2.275, 2.600, 2.925, 0.733, 0.950, 35 & 1.167, 1.383, 1.600, 1.817, 2.033, 2.250 / 36c 37c 38c principal quantum number from atomic number 39c 40 na = 2 41 nb = 2 42 if (atmnum1 .gt. 10) na = 3 43 if (atmnum2 .gt. 10) nb = 3 44c 45c azimuthal quantum number for p-orbitals 46c 47 la = 1 48 lb = 1 49c 50c orbital exponent from stored ideal values 51c 52 za = zeta(atmnum1) 53 zb = zeta(atmnum2) 54c 55c convert interatomic distance to bohrs 56c 57 rbohr = rang / bohr 58c 59c get pi-overlap via generic overlap integral routine 60c 61 call slater (na,la,za,nb,lb,zb,rbohr,s) 62 ovlap = s(2) 63 return 64 end 65c 66c 67c ############################################################### 68c ## ## 69c ## subroutine slater -- find overlap integrals for STO's ## 70c ## ## 71c ############################################################### 72c 73c 74c "slater" is a general routine for computing the overlap 75c integrals between two Slater-type orbitals 76c 77c literature reference: 78c 79c D. B. Cook, "Structures and Approximations for Electrons in 80c Molecules", Ellis Horwood Limited, Sussex, England, 1978, 81c adapted from the code in Chapter 7 82c 83c variables and parameters: 84c 85c na principle quantum number for first orbital 86c la azimuthal quantum number for first orbital 87c za orbital exponent for the first orbital 88c nb principle quantum number for second orbital 89c lb azimuthal quantum number for second orbital 90c zb orbital exponent for the second orbital 91c r interatomic distance in atomic units 92c s vector containing the sigma-sigma, pi-pi 93c and delta-delta overlaps upon output 94c 95c 96 subroutine slater (na,la,za,nb,lb,zb,r,s) 97 implicit none 98 real*8 rmin,eps 99 parameter (rmin=0.000001d0) 100 parameter (eps=0.00000001d0) 101 integer j,k,m,na,nb,la,lb,ja,jb 102 integer nn,max,maxx,novi 103 integer idsga(5),idsgb(5) 104 integer icosa(2),icosb(2) 105 integer isina(4),isinb(4) 106 integer ia(200),ib(200) 107 real*8 an,ana,anb,anr 108 real*8 rhalf,coef,p,pt 109 real*8 r,za,zb,cjkm 110 real*8 s(3),fact(15) 111 real*8 cbase(20),theta(6) 112 real*8 cosa(2),cosb(2) 113 real*8 sinab(4) 114 real*8 dsiga(5),dsigb(5) 115 real*8 a(20),b(20),c(200) 116 logical done 117 save icosa,icosb 118 save cosa,cosb 119 save idsga,idsgb 120 save dsiga,dsigb 121 save isina,isinb 122 save sinab,theta,fact 123 external cjkm 124 data icosa / 0, 1 / 125 data icosb / 0, 1 / 126 data cosa / 1.0d0, 1.0d0 / 127 data cosb / -1.0d0, 1.0d0 / 128 data idsga / 0, 1, 2, 2, 0 / 129 data idsgb / 0, 1, 2, 0, 2 / 130 data dsiga / 3.0d0, 4.0d0, 3.0d0, -1.0d0, -1.0d0 / 131 data dsigb / 3.0d0,-4.0d0, 3.0d0, -1.0d0, -1.0d0 / 132 data isina / 0, 2, 0, 2 / 133 data isinb / 0, 0, 2, 2 / 134 data sinab / -1.0d0, 1.0d0, 1.0d0, -1.0d0 / 135 data theta / 0.7071068d0, 1.2247450d0, 0.8660254d0, 136 & 0.7905694d0, 1.9364916d0, 0.9682458d0 / 137 data fact / 1.0d0, 1.0d0, 2.0d0, 6.0d0, 24.0d0, 120.0d0, 138 & 720.0d0, 5040.0d0, 40320.0d0, 362880.0d0, 139 & 3628800.0d0, 39916800.0d0, 479001600.0d0, 140 & 6227020800.0d0, 87178291200.0d0 / 141c 142c 143c zero out the overlap integrals 144c 145 done = .false. 146 s(1) = 0.0d0 147 s(2) = 0.0d0 148 s(3) = 0.0d0 149 ana = (2.0d0*za)**(2*na+1) / fact(2*na+1) 150 anb = (2.0d0*zb)**(2*nb+1) / fact(2*nb+1) 151c 152c orbitals are on the same atomic center 153c 154 if (r .lt. rmin) then 155 anr = 1.0d0 156 j = na + nb + 1 157 s(1) = fact(j) / ((za+zb)**j) 158 an = sqrt(ana*anb) 159 do novi = 1, 3 160 s(novi) = s(novi) * an * anr 161 end do 162 return 163 end if 164c 165c compute overlap integrals for general case 166c 167 rhalf = 0.5d0 * r 168 p = rhalf * (za+zb) 169 pt = rhalf * (za-zb) 170 nn = na + nb 171 call aset (p,nn,a) 172 call bset (pt,nn,b) 173 k = na - la 174 m = nb - lb 175 max = k + m + 1 176 do j = 1, max 177 ia(j) = j - 1 178 ib(j) = max - j 179 cbase(j) = cjkm(j-1,k,m) 180 c(j) = cbase(j) 181 end do 182 maxx = max 183 if (la .eq. 1) then 184 call polyp (c,ia,ib,maxx,cosa,icosa,icosb,2) 185 else if (la .eq. 2) then 186 call polyp (c,ia,ib,maxx,dsiga,idsga,idsgb,5) 187 end if 188 if (lb .eq. 1) then 189 call polyp (c,ia,ib,maxx,cosb,icosa,icosb,2) 190 else if (lb .eq. 2) then 191 call polyp (c,ia,ib,maxx,dsigb,idsga,idsgb,5) 192 end if 193 novi = 1 194 do while (.not. done) 195 do j = 1, maxx 196 ja = ia(j) + 1 197 jb = ib(j) + 1 198 coef = c(j) 199 if (abs(coef) .ge. eps) then 200 s(novi) = s(novi) + coef*a(ja)*b(jb) 201 end if 202 end do 203 ja = la*(la+1)/2 + novi 204 jb = lb*(lb+1)/2 + novi 205 s(novi) = s(novi) * theta(ja) * theta(jb) 206 if (novi.eq.1 .and. la.ne.0 .and. lb.ne.0) then 207 maxx = max 208 do j = 1, maxx 209 c(j) = cbase(j) 210 end do 211 call polyp (c,ia,ib,maxx,sinab,isina,isinb,4) 212 if (la .eq. 2) then 213 call polyp (c,ia,ib,maxx,cosa,icosa,icosb,2) 214 end if 215 if (lb .eq. 2) then 216 call polyp (c,ia,ib,maxx,cosb,icosa,icosb,2) 217 end if 218 novi = 2 219 else if (novi.eq.2 .and. la.eq.2 .and. lb.eq.2) then 220 maxx = max 221 do j = 1, maxx 222 c(j) = cbase(j) 223 end do 224 call polyp (c,ia,ib,maxx,sinab,isina,isinb,4) 225 call polyp (c,ia,ib,maxx,sinab,isina,isinb,4) 226 novi = 3 227 else 228 anr = rhalf**(na+nb+1) 229 an = sqrt(ana*anb) 230 do novi = 1, 3 231 s(novi) = s(novi) * an * anr 232 end do 233 done = .true. 234 end if 235 end do 236 return 237 end 238c 239c 240c ################################################################ 241c ## ## 242c ## subroutine polyp -- polynomial product for STO overlap ## 243c ## ## 244c ################################################################ 245c 246c 247c "polyp" is a polynomial product routine that multiplies two 248c algebraic forms 249c 250c 251 subroutine polyp (c,ia,ib,max,d,iaa,ibb,n) 252 implicit none 253 integer i,j,k,m,max,n 254 integer ia(200),ib(200) 255 integer iaa(*),ibb(*) 256 real*8 c(200),d(*) 257c 258c 259 do j = 1, max 260 do k = 1, n 261 i = n - k + 1 262 m = (i-1)*max + j 263 c(m) = c(j) * d(i) 264 ia(m) = ia(j) + iaa(i) 265 ib(m) = ib(j) + ibb(i) 266 end do 267 end do 268 max = n * max 269 return 270 end 271c 272c 273c ############################################################## 274c ## ## 275c ## function cjkm -- coefficients of spherical harmonics ## 276c ## ## 277c ############################################################## 278c 279c 280c "cjkm" computes the coefficients of spherical harmonics 281c expressed in prolate spheroidal coordinates 282c 283c 284 function cjkm (j,k,m) 285 implicit none 286 integer i,j,k,m 287 integer min,max 288 integer id,idd,ip1 289 real*8 cjkm,b1,b2,sum 290 real*8 fact(15) 291 save fact 292 data fact / 1.0d0, 1.0d0, 2.0d0, 6.0d0, 24.0d0, 120.0d0, 293 & 720.0d0, 5040.0d0, 40320.0d0, 362880.0d0, 294 & 3628800.0d0, 39916800.0d0, 479001600.0d0, 295 & 6227020800.0d0, 87178291200.0d0 / 296c 297c 298 min = 1 299 if (j .gt. m) min = j - m + 1 300 max = j + 1 301 if (k .lt. j) max = k + 1 302 sum = 0.0d0 303 do ip1 = min, max 304 i = ip1 - 1 305 id = k - i + 1 306 b1 = fact(k+1) / (fact(i+1)*fact(id)) 307 if (j .lt. i) then 308 b2 = 1.0d0 309 else 310 id = m - (j-i) + 1 311 idd = j - i + 1 312 b2 = fact(m+1) / (fact(idd)*fact(id)) 313 end if 314 sum = sum + b1*b2*(-1.0d0)**i 315 end do 316 cjkm = sum * (-1.0d0)**(m-j) 317 return 318 end 319c 320c 321c ########################################################### 322c ## ## 323c ## subroutine aset -- get "A" functions by recursion ## 324c ## ## 325c ########################################################### 326c 327c 328c "aset" computes by recursion the A functions used in the 329c evaluation of Slater-type (STO) overlap integrals 330c 331c 332 subroutine aset (alpha,n,a) 333 implicit none 334 integer i,n 335 real*8 alpha,alp 336 real*8 a(20) 337c 338c 339 alp = 1.0d0 / alpha 340 a(1) = exp(-alpha) * alp 341 do i = 1, n 342 a(i+1) = a(1) + dble(i)*a(i)*alp 343 end do 344 return 345 end 346c 347c 348c ########################################################### 349c ## ## 350c ## subroutine bset -- get "B" functions by recursion ## 351c ## ## 352c ########################################################### 353c 354c 355c "bset" computes by downward recursion the B functions used 356c in the evaluation of Slater-type (STO) overlap integrals 357c 358c 359 subroutine bset (beta,n,b) 360 implicit none 361 real*8 eps 362 parameter (eps=0.000001d0) 363 integer i,j,n 364 real*8 beta,bmax 365 real*8 betam,d1,d2 366 real*8 b(20) 367 external bmax 368c 369c 370 if (abs(beta) .lt. eps) then 371 do i = 1, n+1 372 b(i) = 2.0d0 / dble(i) 373 if ((i/2)*2 .eq. i) b(i) = 0.0d0 374 end do 375 else if (abs(beta) .gt. (dble(n)/2.3d0)) then 376 d1 = exp(beta) 377 d2 = 1.0d0 / d1 378 betam = 1.0d0 / beta 379 b(1) = (d1-d2) * betam 380 do i = 1, n 381 d1 = -d1 382 b(i+1) = (d1-d2+dble(i)*b(i)) * betam 383 end do 384 else 385 b(n+1) = bmax(beta,n) 386 d1 = exp(beta) 387 d2 = 1.0d0 / d1 388 if ((n/2)*2 .ne. n) d1 = -d1 389 do i = 1, n 390 j = n - i + 1 391 d1 = -d1 392 b(j) = (d1+d2+beta*b(j+1)) / dble(j) 393 end do 394 end if 395 return 396 end 397c 398c 399c ############################################################## 400c ## ## 401c ## function bmax -- find maximum order of "B" functions ## 402c ## ## 403c ############################################################## 404c 405c 406c "bmax" computes the maximum order of the B functions needed 407c for evaluation of Slater-type (STO) overlap integrals 408c 409c 410 function bmax (beta,n) 411 implicit none 412 real*8 eps 413 parameter (eps=0.0000001d0) 414 integer n 415 real*8 bmax,beta 416 real*8 b,top,bot 417 real*8 sum,fi 418 real*8 sign,term 419 logical done 420c 421c 422 done = .false. 423 b = beta**2 424 top = dble(n) + 1.0d0 425 sum = 1.0d0 / top 426 fi = 2.0d0 427 sign = 2.0d0 428 if ((n/2)*2 .ne. n) then 429 top = top + 1.0d0 430 sum = beta / top 431 fi = fi + 1.0d0 432 sign = -2.0d0 433 end if 434 term = sum 435 do while (.not. done) 436 bot = top + 2.0d0 437 term = term * b * top / (fi*(fi-1.0d0)*bot) 438 sum = sum + term 439 if (abs(term) .le. eps) then 440 done = .true. 441 else 442 fi = fi + 2.0d0 443 top = bot 444 end if 445 end do 446 bmax = sign * sum 447 return 448 end 449