1 subroutine espod(x,y,z,epot,idebug, 2 & p) 3c THIS IS REALLY espot 4 implicit double precision(a-h,o-z) 5 parameter (numatm=2000) 6 parameter (numprm=1600) 7 parameter (numcex=numprm*3) 8 parameter (numtmp=4000) 9 logical oeerst 10 integer shella,shelln,shellt,shellc,shladf,aos 11 common/b/exx(numcex), 12 & c1(numcex),c2(numcex),c3(numprm),c4(numprm),c5(numprm), 13 & shladf(numprm),gx(numprm),gy(numprm),gz(numprm), 14 & jan(numprm),shella(numprm),shelln(numprm),shellt(numprm), 15 & shellc(numprm),aos(numprm),nshell,maxtyp 16 common /bounds/ ilow(5,5),iupp(5) 17 common /indxe/ jx(35),jy(35),jz(35),ix(20),iy(20),iz(20) 18 common /strt/ oeerst 19 common /moldat/ natoms, norbs, nelecs,nat(numatm) 20 common /pseudo / ipseud,ivale(numatm) 21 common /coord / xyz(3,numatm) 22 common /orbhlp/ mxorb,iuhf,ispd 23 common /extchg/ exchg(3,3),iexchg(3),nexchg 24 dimension tt(4),zz(4),pre(45),ca(35),cb(35), 25 & coeffx(192),coeffy(192),coeffz(192) 26 dimension h(numtmp),f(numtmp),h2x(9),h2y(9),h2z(9), 27 & h3x(16),h3y(16),h3z(16) 28 dimension p(*) 29 data jx/1,2,1,1,3,1,1,2,2,1,4,1,1,2,3,3,2,1,1,2, 30 & 5,1,1,4,4,2,1,2,1,3,3,1,3,2,2/ 31 data jy/1,1,2,1,1,3,1,2,1,2,1,4,1,3,2,1,1,2,3,2, 32 & 1,5,1,2,1,4,4,1,2,3,1,3,2,3,2/ 33 data jz/1,1,1,2,1,1,3,1,2,2,1,1,4,1,1,2,3,3,2,2, 34 & 1,1,5,1,2,1,2,4,4,1,3,3,2,2,3/ 35 data ix/0,4,0,0,8,0,0,4,4,0,12,0,0,4,8,8,4,0,0,4/ 36 data iy/0,0,4,0,0,8,0,4,0,4,0,12,0,8,4,0,0,4,8,4/ 37 data iz/0,0,0,4,0,0,8,0,4,4,0,0,12,0,0,4,8,8,4,4/ 38c iupp 1 4 10 20 35 39c s 3p 6d 10f 15g 40c 41c ilow 0 1 2 3 4 <- shellt 42c--------------------------------------- 43c 0 1(s) 1 (sp) 1 (spd) 1 1 44c 1 1 2 (p) 5 11 21 45c 2 1 1 5 (d) 11 21 46c 3 1 1 5 11 (f) 21 47c 4 1 1 5 11 21 (g) 48c--------------------------------------- 49c ^ 50c | shellc 51c 52 data iupp/1,4,10,20,35/ 53 data ilow/5*1,1,2,5,11,21,1,1,5,11,21,1,1,5,11,21, 54 & 1,1,5,11,21/ 55 56 if (oeerst) then 57 call epint 58 call setcon 59 oeerst = .false. 60 endif 61 62 epot = 0.0d0 63 pi = 4.d0*datan(1.d0) 64c 65c loop over shell pairs. 66c 67 do ishell = 1, nshell 68 xa = gx(ishell) 69 ya = gy(ishell) 70 za = gz(ishell) 71 ifrst = shella(ishell) 72 ilast = ifrst + shelln(ishell) - 1 73 itype = shellt(ishell) 74 itype1 = itype + 1 75 istart = ilow(itype1,shellc(ishell)+1) 76 iend = iupp(itype1) 77 iendt = iend 78 79 do jshell = 1,ishell 80 xb = gx(jshell) 81 yb = gy(jshell) 82 zb = gz(jshell) 83 jfsrt = shella(jshell) 84 jlast = jfsrt + shelln(jshell) - 1 85 jtype = shellt(jshell) 86 jtype1 = jtype + 1 87 jstart = ilow(jtype1,shellc(jshell)+1) 88 jend = iupp(jtype1) 89 ijtyp = itype1 + jtype1 - 1 90 ndim = (iend-istart+1)*(jend-jstart+1) 91 iminj = iabs(ishell - jshell) 92 n = (itype + jtype) / 2 + 1 93 xt = xb - xa 94 yt = yb - ya 95 zt = zb - za 96 rsq = xt*xt + yt*yt + zt*zt 97 if (ndim.gt.numtmp) 98 & print*,'WARNING: exceding array limit in espot' 99 do ii=1,ndim 100 h(ii) = 0.0d0 101 end do 102 jendt = jend 103c 104c loop over primitive gaussians. 105c 106 do igauss = ifrst, ilast 107 ei = exx(igauss) 108 call fcij(itype,ifrst,igauss,shladf(ishell),ca) 109 do jgauss = jfsrt,jlast 110 iend = iendt 111 jend = jendt 112 ej = exx(jgauss) 113 call fcij(jtype,jfsrt,jgauss,shladf(jshell),cb) 114 ep = ei + ej 115 exptmp = ej*ei*rsq / ep 116 if (exptmp.ge.600.0d0) goto 100 117 expp = 2.0d0*dexp(-exptmp) 118 tx = (ei*xa + ej*xb) / ep 119 ty = (ei*ya + ej*yb) / ep 120 tz = (ei*za + ej*zb) / ep 121 xta = tx - xa 122 yta = ty - ya 123 zta = tz - za 124 xtb = tx - xb 125 ytb = ty - yb 126 ztb = tz - zb 127 call coeffs(coeffx,xta,xtb,itype1,jtype1) 128 call coeffs(coeffy,yta,ytb,itype1,jtype1) 129 call coeffs(coeffz,zta,ztb,itype1,jtype1) 130 if (ndim.gt.numtmp) 131 & print*,'WARNING: exceding array limit in espot' 132 do ii=1,ndim 133 f(ii) = 0.0d0 134 end do 135 xct = x - tx 136 yct = y - ty 137 zct = z - tz 138 expont = ep * (xct * xct + yct * yct + zct * zct) 139 call rys(n,expont,tt,zz) 140 epp = 0.5d0 / ep 141 call calct(pre,epp,ijtyp) 142 do ii = 1, n 143 fac1 = (ep + ep)*tt(ii) 144 zf = pi * expp * zz(ii) / ep 145 call twocen(h2x,xct,1.d0,pre,fac1,ijtyp) 146 call twocen(h2y,yct,1.d0,pre,fac1,ijtyp) 147 call twocen(h2z,zct,zf,pre,fac1,ijtyp) 148 call thrcen(h3x,h2x,coeffx,itype1,jtype1) 149 call thrcen(h3y,h2y,coeffy,itype1,jtype1) 150 call thrcen(h3z,h2z,coeffz,itype1,jtype1) 151 152 k = 0 153 do i = istart, iend 154 do j = jstart, jend 155 k = k + 1 156 if (k.gt.numtmp) 157 & print*,'WARNING: exceding array limit in espot' 158 f(k) = f(k)- 159 & (h3x(jx(j)+ix(i))*h3y(jy(j)+iy(i))*h3z(jz(j)+iz(i))) 160 end do 161 end do 162 end do 163 164 k = 0 165 do i = istart, iend 166 do j = jstart,jend 167 k = k + 1 168 if (k.gt.numtmp) 169 & print*,'WARNING: exceding array limit in espot' 170 h(k) = h(k) + f(k)*ca(i)*cb(j) 171 end do 172 end do 173 kend = k 174 175100 continue 176 end do 177 end do 178 179 call purdf(itype,jtype,istart,jstart,iend,jend,h) 180c 181c Density matrix contribution 182c 183 184 ii = 0 185 ist = aos(ishell) - 1 186 jst = aos(jshell) - 1 187 do i = istart,iend 188 do j = jstart,jend 189 ii = ii + 1 190 if (iminj.eq.0) then 191 a1or2 = 1.0d0 192 if (i-j.lt.0) then 193 pt = p((ist+i-1)*mxorb + (jst+j)) 194 else 195 pt = p((jst+j-1)*mxorb + (ist+i)) 196 endif 197 else 198 a1or2 = 2.0d0 199 pt = p((jst+j-1)*mxorb + (ist+i)) 200 endif 201 phelp = pt*a1or2 202 epot = epot + phelp*h(ii) 203 if (idebug.eq.1) print*,'v(',ist+i,',',jst+j,')=',h(ii) 204 end do 205 end do 206 207 jend = jendt 208 iend = iendt 209 end do 210 end do 211 212c sum in nuclear contribution 213 214c 215 do ii = 1 ,natoms 216 xt = xyz(1,ii) - x 217 yt = xyz(2,ii) - y 218 zt = xyz(3,ii) - z 219 rsq = xt*xt + yt*yt + zt*zt 220 if (rsq.ge.1.0d-8) then 221 if (ipseud.eq.1) then 222 epot = epot + dfloat(ivale(ii))*dsqrt(1.0d0 / rsq) 223 else 224 epot = epot + dfloat(nat(ii))*dsqrt(1.0d0 / rsq) 225 endif 226 endif 227 end do 228 229 if (nexchg.ne.0) then 230 do ii=1,nexchg 231 xt = exchg(1,ii) - x 232 yt = exchg(2,ii) - y 233 zt = exchg(3,ii) - z 234 rsq = xt*xt + yt*yt + zt*zt 235 if (rsq.ge.1.0d-8) then 236 if (iexchg(ii).eq.1) then 237 epot = epot + dsqrt(1.0d0 / rsq) 238 else 239 epot = epot - dsqrt(1.0d0 / rsq) 240 endif 241 endif 242 end do 243 endif 244 245 return 246 end 247 248 subroutine espgrdd(npts1,npts2,npts3,idebug, 249 & denn,p,xden,yden,zden) 250c THIS IS REALLY espgrd 251 implicit double precision(a-h,o-z) 252 parameter (numatm=2000) 253 parameter (numprm=1600) 254 parameter (numcex=numprm*3) 255 parameter (nmcex2=numcex*numcex) 256 parameter (numtmp=4000) 257 character monstr*100 258 character*2 gstr 259 logical oeerst 260 integer shella,shelln,shellt,shellc,shladf,aos 261 common/b/exx(numcex), 262 & c1(numcex),c2(numcex),c3(numprm),c4(numprm),c5(numprm), 263 & shladf(numprm),gx(numprm),gy(numprm),gz(numprm), 264 & jan(numprm),shella(numprm),shelln(numprm),shellt(numprm), 265 & shellc(numprm),aos(numprm),nshell,maxtyp 266 common /bounds/ ilow(5,5),iupp(5) 267 common /indxe/ jx(35),jy(35),jz(35),ix(20),iy(20),iz(20) 268 common /strt/ oeerst 269 common /moldat/ natoms, norbs, nelecs,nat(numatm) 270 common /pseudo / ipseud,ivale(numatm) 271 common /coord / xyz(3,numatm) 272 common /orbhlp/ mxorb,iuhf,ispd 273 common /extchg/ exchg(3,3),iexchg(3),nexchg 274 common /plane/ px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat 275 common /grdhlp/ mx3d,mx3d2 276 dimension tt(4),zz(4),pre(45),ca(35),cb(35), 277 & coeffx(192),coeffy(192),coeffz(192) 278 dimension h(numtmp),f(numtmp),h2x(9),h2y(9),h2z(9), 279 & h3x(16),h3y(16),h3z(16),v3(3) 280 dimension p(*),denn(*),xden(*),yden(*),zden(*) 281 282 if (oeerst) then 283 call epint 284 call setcon 285 oeerst = .false. 286 do i=1,192 287 coeffx(i) = 0.d0 288 end do 289 endif 290 291 292 call pregrd(npts1,npts2,npts3,xden,yden,zden) 293 294 do kc=1,npts3 295 do ic=1,npts1 296 do jc=1,npts2 297 iadrs = (kc-1)*mx3d2 + ij 298 denn(iadrs) = 0.0d0 299 end do 300 end do 301 end do 302 303 pi = 4.d0*datan(1.d0) 304c 305c loop over shell pairs. 306c 307 nn = 0 308 iup = 0 309 do ishell = 1, nshell 310 311 xa = gx(ishell) 312 ya = gy(ishell) 313 za = gz(ishell) 314 315 ifrst = shella(ishell) 316 ilast = ifrst + shelln(ishell) - 1 317 itype = shellt(ishell) 318 itype1 = itype + 1 319 istart = ilow(itype1,shellc(ishell)+1) 320 iend = iupp(itype1) 321 iendt = iend 322 323 do jshell = 1,ishell 324 325 xb = gx(jshell) 326 yb = gy(jshell) 327 zb = gz(jshell) 328 329 jfsrt = shella(jshell) 330 jlast = jfsrt + shelln(jshell) - 1 331 jtype = shellt(jshell) 332 jtype1 = jtype + 1 333 jstart = ilow(jtype1,shellc(jshell)+1) 334 jend = iupp(jtype1) 335 ijtyp = itype1 + jtype1 - 1 336 ndim = (iend-istart+1)*(jend-jstart+1) 337 iminj = iabs(ishell - jshell) 338 339 n = (itype + jtype) / 2 + 1 340 341 xt = xb - xa 342 yt = yb - ya 343 zt = zb - za 344 345 rsq = xt*xt + yt*yt + zt*zt 346 347 if (ndim.gt.numtmp) 348 & print*,'WARNING: exceding array limit in espot' 349 350c deze jendt moet niet in binnenste loop ?? 351 jendt = jend 352 353 do kc=1,npts3 354 355 monstr = 'Progress Monitor [1-' // gstr(npts3)//'] ' // 356 & gstr(kc) 357 ij = 0 358 359 do ic=1,npts1 360 361 do jc=1,npts2 362 363 364 ij = ij + 1 365 iadrs = (kc-1)*mx3d2 + ij 366 367 x = xden(iadrs) 368 y = yden(iadrs) 369 z = zden(iadrs) 370 371 372 do ii=1,ndim 373 h(ii) = 0.0d0 374 end do 375 376 377c loop over primitive gaussians. 378 379 do igauss = ifrst, ilast 380 381 ei = exx(igauss) 382 call fcij(itype,ifrst,igauss,shladf(ishell),ca) 383 384 do jgauss = jfsrt,jlast 385 386 iend = iendt 387 jend = jendt 388 389 ej = exx(jgauss) 390 call fcij(jtype,jfsrt,jgauss,shladf(jshell),cb) 391 392 ep = ei + ej 393 epr = 1.0d0 / ep 394 395 exptmp = ej*ei*rsq*epr 396 397 if (exptmp.ge.600.0d0) goto 100 398 399 expp = 2.0d0*dexp(-exptmp) 400 tx = (ei*xa + ej*xb) * epr 401 ty = (ei*ya + ej*yb) * epr 402 tz = (ei*za + ej*zb) * epr 403 404 xta = tx - xa 405 yta = ty - ya 406 zta = tz - za 407 408 xtb = tx - xb 409 ytb = ty - yb 410 ztb = tz - zb 411 412 call coeffs(coeffx,xta,xtb,itype1,jtype1) 413 call coeffs(coeffy,yta,ytb,itype1,jtype1) 414 call coeffs(coeffz,zta,ztb,itype1,jtype1) 415 416 if (ndim.gt.numtmp) 417 & print*,'WARNING: exceding array limit in espot' 418 419 do ii=1,ndim 420 f(ii) = 0.0d0 421 end do 422 423 xct = x - tx 424 yct = y - ty 425 zct = z - tz 426 427 expont = ep * (xct * xct + yct * yct + zct * zct) 428 429 call rys(n,expont,tt,zz) 430 epp = 0.5d0 * epr 431 call calct(pre,epp,ijtyp) 432 ep2 = ep + ep 433 pee = pi*expp*epr 434 435 do ii = 1, n 436 fac1 = ep2 * tt(ii) 437 zf = pee * zz(ii) 438 call twocen(h2x,xct,1.d0,pre,fac1,ijtyp) 439 call twocen(h2y,yct,1.d0,pre,fac1,ijtyp) 440 call twocen(h2z,zct,zf,pre,fac1,ijtyp) 441 call thrcen(h3x,h2x,coeffx,itype1,jtype1) 442 call thrcen(h3y,h2y,coeffy,itype1,jtype1) 443 call thrcen(h3z,h2z,coeffz,itype1,jtype1) 444 445 k = 0 446 do i = istart, iend 447 do j = jstart, jend 448 k = k + 1 449 if (k.gt.numtmp) print*, 450 & 'WARNING: exceding array limit in espot' 451 f(k) = f(k) - (h3x(jx(j)+ix(i))* 452 & h3y(jy(j)+iy(i))* 453 & h3z(jz(j)+iz(i))) 454 455 end do 456 end do 457 458 end do 459 460 k = 0 461 do i = istart, iend 462 do j = jstart,jend 463 k = k + 1 464 465 if (k.gt.numtmp) print*, 466 & 'WARNING: exceding array limit in espot' 467 468 h(k) = h(k) + f(k)*ca(i)*cb(j) 469 470 end do 471 end do 472 473 474100 continue 475 476 end do 477 end do 478 479 call purdf(itype,jtype,istart,jstart,iend,jend,h) 480c 481c Density matrix contribution 482c 483 484 ii = 0 485 ist = aos(ishell) - 1 486 jst = aos(jshell) - 1 487 488 do i = istart,iend 489 do j = jstart,jend 490 ii = ii + 1 491 if (iminj.eq.0) then 492 a1or2 = 1.0d0 493 if (i-j.lt.0) then 494 pt = p((ist+i-1)*mxorb + (jst+j)) 495 else 496 pt = p((jst+j-1)*mxorb + (ist+i)) 497 endif 498 else 499 a1or2 = 2.0d0 500 pt = p((jst+j-1)*mxorb + (ist+i)) 501 endif 502 phelp = pt*a1or2 503 denn(iadrs) = denn(iadrs) + phelp*h(ii) 504 505 if (idebug.eq.1) 506 & print*,'v(',ist+i,',',jst+j,')=',h(ii) 507 508 end do 509 end do 510 511 jend = jendt 512 iend = iendt 513 514c end do's kc,ic,jc 515 516 end do 517 end do 518 end do 519 520c end do's ishell, jshell 521 522 end do 523 end do 524 525c sum in nuclear contribution 526 527 do kc=1,npts3 528 ij = 0 529 do ic=1,npts1 530 do jc=1,npts2 531 532 ij = ij + 1 533 iadrs = (kc-1)*mx3d2 + ij 534 535 x = xden(iadrs) 536 y = yden(iadrs) 537 z = zden(iadrs) 538 539 do ii = 1,natoms 540 541 xt = xyz(1,ii) - x 542 yt = xyz(2,ii) - y 543 zt = xyz(3,ii) - z 544 545 rsq = xt*xt + yt*yt + zt*zt 546 547 if (rsq.ge.1.0d-8) then 548 if (ipseud.eq.1) then 549 denn(iadrs) = denn(iadrs) + 550 & dfloat(ivale(ii))*dsqrt(1.0d0 / rsq) 551 else 552 denn(iadrs) = denn(iadrs) + 553 & dfloat(nat(ii))*dsqrt(1.0d0 / rsq) 554 endif 555 endif 556 end do 557 558 if (nexchg.ne.0) then 559 do ii=1,nexchg 560 xt = exchg(1,ii) - x 561 yt = exchg(2,ii) - y 562 zt = exchg(3,ii) - z 563 rsq = xt*xt + yt*yt + zt*zt 564 if (rsq.ge.1.0d-8) then 565 if (iexchg(ii).eq.1) then 566 denn(iadrs) = denn(iadrs) + 567 & dsqrt(1.0d0 / rsq) 568 else 569 denn(iadrs) = denn(iadrs) - 570 & dsqrt(1.0d0 / rsq) 571 endif 572 endif 573 end do 574 endif 575 576 end do 577 end do 578 end do 579 580 return 581 end 582 583 subroutine purdf(itype,jtype,istart,jstart,iend,jend,h) 584c 585c only iend and jend are really nescessary to return 586c 587c iupp 1 4 10 20 35 588c s 3p 6d 10f 15g 589c 590c ilow 0 1 2 3 4 <- shellt 591c--------------------------------------- 592c 0 1(s) 1 (sp) 1 (spd) 1 1 593c 1 1 2 (p) 5 11 21 594c 2 1 1 5 (d) 11 21 595c 3 1 1 5 11 (f) 21 596c 4 1 1 5 11 21 (g) 597c--------------------------------------- 598c ^ 599c | shellc 600c 601c question does h come filled in shell form ?, eg d is really spd 602c with start at 5 (instead of 1) (we treat d such but not f?) 603c 604c xxxx,yyyy,zzzz *1 605c xxxy ... *d7 606c xxyy ... *d5*d7/d3 607c xxyz ... *d5*d7 608c 609c via fcij ca and cb have these, and therfor h() has them 610c 611 implicit double precision (a-h,o-z) 612 common /slagau/ ihasd,isgau,ido5d,ido7f,ido9g,ihasg 613 common /intcon/ d3,d5,d7,r1,r2,r3,r4,r3ov2,z1,z2,z3, 614 & g1,g2,g3,g4,g5,g6 615 logical debug 616 dimension h(*),inc(14) 617 debug = .false. 618c 619c convert gaussians to pure angular functions. 620c (6D -> 5D, 10F -> 7F, 15G -> 9G) 621c 622 623 b1 = 1.0d0/d7 624 b2 = d3/(d5*d7) 625 b3 = 1.0d0/(d5*d7) 626 627 idim = iend - istart + 1 628 jdim = jend - jstart + 1 629 630 if (ido5d.eq.1.and.jtype.eq.2) then 631c 632c The row side of the matrix: pure d 633c 634 635 i1 = 5 - jstart + 1 636 637c when d shell jstart = 5 , i1 = 1, jend = 10, jdim = 6 638c when spd shell jstart = 1 , i1 = 5, jend = 10, jdim = 10 639c single d shell is stored in H() as: 640c 641c d1 d2 d3 d4 d5 d6 NOT s px py pz d1 d2 etc 642 643 do i=1,idim 644 dz2 = h(i1+2) - 0.5d0*(h(i1) + h(i1+1)) 645 dx2y2 = r3ov2*(h(i1) - h(i1+1)) 646 h(i1 ) = dz2 647 h(i1+1) = h(i1+4) 648 h(i1+2) = h(i1+5) 649 h(i1+4) = h(i1+3) 650 h(i1+3) = dx2y2 651 i1 = i1 + jdim 652 end do 653 654 endif 655 656 if (ido7f.eq.1.and.jtype.eq.3) then 657c 658c The row side af the matrix: pure f 659c 660c when f shell jstart = 11 , jend = 20, jdim = 10 661c 662 i1 = 0 663 do i=1,idim 664 665 f0 = h(i1+3) - r2*(h(i1+6) + h(i1+9)) 666 f1p = r4*(z1*h(i1+7) - h(i1+1) - z2*h(i1+4)) 667 f1m = r4*(z1*h(i1+8) - h(i1+2) - z2*h(i1+5)) 668 f2p = r3ov2*(h(i1+6) - h(i1+9)) 669 f2m = h(i1+10) 670 f3p = r1*(h(i1+1) - z3*h(i1+4)) 671 f3m = r1*(z3*h(i1+5) - h(i1+2)) 672 673 h(i1+1) = f0 674 h(i1+2) = f1p 675 h(i1+3) = f1m 676 h(i1+4) = f2p 677 h(i1+5) = f2m 678 h(i1+6) = f3p 679 h(i1+7) = f3m 680 681 i1 = i1 + jdim 682 683 end do 684 endif 685 686 if (ido9g.eq.1.and.jtype.eq.4) then 687c 688c The row side af the matrix: pure g 689c 690c when g shell jstart = 21 , jend = 35, jdim = 15 691c 692 i1 = 0 693 do i=1,idim 694 695 g0 = 696 & 0.375d0*h(i1+1) + 0.375d0* h(i1+2) + h(i1+3) 697 & -3.0d0*b2*h(i1+11) - 3.0d0*b2*h(i1+12) + 0.75d0*b2*h(i1+10) 698 699 g1p = 700 & g1*(4.0d0*b1*h(i1+8) - 3.0d0*b3*h(i1+14) - 3.0d0*b1* h(i1+5)) 701 702 g1m = 703 & g1*(4.0d0*b1*h(i1+9) - 3.0d0*b3*h(i1+13) - 3.0d0*b1* h(i1+7)) 704 705 g2p = 706 & g2*(6.0d0*b2*h(i1+11) - 6.0d0*b2*h(i1+12) - h(i1+1) 707 & + h(i1+2)) 708 709 g2m = 710 & g3*(6.0d0*b3*h(i1+15) - b1* h(i1+4) - b1* h(i1+6)) 711 712 g3p = 713 & g4*( b1* h(i1+5) - 3.0d0*b3*h(i1+14)) 714 715 g3m = 716 & g4*(3.0d0*b3*h(i1+13) - b1* h(i1+7)) 717 718 g4p = 719 & g5*( h(i1+1) - 6.0d0*b2*h(i1+10) + h(i1+2)) 720 721 g4m = 722 & g6*b1*( h(i1+4) - h(i1+6)) 723 724 h(i1+1) = g0 725 h(i1+2) = g1p 726 h(i1+3) = g1m 727 h(i1+4) = g2p 728 h(i1+5) = g2m 729 h(i1+6) = g3p 730 h(i1+7) = g3m 731 h(i1+8) = g4p 732 h(i1+9) = g4m 733 734 i1 = i1 + jdim 735 736 end do 737 endif 738c 739c the row size of the matrix has changed, so get rid of 740c the superflous functions. 741c 742 if ((ido5d.eq.1.and.jtype.eq.2).or. 743 & (ido7f.eq.1.and.jtype.eq.3).or. 744 & (ido9g.eq.1.and.jtype.eq.4)) then 745 746 if (jtype.eq.2) then 747 jendp = 9 748 else if (jtype.eq.3) then 749c when f shell jstart = 11 , jend = 17, jpure = jdim = 7, jend = 17 750 jendp = 17 751 else if (jtype.eq.4) then 752c when g shell jstart = 21 , jendp = 29, jpure = jdim = 9, jend = 29 753 jendp = 29 754 endif 755 if (debug) print*,'row.old.new ',itype,jtype,jend,jendp 756 jrpure = jendp - jstart + 1 757 758 i1 = 0 759 i2 = 0 760 761 do i=1,idim 762 do j=1,jrpure 763 h(i2+j) = h(i1+j) 764 end do 765 i1 = i1 + jdim 766 i2 = i2 + jrpure 767 end do 768 jend = jendp 769 jdim = jrpure 770 771 endif 772c 773c transformation at row side complete, start of column side 774c 775 776 if (ido5d.eq.1.and.itype.eq.2) then 777c 778c The column side of the matrix: pure d 779c 780 i1 = (5-istart)*jdim + 1 781 do i=1,5 782 inc(i) = i*jdim 783 end do 784 785 do j=1,jdim 786 dz2 = h(i1+inc(2)) - 0.5d0*(h(i1) + h(i1+inc(1))) 787 dx2y2 = r3ov2*(h(i1) - h(i1+inc(1))) 788 h(i1 ) = dz2 789 h(i1+inc(1)) = h(i1+inc(4)) 790 h(i1+inc(2)) = h(i1+inc(5)) 791 h(i1+inc(4)) = h(i1+inc(3)) 792 h(i1+inc(3)) = dx2y2 793 i1 = i1 + 1 794 end do 795 796 endif 797 798 if (ido7f.eq.1.and.itype.eq.3) then 799c 800c The column side af the matrix: pure F 801c 802 i1 = 1 803 do i=1,9 804 inc(i) = i*jdim 805 end do 806 807 do j=1,jdim 808 809 f0 = h(i1+inc(2)) - r2*(h(i1+inc(5)) + h(i1+inc(8))) 810 f1p = r4*(z1*h(i1+inc(6)) - h(i1) - z2*h(i1+inc(3))) 811 f1m = r4*(z1*h(i1+inc(7)) - h(i1+inc(1))-z2*h(i1+inc(4))) 812 f2p = r3ov2*(h(i1+inc(5)) - h(i1+inc(8))) 813 f2m = h(i1+inc(9)) 814 f3p = r1*(h(i1) - z3*h(i1+inc(3))) 815 f3m = r1*(z3*h(i1+inc(4)) - h(i1+inc(1))) 816 817 h(i1 ) = f0 818 h(i1+inc(1)) = f1p 819 h(i1+inc(2)) = f1m 820 h(i1+inc(3)) = f2p 821 h(i1+inc(4)) = f2m 822 h(i1+inc(5)) = f3p 823 h(i1+inc(6)) = f3m 824 825 i1 = i1 + 1 826 end do 827 828 endif 829 830 if (ido9g.eq.1.and.itype.eq.4) then 831c 832c The column side af the matrix: pure G 833c 834 i1 = 1 835 do i=1,14 836 inc(i) = i*jdim 837 end do 838 839 do j=1,jdim 840 841 g0 = 842 & 0.375d0*h(i1) + 0.375d0* h(i1+inc(1)) + 843 & h(i1+inc(2)) - 3.0d0*b2* h(i1+inc(10)) - 844 & 3.0d0*b2*h(i1+inc(11)) + 0.75d0*b2*h(i1+inc(9)) 845 846 g1p = 847 & g1*(4.0d0*b1*h(i1+inc(7)) - 3.0d0*b3* h(i1+inc(13)) - 848 & 3.0d0*b1*h(i1+inc(4))) 849 850 g1m = 851 & g1*(4.0d0*b1*h(i1+inc(8)) - 3.0d0*b3* h(i1+inc(12)) - 852 & 3.0d0*b1*h(i1+inc(6))) 853 854 g2p = 855 & g2*(6.0d0*b2*h(i1+inc(10)) - 6.0d0*b2* h(i1+inc(11)) - 856 & h(i1) + h(i1+inc(1))) 857 858 g2m = 859 & g3*(6.0d0*b3*h(i1+inc(14)) - b1* h(i1+inc(3)) - 860 & b1*h(i1+inc(5))) 861 862 g3p = 863 & g4*(b1* h(i1+inc(4)) - 3.0d0*b3* h(i1+inc(13))) 864 865 g3m = 866 & g4*(3.0d0*b3*h(i1+inc(12)) - b1* h(i1+inc(6))) 867 868 g4p = 869 & g5*( h(i1) - 6.0d0*b2* h(i1+inc(9)) + 870 & h(i1+inc(1))) 871 872 g4m = 873 & g6*b1*( h(i1+inc(3)) - h(i1+inc(5)) ) 874 875 h(i1 ) = g0 876 h(i1+inc(1)) = g1p 877 h(i1+inc(2)) = g1m 878 h(i1+inc(3)) = g2p 879 h(i1+inc(4)) = g2m 880 h(i1+inc(5)) = g3p 881 h(i1+inc(6)) = g3m 882 h(i1+inc(7)) = g4p 883 h(i1+inc(8)) = g4m 884 885 i1 = i1 + 1 886 end do 887 888 endif 889 if ((ido5d.eq.1.and.itype.eq.2) .or. 890 & (ido7f.eq.1.and.itype.eq.3) .or. 891 & (ido9g.eq.1.and.itype.eq.4)) then 892 893 if (debug) print*,'colm.old ',itype,jtype,iend 894 if (itype.eq.2) then 895 iend = 9 896 else if (itype.eq.3) then 897 iend = 17 898 else if (itype.eq.4) then 899 iend = 29 900 endif 901 if (debug) print*,'colm.new ',itype,jtype,iend 902 903 endif 904 905 return 906 end 907 908 subroutine setcon 909 implicit double precision (a-h,o-z) 910 common /intcon/ d3,d5,d7,r1,r2,r3,r4,r3ov2,z1,z2,z3, 911 & g1,g2,g3,g4,g5,g6 912 d3 = dsqrt(3.0d0) 913 d5 = dsqrt(5.0d0) 914 d7 = dsqrt(7.0d0) 915 z1 = 4.0d0/d5 916 z2 = 1.0d0/d5 917 z3 = 3.0d0/d5 918 r1 = 0.5d0*dsqrt(5.0d0/2.0d0) 919 r2 = 1.5d0/d5 920 r3ov2 = 0.5d0*dsqrt(3.0d0) 921 r3 = r3ov2 922 r4 = 0.5d0*dsqrt(3.0d0/2.0d0) 923 924 g1 = dsqrt(5.0d0/8.0d0) 925 g2 = dsqrt(5.0d0/16.0d0) 926 g3 = dsqrt(5.0d0/4.0d0) 927 g4 = dsqrt(35.0d0/8.0d0) 928 g5 = dsqrt(35.0d0/64.0d0) 929 g6 = dsqrt(35.0d0/4.0d0) 930 931 return 932 end 933 934 subroutine dipold(ntt,nor,dmao,focc,focb,vectrs,vectrb,p) 935 implicit double precision (a-h,o-z), integer (i-n) 936 parameter (numatm=2000) 937 parameter (numprm=1600) 938 parameter (numcex=numprm*3) 939 940c dipole integrals 941 942 logical debug 943 common /moldat/ natoms, norbs, nelecs,nat(numatm) 944 common /charge/ dipo(3),ihasq,ihsdp,iqon,idipon 945 common /bounds/ ilow(5,5),iupp(5) 946 common /orbhlp/ mxorb,iuhf,ispd 947 common /coord / xyz(3,numatm) 948 949 integer shella,shelln,shellt,shellc,shladf,aos 950 common/b/exx(numcex), 951 & c1(numcex),c2(numcex),c3(numprm),c4(numprm),c5(numprm), 952 & shladf(numprm),gx(numprm),gy(numprm),gz(numprm), 953 & jan(numprm),shella(numprm),shelln(numprm),shellt(numprm), 954 & shellc(numprm),aos(numprm),nshell,maxtyp 955 956 common /indxe/ jx(35),jy(35),jz(35),ix(20),iy(20),iz(20) 957 common /slagau/ ihasd,isgau,ido5d,ido7f,ido9g,ihasg 958 double precision imprd 959 integer genaos 960 961 dimension focc(*),focb(*),vectrs(*),vectrb(*),nor(*) 962 dimension ca(35),cb(35),sx(36),sy(36),sz(36),dd(3,100), 963 & ccx(330),ccy(330),ccz(330),p(*),dmao(ntt,3), 964 & d12(100),dipon(3) 965 966 lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j)) 967 968c iupp 1 4 10 20 35 969c s 3p 6d 10f 15g 970c 971c ilow 0 1 2 3 4 <- shellt 972c--------------------------------------- 973c 0 1(s) 1 (sp) 1 (spd) 1 1 974c 1 1 2 (p) 5 11 21 975c 2 1 1 5 (d) 11 21 976c 3 1 1 5 11 (f) 21 977c 4 1 1 5 11 21 (g) 978c--------------------------------------- 979c ^ 980c | shellc 981c 982 983 debug = .false. 984 cf = 2.5417463d0 985 idum = genaos(.false.) 986 987 pi = 4.d0*datan(1.d0) 988 989 do i=1,3 990 dipo(i) = 0.d0 991 dipon(i) = 0.d0 992 end do 993 994 do i=1,natoms 995 do j=1,3 996 dipon(j) = dipon(j) + dble(nat(i))*xyz(j,i) 997 end do 998 end do 999 1000 call denini 1001 1002 do ishell=1,nshell 1003 1004 xi = gx(ishell) 1005 yi = gy(ishell) 1006 zi = gz(ishell) 1007 1008 ifrst = shella(ishell) 1009 ilast = ifrst + shelln(ishell) - 1 1010 itype = shellt(ishell) 1011 1012 itype1 = itype + 1 1013 lidim = itype1 + 2 1014 1015 istart = ilow(itype1,shellc(ishell)+1) 1016 iend = iupp(itype1) 1017 iendt = iend 1018 1019 do jshell=1,ishell 1020 1021 xj = gx(jshell) 1022 yj = gy(jshell) 1023 zj = gz(jshell) 1024 1025 jfrst = shella(jshell) 1026 jlast = jfrst + shelln(jshell) - 1 1027 jtype = shellt(jshell) 1028 1029 jtype1 = jtype+1 1030 ljdim = jtype+2 1031 1032 jstart = ilow(jtype1,shellc(jshell)+1) 1033 jend = iupp(jtype1) 1034 jendt = jend 1035 lpdim = lidim + ljdim 1036 1037 xt = xi - xj 1038 yt = yi - yj 1039 zt = zi - zj 1040 1041 r2 = xt*xt + yt*yt + zt*zt 1042 1043 do l=1,100 1044 do k=1,3 1045 dd(k,l) = 0.d0 1046 end do 1047 end do 1048 1049 irp = iendt - istart + 1 1050 jrp = jendt - jstart + 1 1051 lentqp = irp*jrp 1052 1053 if (ishell.eq.jshell) lentqp = (irp*(irp+1))/2 1054 1055 do igauss=ifrst,ilast 1056 1057 ei = exx(igauss) 1058 twoei = ei + ei 1059 1060 call fcij(itype,ifrst,igauss,shladf(ishell),ca) 1061 1062 do jgauss=jfrst,jlast 1063 1064 ej = exx(jgauss) 1065 twoej = ej + ej 1066 1067 call fcij(jtype,jfrst,jgauss,shladf(jshell),cb) 1068 1069 ep = ei + ej 1070 1071 tx = (ei*xi + ej*xj) / ep 1072 ty = (ei*yi + ej*yj) / ep 1073 tz = (ei*zi + ej*zj) / ep 1074 1075 xit = tx - xi 1076 yit = ty - yi 1077 zit = tz - zi 1078 1079 xjt = tx - xj 1080 yjt = ty - yj 1081 zjt = tz - zj 1082 1083 expont = (ej*ei*r2) / ep 1084 pexp = 0.d0 1085 1086 if (expont.lt.600.d0) pexp = dexp(-expont) 1087 1088 call dipint(lidim,ljdim,lpdim, 1089 & xit,xjt,ccx, yit,yjt,ccy, zit,zjt,ccz, 1090 & ei,ej,pi,pexp,sx,sy,sz) 1091 1092 call dipao((ishell.eq.jshell),istart,iend,jstart,jend, 1093 & lidim,ljdim,ca,cb, 1094 & sx,sy,sz,xi,yi,zi,dd) 1095 1096 end do 1097 end do 1098 1099 do l=1,3 1100 call purdf(itype,jtype,istart,jstart,iend,jend,dd(l,1)) 1101 end do 1102 1103 ii = 0 1104 ist = aos(ishell) - 1 1105 jst = aos(jshell) - 1 1106 1107 do i=istart,iend 1108 jnd = jend 1109 if (ishell.eq.jshell) jnd = i 1110 1111 do j=jstart,jnd 1112 ii = ii + 1 1113 d12(ii) = p((ist+i-1)*mxorb + (jst+j)) 1114 if (ist+i.ne.jst+j) d12(ii) = 2.0d0*d12(ii) 1115 ll = lind(ist+i,jst+j) 1116 do l=1,3 1117 dmao(ll,l) = dd(l,ii) 1118 end do 1119 end do 1120 1121 end do 1122 1123 do l=1,3 1124 dipo(l) = dipo(l) - imprd(lentqp,d12,dd,l) 1125 end do 1126 1127 end do 1128 end do 1129 1130 if (debug) then 1131 write(6,'(a)') 'Dipole moment vector (Debye)' 1132 write(6,'(a,3f10.4)') 'elec dipole ',(cf*dipo(i),i=1,3) 1133 write(6,'(a,3f10.4)') 'nuc dipole ',(cf*dipon(i),i=1,3) 1134 1135 dipt = 0.d0 1136 do l=1,3 1137 dipo(l) = dipo(l) + dipon(l) 1138 dipt = dipt + dipo(l)*dipo(l) 1139 end do 1140 dipt = dsqrt(dipt) 1141 write(6,'(a,3f10.4)') 'tot dipole ',(cf*dipo(i),i=1,3) 1142 write(6,'(a,f10.4)') 'total dipole scalar ',cf*dipt 1143 endif 1144 1145 call setnor(nocc,norbs,nor,focc) 1146 call boys(nocc,nor,dmao,vectrs) 1147 1148 if (iuhf.ne.0) then 1149 call setnor(nocc,norbs,nor,focb) 1150 call boys(nocc,nor,dmao,vectrb) 1151 endif 1152 1153 return 1154 end 1155 1156 subroutine dipint(l1max,l2max,l12mxp, 1157 & ax,bx,ccx, ay,by,ccy, az,bz,ccz, 1158 & as,bs,pi,pexp,sx,sy,sz) 1159 implicit double precision (a-h,o-z), integer (i-n) 1160 dimension ccx(l12mxp,l2max,l1max) 1161 dimension ccy(l12mxp,l2max,l1max) 1162 dimension ccz(l12mxp,l2max,l1max) 1163 dimension sx(l2max,l1max),sy(l2max,l1max),sz(l2max,l1max) 1164 dimension s1(21) 1165 1166 do l1=1,l1max 1167 do l2=1,l2max 1168 do l12=1,l12mxp 1169 ccx(l12,l2,l1) = 0.d0 1170 ccy(l12,l2,l1) = 0.d0 1171 ccz(l12,l2,l1) = 0.d0 1172 end do 1173 end do 1174 end do 1175 1176 ccx(2,1,1) = 1.d0 1177 ccy(2,1,1) = 1.d0 1178 ccz(2,1,1) = 1.d0 1179 1180 do l1=1,l1max 1181 do l2=1,l2max 1182 iwmax = l1 + l2 1183 do iw=2,iwmax 1184 1185 if (l1.gt.1) then 1186 1187 ccx(iw,l2,l1) = 1188 & ax*ccx(iw,l2,l1-1) + ccx(iw-1,l2,l1-1) 1189 ccy(iw,l2,l1) = 1190 & ay*ccy(iw,l2,l1-1) + ccy(iw-1,l2,l1-1) 1191 ccz(iw,l2,l1) = 1192 & az*ccz(iw,l2,l1-1) + ccz(iw-1,l2,l1-1) 1193 1194 elseif (l2.gt.1) then 1195 1196 ccx(iw,l2,l1) = 1197 & bx*ccx(iw,l2-1,l1) + ccx(iw-1,l2-1,l1) 1198 ccy(iw,l2,l1) = 1199 & by*ccy(iw,l2-1,l1) + ccy(iw-1,l2-1,l1) 1200 ccz(iw,l2,l1) = 1201 & bz*ccz(iw,l2-1,l1) + ccz(iw-1,l2-1,l1) 1202 1203 endif 1204 1205 end do 1206 end do 1207 end do 1208 1209 temp = 1.d0/(2.d0*(as+bs)) 1210 s1(1) = dsqrt(pi/(as+bs)) 1211 1212 lim = l1max + l2max - 1 1213 1214 if (lim.ge.2) then 1215 1216 fact = 1.d0 1217 1218 do i=2,lim 1219 if (mod(i,2).eq.0) then 1220 s1(i) = 0.d0 1221 else 1222 s1(i) = s1(i-2)*temp*fact 1223 fact = fact + 2.d0 1224 endif 1225 end do 1226 1227 endif 1228 1229 do l1=1,l1max 1230 do l2=1,l2max 1231 1232 x = 0.d0 1233 y = 0.d0 1234 z = 0.d0 1235 iwmax = l1 + l2 - 1 1236 1237 do iw=1,iwmax 1238 x = x + ccx(iw+1,l2,l1)*s1(iw) 1239 y = y + ccy(iw+1,l2,l1)*s1(iw) 1240 z = z + ccz(iw+1,l2,l1)*s1(iw) 1241 end do 1242 1243 sx(l2,l1) = x 1244 sy(l2,l1) = y 1245 sz(l2,l1) = z*pexp 1246 1247 end do 1248 end do 1249 1250 return 1251 end 1252 1253 subroutine dipao(ieqj,istart,iend,jstart,jend, 1254 & l1max,l2max,ca,cb,sx,sy,sz,xa,ya,za,dd) 1255 implicit double precision (a-h,o-z), integer (i-n) 1256 logical ieqj 1257 common /indxe/ lx(35),ly(35),lz(35),kx(20),ky(20),kz(20) 1258 dimension ca(*),cb(*),dd(3,*) 1259 dimension sx(l2max,l1max),sy(l2max,l1max),sz(l2max,l1max) 1260 dimension rxyz(3) 1261 1262 ii = 0 1263 1264 do i=istart,iend 1265 1266 ix = lx(i) 1267 iy = ly(i) 1268 iz = lz(i) 1269 1270 jnd = jend 1271 1272 if (ieqj) jnd = i 1273 1274 do j=jstart,jnd 1275 1276 jx = lx(j) 1277 jy = ly(j) 1278 jz = lz(j) 1279 1280 ii = ii + 1 1281 coef = ca(i)*cb(j) 1282 1283 rxyz(1) = (xa*sx(jx,ix) + sx(jx,ix+1)) 1284 & *sy(jy,iy)*sz(jz,iz) 1285 rxyz(2) = (ya*sy(jy,iy) + sy(jy,iy+1)) 1286 & *sx(jx,ix)*sz(jz,iz) 1287 rxyz(3) = (za*sz(jz,iz) + sz(jz,iz+1)) 1288 & *sx(jx,ix)*sy(jy,iy) 1289 1290 do l=1,3 1291 dd(l,ii) = dd(l,ii) + coef*rxyz(l) 1292 end do 1293 1294 end do 1295 end do 1296 1297 return 1298 end 1299 1300 double precision function imprd(n,a,b,l) 1301 implicit double precision (a-h,o-z), integer (i-n) 1302c improduct of a and b. 1303 dimension a(*), b(3,*) 1304 1305 imprd = 0.d0 1306 1307 if (n.lt.1) return 1308 1309 do i=1,n 1310 imprd = imprd + a(i)*b(l,i) 1311 end do 1312 1313 return 1314 end 1315 1316 subroutine boyd(norb,nor,dmao,cc,nbasis,ntt,cl, 1317 & iord,iir,rij,qpix,qpjx) 1318 implicit double precision (a-h,o-z), integer (i-n) 1319 parameter (numatm=2000) 1320 logical debug 1321 1322C Boys Localization. After QCPE program 354 1323C 1324C NOrb ... Number of orbitals to localize. 1325C CC ... Input and output orbitals (NBasis,). 1326C CL ... Scratch matrix (NBasis,NOrb). 1327 1328 common /orbhlp/ mxorb,iuhf,ispd 1329 common /moldat/ natoms, norbs, nelecs,nat(numatm) 1330 dimension cc(*), cl(nbasis,*), dmao(ntt,3), nor(nbasis,nbasis), 1331 & iord(*), iir(*), rij(ntt,3), qpix(*), qpjx(*) 1332 1333 lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j)) 1334 kind(j,i) = ((i-1)*mxorb + j) 1335 1336c essential that in above definition i and j have swapped places !!! 1337 1338 1339 debug = .false. 1340 cf = 2.5417463d0 1341 1342 do itry=1,3 1343 1344 do i=1,norb 1345 ii = i 1346 1347 do j=1,i 1348 jj = j 1349 sumx = 0.0d0 1350 sumy = 0.0d0 1351 sumz = 0.0d0 1352 1353 do k=1,nbasis 1354 do l=1,nbasis 1355 cpr = cc(kind(k,ii))*cc(kind(l,jj)) 1356 sumx = sumx + cpr*dmao(lind(k,l),1) 1357 sumy = sumy + cpr*dmao(lind(k,l),2) 1358 sumz = sumz + cpr*dmao(lind(k,l),3) 1359 end do 1360 end do 1361 1362 rij(lind(i,j),1) = sumx 1363 rij(lind(i,j),2) = sumy 1364 rij(lind(i,j),3) = sumz 1365 1366 end do 1367 1368 end do 1369 1370 1371 call caltra(norb,ntt,rij,tracia) 1372 1373 if (debug) write(6,'(a,f15.8)') 1374 & ' Initial TraceA=',tracia*cf*cf 1375 1376 do i=1,nbasis 1377 do j=1,nbasis 1378 cl(i,j) = 0.0d0 1379 end do 1380 end do 1381 1382 do i=1,norb 1383 cl(i,i) = 1.0d0 1384 end do 1385 1386 iter = 0 1387 shift = datan(1.0d0) 1388 1389 change = 1.0d0 1390 1391 do while (iter.lt.150.and.change.gt.1.d-6) 1392 1393 if (change.ge.1.d-10) then 1394 change = 0.0d0 1395 1396 iter = iter + 1 1397 1398 do i=1,norb 1399 iir(i) = i 1400 end do 1401 1402 nnn = norb 1403 1404 do i=1,norb 1405 iii = irnumb(nnn) 1406 iord(i) = iir(iii) 1407 iir(iii) = iir(nnn) 1408 nnn = nnn - 1 1409 end do 1410 1411 endif 1412 1413 do iii=1,norb 1414 1415 i = iord(iii) 1416 if (nor(i,i).eq.0) then 1417 1418 li = lind(i,i) 1419 1420 jm = 1 1421 rm = 0.0d0 1422 tm = 0.0d0 1423 sm = 0.0d0 1424 cm = 1.0d0 1425 1426 do j=1,norb 1427 1428 if (i.ne.j.and.nor(i,j).eq.0) then 1429 1430 lj = lind(j,j) 1431 ij = lind(i,j) 1432 t = 0.0d0 1433 tx = 0.0d0 1434 1435 do kk=1,3 1436 t = t + 4.0d0*rij(ij,kk)**2 - rij(li,kk)**2 1437 & - rij(lj,kk)**2 1438 & + 2.0d0*rij(li,kk)*rij(lj,kk) 1439 tx = tx + rij(ij,kk)*(rij(lj,kk)-rij(li,kk)) 1440 end do 1441 1442 if (dabs(t).gt.1.d-10.or.dabs(tx).gt.1.d-10) then 1443 1444 tx = 4.0d0*tx 1445 t = datan2(tx,t)/4.0d0 1446 sign = 1.0d0 1447 if (t.gt.0.0d0) sign = -1.0d0 1448 t = t + sign*shift 1449 itim = 0 1450 1451 s = dsin(t) 1452 c = dcos(t) 1453 rin = 0.0d0 1454 itim = itim + 1 1455 1456 do kk=1,3 1457 qpi = c*c*rij(li,kk) + s*s*rij(lj,kk) 1458 & + 2.0d0*c*s*rij(ij,kk) 1459 qpj = c*c*rij(lj,kk)+s*s*rij(li,kk) 1460 & - 2.0d0*c*s*rij(ij,kk) 1461 rin = rin + qpi*qpi+qpj*qpj - rij(li,kk)**2 1462 & - rij(lj,kk)**2 1463 end do 1464 1465 ttest = dabs(t) - shift 1466 1467 if (dabs(t).le.1.d-8.or.dabs(ttest).lt.1.d-8.or. 1468 & rin.ge.(-1.d-8)) then 1469 if (rin.gt.rm) then 1470 rm = rin 1471 jm = j 1472 sm = s 1473 cm = c 1474 tm = t 1475 endif 1476 else 1477 1478 if (debug) then 1479 write(6,'(a,i4,a,i4,a,f15.8,a,f15.8)') 1480 & ' BoyLoc: No rotation increases integrals for I=', 1481 & i,' J=',j,' Theta=',t,' Change=',rin 1482 endif 1483 1484 endif 1485 1486 endif 1487 1488 endif 1489 1490 end do 1491 1492c end loop over j 1493 1494 rin = rm 1495 j = jm 1496 s = sm 1497 c = cm 1498 t = tm 1499 1500 ij = lind(i,j) 1501 lj = lind(j,j) 1502 1503 if (nor(i,j).eq.0) then 1504 1505 change = change + t*t 1506 1507 do kk=1,3 1508 qpi = c*c*rij(li,kk) + s*s*rij(lj,kk) 1509 & + 2.0d0*c*s*rij(ij,kk) 1510 qpj = c*c*rij(lj,kk) + s*s*rij(li,kk) 1511 & - 2.0d0*c*s*rij(ij,kk) 1512 qpij = (c*c-s*s)*rij(ij,kk) + 1513 & c*s*(rij(lj,kk)-rij(li,kk)) 1514 1515 do k=1,norb 1516 if (i.ne.k.and.j.ne.k) then 1517 ik = lind(i,k) 1518 jk = lind(j,k) 1519 qpix(k) = c*rij(ik,kk) + s*rij(jk,kk) 1520 qpjx(k) = c*rij(jk,kk) - s*rij(ik,kk) 1521 endif 1522 end do 1523 1524 do k=1,norb 1525 if (i.ne.k.and.j.ne.k) then 1526 ik = lind(i,k) 1527 jk = lind(j,k) 1528 rij(ik,kk) = qpix(k) 1529 rij(jk,kk) = qpjx(k) 1530 endif 1531 end do 1532 1533 rin = rin + qpi + qpj - rij(li,kk) - rij(lj,kk) 1534 rij(li,kk) = qpi 1535 rij(lj,kk) = qpj 1536 rij(ij,kk) = qpij 1537 1538 end do 1539 1540 do k=1,norb 1541 c1 = c*cl(k,i) + s*cl(k,j) 1542 c2 = -s*cl(k,i) + c*cl(k,j) 1543 cl(k,i) = c1 1544 cl(k,j) = c2 1545 end do 1546 1547 endif 1548 1549 endif 1550 end do 1551 1552c end loop over iii 1553 1554 change = dsqrt(2.0d0*change/dble(norb*(norb-1))) 1555 1556 1557 if (debug) then 1558 1559 write(6,'(a,i4,a,f15.6)') 1560 & ' BoyLoc: Iteration',iter,' Change=',change 1561 1562 1563 endif 1564 1565 end do 1566 1567 end do 1568 1569 if (iter.gt.150) then 1570 write(6,'(a)') 1571 & ' Localization failed. after 3 tries of 150 iterations.' 1572 else 1573 if (change.le.1.d-6) then 1574 write(6,'(a,i4,a)') 1575 & ' Localization complete after',iter,' iterations.' 1576 endif 1577 endif 1578 1579 do i=1,norb 1580 1581 do k=1,nbasis 1582 qpix(k) = 0.0d0 1583 end do 1584 1585 do j=1,norb 1586 jj = j 1587 1588 do k=1,nbasis 1589 qpix(k) = qpix(k) + cc(kind(k,jj))*cl(j,i) 1590 end do 1591 1592 end do 1593 1594 do k=1,nbasis 1595 cl(k,i) = qpix(k) 1596 end do 1597 1598 end do 1599 1600c put rotated orbitals back 1601 1602 do i=1,norb 1603 ii = i 1604 1605 do k=1,nbasis 1606 cc(kind(k,ii)) = cl(k,i) 1607 end do 1608 1609 end do 1610 1611 if (debug) then 1612 1613 print*,"rotated orbitals" 1614 print*," " 1615 call prev(cc,nbasis,nbasis,mxorb) 1616 1617 call caltra(norb,ntt,rij,tracfa) 1618 write(6,'(a,f15.8)') ' Final TraceA=',tracfa 1619 1620 endif 1621 1622 call caldrv(norb,ntt,rij,dmax,imax,jmax) 1623 1624 if (debug) then 1625 write(6,'(a,i3,a,i3,a,f15.8)') 1626 & ' Largest derivative: DDelta(',imax,',',jmax, 1627 & ')=', dmax 1628 endif 1629 1630 if (dmax.gt.1.d-8) then 1631 write(6,'(a,i3,a,i3,a,f15.8)') 1632 & ' Internal error in BoyLoc, DDelta(',imax,',',jmax, 1633 & ')=',dmax 1634 endif 1635 1636 return 1637 end 1638 1639 integer function irnumb(max) 1640 implicit double precision (a-h,o-z), integer (i-n) 1641 1642 irnumb = int(random()*dble(max)+1.0d0) 1643 1644 return 1645 end 1646 1647 subroutine caltra(norb,ntt,rij,traca) 1648 implicit double precision (a-h,o-z), integer (i-n) 1649 dimension rij(ntt,3) 1650 1651 lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j)) 1652 1653 traca = 0.0d0 1654 1655 ii = 0 1656 do i=1,norb 1657 ii = ii + i 1658 if (i.ne.1) then 1659 traca = traca + rij(ii,1)**2 + rij(ii,2)**2 + rij(ii,3)**2 1660 endif 1661 end do 1662 1663 return 1664 end 1665 1666 subroutine caldrv(norb,ntt,rij,dmax,imax,jmax) 1667 implicit double precision (a-h,o-z), integer (i-n) 1668 dimension rij(ntt,3) 1669 1670 lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j)) 1671 1672 dmax = -1.0d0 1673 1674 do i=2,norb 1675 1676 li = lind(i,i) 1677 1678 do j=1,i-1 1679 1680 lj = lind(j,j) 1681 ij = lind(i,j) 1682 dt = 0.0d0 1683 1684 do kk=1,3 1685 dt = dt + rij(ij,kk)*(rij(li,kk)-rij(lj,kk)) 1686 end do 1687 1688 if (dabs(dt).gt.dmax) then 1689 dmax = dabs(dt) 1690 imax = i 1691 jmax = j 1692 endif 1693 1694 end do 1695 end do 1696 1697 dmax = 2.0d0*dmax 1698 1699 return 1700 end 1701 1702 double precision function random() 1703 implicit double precision (a-h,o-z), integer (i-n) 1704 parameter (mplier=16807,modlus=2147483647,mobymp=127773, 1705 & momdmp=2836) 1706 common /seed/jseed,ifrst,nextn 1707c mseed comes from alloc.c 1708 integer hvlue, lvlue, testv, nextn, mseed 1709 1710 if (ifrst.eq.0) then 1711 jseed = mseed() 1712 nextn = jseed 1713 ifrst = 1 1714 endif 1715 1716 hvlue = nextn / mobymp 1717 lvlue = mod(nextn, mobymp) 1718 testv = mplier*lvlue - momdmp*hvlue 1719 1720 if (testv.gt.0) then 1721 nextn = testv 1722 else 1723 nextn = testv + modlus 1724 endif 1725 1726 random = dble(nextn)/dble(modlus) 1727 1728 return 1729 end 1730 1731 integer function ncoree() 1732 implicit double precision (a-h,o-z), integer (i-n) 1733 parameter (numatm=2000) 1734 common /moldat/ natoms, norbs, nelecs,nat(numatm) 1735 dimension ncore(103) 1736 data ncore /2*0,2*1,6*1,2*5,6*5,2*9,10*9,9,5*14, 1737 & 2*18,10*18,18,5*23,2*27,14*27,10*34, 1738 & 34,5*39,2*43,14*43,50/ 1739 1740 ncoree = 0 1741 1742 do i=1,natoms 1743 ncoree = ncoree + ncore(nat(i)) 1744 end do 1745 1746 return 1747 end 1748 1749 subroutine setnor(nocc,nbasis,nor,focc) 1750 implicit double precision (a-h,o-z), integer (i-n) 1751 dimension focc(*),nor(nbasis,nbasis) 1752 1753 nocc = 0 1754 do i=1,nbasis 1755 if (focc(i).ne.0.d0) then 1756 nocc = nocc + 1 1757 endif 1758 end do 1759 1760 do i=1,nbasis 1761 do j=1,nbasis 1762 nor(i,j) = 1 1763 end do 1764 end do 1765 1766 mcore = ncoree() 1767 do i=mcore+1,nocc 1768 do j=mcore+1,nocc 1769 nor(i,j) = 0 1770 end do 1771 end do 1772 1773 return 1774 end 1775 1776 subroutine pregrd(npts1,npts2,npts3,xden,yden,zden) 1777 implicit double precision(a-h,o-z) 1778 common /plane/ px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat 1779 common /grdhlp/ mx3d,mx3d2 1780 dimension v3(3),xden(*),yden(*),zden(*) 1781 1782 v3(1) = cx 1783 v3(2) = cy 1784 v3(3) = cz 1785 1786 call vnrm(v3) 1787 1788 radius1 = 0.5d0 * r(1) 1789 radius2 = 0.5d0 * r(2) 1790 radius3 = 0.5d0 * r(3) 1791 rf1 = r(1) / dble(npts1-1) 1792 rf2 = r(2) / dble(npts2-1) 1793 rf3 = r(3) / dble(npts3-1) 1794 1795 do kc=1,npts3 1796 z2 = -radius3+dble(kc-1)*rf3 1797 do ic=1,npts1 1798 x2 = -radius1+dble(ic-1)*rf1 1799 do jc=1,npts2 1800 y2 = -radius2+dble(jc-1)*rf2 1801 1802 iadrs = (kc-1)*mx3d2 + ij 1803 1804 xden(iadrs) = v1(1)*x2 + v2(1)*y2 + px - v3(1)*z2 1805 yden(iadrs) = v1(2)*x2 + v2(2)*y2 + py - v3(2)*z2 1806 zden(iadrs) = v1(3)*x2 + v2(3)*y2 + pz - v3(3)*z2 1807 1808 end do 1809 end do 1810 end do 1811 1812 return 1813 end 1814 1815