1c 2c 3c ################################################################# 4c ## COPYRIGHT (C) 2007 by Alexey Kaledin & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################################# 7c 8c ################################################################ 9c ## ## 10c ## program vibbig -- block iterative vibrational analysis ## 11c ## ## 12c ################################################################ 13c 14c 15c "vibbig" performs large-scale vibrational mode analysis using 16c only vector storage and gradient evaluations; preconditioning 17c is via an approximate inverse from a block diagonal Hessian, 18c and a sliding block method is used to converge any number of 19c eigenvectors starting from either lowest or highest frequency 20c 21c literature references: 22c 23c C. Murray, S. C. Racine and E. R. Davidson, "Improved Algorithms 24c for the Lowest Few Eigenvalues and Associated Eigenvectors of 25c Large Matrices", Journal of Computational Physics, 103, 382-389 26c (1992) 27c 28c A. L. Kaledin, "Gradient-Based Direct Normal-Mode Analysis", 29c Journal of Chemical Physics, 122, 184106 (2005) 30c 31c A. L. Kaledin, M. Kaledin and J. M. Bowman, "All-Atom Calculation 32c of the Normal Modes of Bacteriorhodopsin Using a Sliding Block 33c Iterative Diagonalization Method", Journal of Chemical Theory 34c and Computation, 2, 166-174 (2006) 35c 36c 37 program vibbig 38 use atomid 39 use atoms 40 use files 41 use inform 42 use iounit 43 use keys 44 use units 45 use vibs 46 implicit none 47 integer i,j,k,ii,next 48 integer i1,i2,k0,k1,k2 49 integer ivib,ivb1,ivb2 50 integer iblock,iconv 51 integer iter,isave 52 integer nvar,nblk 53 integer nroot,nbasis 54 integer nr,npair 55 integer nlock,nconv 56 integer irange,ifactor 57 integer maxroot,maxiter 58 integer maxhess 59 integer freeunit 60 integer, allocatable :: iblk(:) 61 real*8 fmax,funit 62 real*8 wtol,factor 63 real*8 size,sizmax 64 real*8 space,sum 65 real*8 dfreq,rnorm,rcomp 66 real*8 ratio,shift 67 real*8 uku_min,uku_max 68 real*8, allocatable :: xe(:) 69 real*8, allocatable :: xm(:) 70 real*8, allocatable :: r(:) 71 real*8, allocatable :: rk(:) 72 real*8, allocatable :: hmin(:) 73 real*8, allocatable :: uku(:) 74 real*8, allocatable :: uku0(:) 75 real*8, allocatable :: uu(:) 76 real*8, allocatable :: freq(:) 77 real*8, allocatable :: freqold(:) 78 real*8, allocatable :: tmp1(:) 79 real*8, allocatable :: tmp2(:) 80 real*8, allocatable :: u(:,:) 81 real*8, allocatable :: ur(:,:) 82 real*8, allocatable :: h(:,:) 83 real*8, allocatable :: c(:,:) 84 character*1 answer 85 character*20 keyword 86 character*240 record 87 character*240 string 88 character*240 datafile 89 character*240 blockfile 90 logical exist,restart 91 logical header,done 92c 93c 94c set up the structure and mechanics calculation 95c 96 call initial 97 call getxyz 98 call mechanic 99c 100c set default parameters for the normal mode computation 101c 102 nvar = 3 * n 103 maxroot = 50 104 nr = 6 105 iter = 0 106 isave = 10 107 maxhess = nvar * (nvar-1) / 2 108 maxiter = 100000 109 wtol = 0.00001d0 110 sizmax = 500.0d0 111 header = .true. 112c 113c search the keywords for normal mode parameters 114c 115 do i = 1, nkey 116 next = 1 117 record = keyline(i) 118 call gettext (record,keyword,next) 119 call upcase (keyword) 120 string = record(next:240) 121 if (keyword(1:8) .eq. 'MAXITER ') then 122 read (string,*,err=10,end=10) maxiter 123 else if (keyword(1:11) .eq. 'SAVE-VECTS ') then 124 read (string,*,err=10,end=10) isave 125 else if (keyword(1:10) .eq. 'VIB-ROOTS ') then 126 read (string,*,err=10,end=10) nroot 127 nroot = min(nroot,maxroot) 128 else if (keyword(1:14) .eq. 'VIB-TOLERANCE ') then 129 read (string,*,err=10,end=10) wtol 130 end if 131 10 continue 132 end do 133c 134c find either the lowest or highest normal modes 135c 136 factor = 1.0d0 137 call nextarg (answer,exist) 138 if (.not. exist) then 139 answer = 'L' 140 write (iout,20) answer 141 20 format (/,' Start at Lowest or Highest Frequency', 142 & ' Normal Mode [',a1,'] : ',$) 143 read (input,30) record 144 30 format (a240) 145 next = 1 146 call gettext (record,answer,next) 147 end if 148 call upcase (answer) 149 if (answer .eq. 'H') factor = -1.0d0 150c 151c find cutoff value for desired extreme frequency 152c 153 fmax = -1.0d0 154 call nextarg (string,exist) 155 if (exist) read (string,*,err=40,end=40) fmax 156 40 continue 157 if (fmax .le. 0.0d0) then 158 write (iout,50) 159 50 format (/,' Enter Desired Frequency Cutoff in cm-1', 160 & ' [0.0] : ',$) 161 read (input,60) fmax 162 60 format (f20.0) 163 end if 164 if (fmax .le. 0.0d0) fmax = 0.0d0 165c 166c set default values for some additional parameters 167c 168 funit = factor * efreq * emass 169 ifactor = int(factor) 170 irange = (nvar-nr+1) * max((1-ifactor)/2,0) 171 npair = 2 * nroot 172 nbasis = 3 * nroot 173c 174c open or create eigenvector file for use during restarts 175c 176 ivb1 = freeunit () 177 datafile = filename(1:leng)//'.vb1' 178 call version (datafile,'old') 179 inquire (file=datafile,exist=exist) 180 if (exist) then 181 open (unit=ivb1,file=datafile,status='old',form='unformatted') 182 else 183 open (unit=ivb1,file=datafile,status='new',form='unformatted') 184 end if 185c 186c open or create basis vector file for use during restarts 187c 188 ivb2 = freeunit () 189 datafile = filename(1:leng)//'.vb2' 190 call version (datafile,'old') 191 inquire (file=datafile,exist=exist) 192 if (exist) then 193 restart = .true. 194 open (unit=ivb2,file=datafile,status='old',form='unformatted') 195 else 196 restart = .false. 197 open (unit=ivb2,file=datafile,status='new',form='unformatted') 198 end if 199c 200c perform dynamic allocation of some local arrays 201c 202 allocate (iblk(n)) 203 allocate (xe(nvar)) 204 allocate (xm(nvar)) 205 allocate (r(nvar)) 206 allocate (rk(nvar)) 207 allocate (hmin(nvar)) 208 allocate (uku(nvar)) 209 allocate (uku0(nvar)) 210 allocate (uu(maxhess)) 211 allocate (u(nvar,6)) 212 allocate (ur(nvar,3)) 213c 214c perform dynamic allocation of some global arrays 215c 216 allocate (rho(nvar,nbasis)) 217 allocate (rhok(nvar,nbasis)) 218 allocate (rwork(nvar,nbasis)) 219c 220c store a coordinate vector for each atom 221c 222 do i = 1, n 223 xe(3*i-2) = x(i) / bohr 224 xe(3*i-1) = y(i) / bohr 225 xe(3*i) = z(i) / bohr 226 end do 227c 228c store atomic mass for each coordinate component 229c 230 k = 0 231 do i = 1, n 232 mass(i) = mass(i) * emass 233 do j = 1, 3 234 k = k + 1 235 xm(k) = mass(i) 236 end do 237 end do 238c 239c remove pure translational and rotational modes 240c 241 call trbasis (nvar,nr,xe,u,ur) 242c 243c set number and size of blocks based on storage space 244c 245 space = 0.9d0 * dble(maxhess) 246 do i = 1, n 247 size = 9.0d0 * (dble(n))**2 / dble(i) 248 if (size .lt. space) then 249 nblk = i 250 goto 70 251 end if 252 end do 253 70 continue 254 nblk = max(3,nblk) 255 size = dble(n) / dble(nblk) 256 size = min(size,sizmax) 257 do i = 1, nblk 258 iblk(i) = nint(dble(i)*size) 259 end do 260 do i = nblk, 2, -1 261 iblk(i) = iblk(i) - iblk(i-1) 262 end do 263c 264c get number and size of blocks from an external file 265c 266 iblock = freeunit () 267 blockfile = filename(1:leng)//'.blk' 268 call version (blockfile,'old') 269 inquire (file=blockfile,exist=exist) 270 if (exist) then 271 open (unit=iblock,file=blockfile,status='old') 272 i = 0 273 do while (.true.) 274 i = i + 1 275 read (iblock,*,err=80,end=80) iblk(i) 276 end do 277 80 continue 278 nblk = i - 1 279 close (unit=iblock) 280 end if 281c 282c print info about the atom blocks and preconditioning 283c 284 write (iout,90) 285 90 format (/,' Atom Blocks Used to Subdivide the System :',/) 286 k = 0 287 do i = 1, nblk 288 write (iout,100) i,iblk(i),k+1,k+iblk(i) 289 100 format (' Block :',i7,9x,'Size :',i7,9x,'Atoms :',i7,' to',i7) 290 k = k + iblk(i) 291 end do 292 k = 0 293 do i = 1, nblk 294 k = k + 9*iblk(i)**2 295 end do 296 write (iout,110) k 297 110 format (/,' Storage for Preconditioning Array :',5x,i12) 298c 299c determine number of prior modes available at restart 300c 301 nlock = 0 302 do while (.true.) 303 read (ivb1,err=120,end=120) (r(k),k=1,nvar) 304 nlock = nlock + 1 305 end do 306 120 continue 307 rewind (unit=ivb1) 308 if (nlock .ne. 0) then 309 write (iout,130) nlock 310 130 format (/,' Prior Normal Modes Available at Restart :',i11) 311 end if 312 nconv = nlock 313c 314c compute and diagonalize the Hessian for each block 315c 316 k0 = 0 317 i1 = 1 318 do i = 1, nblk 319 if (i .gt. 1) then 320 k0 = k0 + 9*iblk(i-1)**2 321 i1 = i1 + iblk(i-1) 322 end if 323 i2 = i1 + iblk(i) - 1 324 k1 = 3*i1 - 2 325 k2 = 3*i2 326 call hessblk (mass,k0,i1,i2,uu) 327 call diagblk (k0,k1,3*iblk(i),uu,uku) 328 end do 329c 330c use negative of eigenvalues if doing high frequencies 331c 332 do k = 1, nvar 333 uku(k) = factor * uku(k) 334 uku0(k) = uku(k) 335 end do 336 uku_max = uku(1) 337 uku_min = uku(1) 338 do k = 2, nvar 339 if (uku(k) .gt. uku_max) uku_max = uku(k) 340 if (uku(k) .lt. uku_min) uku_min = uku(k) 341 end do 342c 343c perform dynamic allocation of some local arrays 344c 345 allocate (freq(nbasis)) 346 allocate (freqold(nbasis)) 347 allocate (tmp1(nbasis)) 348 allocate (tmp2(nbasis)) 349 allocate (h(nbasis,nbasis)) 350 allocate (c(nbasis,nbasis)) 351c 352c if restarting, read trial vectors and estimate eigenvalues 353c 354 if (restart) then 355 do i = 1, npair 356 read (ivb2) (rho(k,i),k=1,nvar) 357 read (ivb2) (rhok(k,i),k=1,nvar) 358 end do 359 do i = 1, nroot 360 h(i,i) = 0.0d0 361 do k = 1, nvar 362 h(i,i) = h(i,i) + rhok(k,i)*rho(k,i) 363 end do 364 freqold(i) = sign(1.0d0,h(i,i)) * sqrt(abs(h(i,i))) 365 end do 366 goto 140 367 end if 368c 369c if not restarting, generate initial guess eigenvectors 370c 371 do i = 1, nroot 372 call trigger (nvar,nbasis,nr,ifactor,nblk,iblk,u,uu,r) 373 do k = 1, nvar 374 rho(k,i) = r(k) 375 end do 376 end do 377c 378c project out locked roots from components of rho 379c 380 call project (nvar,nconv,ivb1,nroot,0) 381 call projectk (nvar,nconv,ivb1,nroot,0) 382c 383c reload and make vector orthonormal to existing basis 384c 385 do i = 1, nroot 386 do k = 1, nvar 387 r(k) = rho(k,i) 388 end do 389 if (i .eq. 1) then 390 sum = 0.0d0 391 do k = 1, nvar 392 sum = sum + r(k)*r(k) 393 end do 394 sum = sqrt(sum) 395 do k = 1, nvar 396 r(k) = r(k) / sum 397 end do 398 else 399 call gsort (nvar,i-1,r) 400 end if 401 do k = 1, nvar 402 rho(k,i) = r(k) 403 end do 404 end do 405c 406c store K on rho 407c 408 do i = 1, nroot 409 do k = 1, nvar 410 r(k) = rho(k,i) 411 end do 412 call konvec (nvar,xm,xe,r,rk) 413 do k = 1, nvar 414 rhok(k,i) = factor * rk(k) 415 end do 416 end do 417c 418c make nroot-by-nroot CI matrix 419c 420 do i = 1, nroot 421 do j = i, nroot 422 h(i,j) = 0.0d0 423 do k = 1, nvar 424 h(i,j) = h(i,j) + rhok(k,i)*rho(k,j) 425 end do 426 h(j,i) = h(i,j) 427 end do 428 end do 429c 430c diagonalize and use first nroot solutions as starting basis 431c 432 call transform (nroot,nbasis,h,c) 433c 434c fill up arrays 435c 436 do k = 1, nvar 437 do j = 1, nroot 438 tmp1(j) = 0.0d0 439 tmp2(j) = 0.0d0 440 do i = 1, nroot 441 tmp1(j) = tmp1(j) + c(i,j)*rho(k,i) 442 tmp2(j) = tmp2(j) + c(i,j)*rhok(k,i) 443 end do 444 end do 445 do j = 1, nroot 446 rho(k,j) = tmp1(j) 447 rhok(k,j) = tmp2(j) 448 end do 449 end do 450c 451c residues of guesses 452c 453 do i = 1, nroot 454 freq(i) = funit * sign(1.0d0,h(i,i)) * sqrt(abs(h(i,i))) 455 freqold(i) = freq(i) 456 do k = 1, nvar 457 rk(k) = rhok(k,i) - h(i,i)*rho(k,i) 458 end do 459c 460c use Davidson preconditioner if finding low frequencies 461c 462 if (factor .gt. 0.0d0) then 463 call preconblk (nvar,nblk,iblk,uku,uu,h(i,i),hmin(i),rk) 464 end if 465c 466c project residual onto P-space 467c 468 call qonvec (nvar,nr,u,rk,r) 469 do k = 1, nvar 470 rho(k,i+nroot) = r(k) 471 end do 472 end do 473c 474c project out locked roots from components of rho 475c 476 call project (nvar,nconv,ivb1,nroot,nroot) 477 call projectk (nvar,nconv,ivb1,nroot,nroot) 478c 479c reload and make vector orthonormal to existing basis 480c 481 do i = 1, nroot 482 do k = 1, nvar 483 r(k) = rho(k,i+nroot) 484 end do 485 call gsort (nvar,nroot+i-1,r) 486 do k = 1, nvar 487 rho(k,i+nroot) = r(k) 488 end do 489 end do 490c 491c store K on rho 492c 493 do i = 1, nroot 494 do k = 1, nvar 495 r(k) = rho(k,i+nroot) 496 end do 497 call konvec (nvar,xm,xe,r,rk) 498 do k = 1, nvar 499 rhok(k,i+nroot) = factor * rk(k) 500 end do 501 end do 502c 503c make npair-by-npair CI matrix 504c 505 do i = 1, npair 506 do j = i, npair 507 h(i,j) = 0.0d0 508 do k = 1, nvar 509 h(i,j) = h(i,j) + rhok(k,i)*rho(k,j) 510 end do 511 h(j,i) = h(i,j) 512 end do 513 end do 514c 515c diagonalize and use first nroot solutions as new guess 516c 517 call transform (npair,nbasis,h,c) 518 do k = 1, nvar 519 do j = 1, nroot 520 tmp1(j) = 0.0d0 521 tmp2(j) = 0.0d0 522 do i = 1, npair 523 tmp1(j) = tmp1(j) + c(i,j)*rho(k,i) 524 tmp2(j) = tmp2(j) + c(i,j)*rhok(k,i) 525 end do 526 end do 527c 528c old solution fills up 2nd block 529c 530 do j = 1, nroot 531 rho(k,j+nroot) = rho(k,j) 532 rhok(k,j+nroot) = rhok(k,j) 533 end do 534c 535c new solution fills up 1st block 536c 537 do j = 1, nroot 538 rho(k,j) = tmp1(j) 539 rhok(k,j) = tmp2(j) 540 end do 541c 542c orthogonalize 2nd block to 1st 543c 544 do j = 1, nroot 545 do i = 1, nroot 546 rho(k,j+nroot) = rho(k,j+nroot) - c(j,i)*rho(k,i) 547 rhok(k,j+nroot) = rhok(k,j+nroot) - c(j,i)*rhok(k,i) 548 end do 549 end do 550 end do 551c 552c orthogonalize 2nd block on itself 553c 554 sum = 0.0d0 555 do k = 1, nvar 556 sum = sum + rho(k,nroot+1)*rho(k,nroot+1) 557 end do 558 sum = sqrt(sum) 559c 560c normalize leading vector 561c 562 do k = 1, nvar 563 rho(k,nroot+1) = rho(k,nroot+1) / sum 564 rhok(k,nroot+1) = rhok(k,nroot+1) / sum 565 end do 566c 567c orthogonalize the rest one-by-one 568c 569 if (nroot .gt. 1) then 570 do i = 2, max(2,nroot) 571 do j = 1, i-1 572 sum = 0.0d0 573 do k = 1, nvar 574 sum = sum + rho(k,i+nroot)*rho(k,j+nroot) 575 end do 576 do k = 1, nvar 577 rho(k,i+nroot) = rho(k,i+nroot)-sum*rho(k,j+nroot) 578 rhok(k,i+nroot) = rhok(k,i+nroot)-sum*rhok(k,j+nroot) 579 end do 580 end do 581 sum = 0.0d0 582 do k = 1, nvar 583 sum = sum + rho(k,i+nroot)*rho(k,i+nroot) 584 end do 585 sum = sqrt(sum) 586 do k = 1, nvar 587 rho(k,i+nroot) = rho(k,i+nroot) / sum 588 rhok(k,i+nroot) = rhok(k,i+nroot) / sum 589 end do 590 end do 591 end if 592c 593c residue of new solution (if restarting, begin here) 594c 595 140 continue 596 do i = 1, nroot 597 freq(i) = funit * sign(1.0d0,h(i,i)) * sqrt(abs(h(i,i))) 598 freq(i+nroot) = funit * sign(1.0d0,h(i+nroot,i+nroot)) 599 & * sqrt(abs(h(i+nroot,i+nroot))) 600 freq(i+npair) = funit * sign(1.0d0,h(i+npair,i+npair)) 601 & * sqrt(abs(h(i+npair,i+npair))) 602 do k = 1, nvar 603 rk(k) = rhok(k,i) - h(i,i)*rho(k,i) 604 end do 605c 606c use Davidson preconditioner if finding low frequencies 607c 608 if (factor .gt. 0.0d0) then 609 call preconblk (nvar,nblk,iblk,uku,uu,h(i,i),hmin(i),rk) 610 end if 611c 612c project onto P-space 613c 614 call qonvec (nvar,nr,u,rk,r) 615 do k = 1, nvar 616 rho(k,i+npair) = r(k) 617 end do 618 end do 619c 620c project out locked roots from components of rho 621c 622 call project (nvar,nconv,ivb1,nroot,npair) 623 call projectk (nvar,nconv,ivb1,nroot,npair) 624c 625c reload and orthogonalize to 1st + 2nd 626c 627 do i = 1, nroot 628 do k = 1, nvar 629 r(k) = rho(k,i+npair) 630 end do 631 call gsort (nvar,npair+i-1,r) 632 do k = 1, nvar 633 rho(k,i+npair) = r(k) 634 end do 635 end do 636c 637c store K on rho 638c 639 do i = 1, nroot 640 do k = 1, nvar 641 r(k) = rho(k,i+npair) 642 end do 643 call konvec (nvar,xm,xe,r,rk) 644 do k = 1, nvar 645 rhok(k,i+npair) = factor * rk(k) 646 end do 647 end do 648c 649c beginning of iterations 650c 651 iconv = 0 652 150 continue 653 done = .false. 654 iter = iter + 1 655c 656c make nbasis-by-nbasis CI matrix 657c 658 do i = 1, nbasis 659 do j = i, nbasis 660 h(i,j) = 0.0d0 661 do k = 1, nvar 662 h(i,j) = h(i,j) + rhok(k,i)*rho(k,j) 663 end do 664 h(j,i) = h(i,j) 665 end do 666 end do 667c 668c list of previous frequencies 669c 670 do i = 1, npair 671 freqold(i) = freq(i) 672 end do 673c 674c diagonalize and use first nroot solutions as new guess 675c 676 call transform (nbasis,nbasis,h,c) 677c 678c check for collapse based on leading component of ground state 679c 680 if (iconv.eq.0 .and. nconv.gt.0) then 681 sum = sqrt(1.0d0-c(1,1)**2) 682 if (sum .gt. 0.9d0) then 683 write (iout,160) nconv-nlock 684 160 format (/,' Number of Converged Normal Modes :',6x,i12) 685 write (iout,170) 686 170 format (/,' VIBBIG -- Loss of Root Identity; Please', 687 & ' Try to Restart') 688 close (unit=ivb2,status='delete') 689 goto 270 690 end if 691 end if 692c 693c list of new frequencies 694c 695 do i = 1, npair 696 freq(i) = funit * sign(1.0d0,h(i,i)) * sqrt(abs(h(i,i))) 697 end do 698c 699c check if first few have converged 700c 701 iconv = 0 702 180 continue 703 dfreq = freqold(iconv+1) - freq(iconv+1) 704 if (dfreq*factor.gt.0.0d0 .and. dfreq*factor.lt.wtol) then 705 iconv = iconv + 1 706 goto 180 707 end if 708c 709c shift levels of preconditioner matrix; since the Hessian 710c is gradually deflated, reduce effect of the preconditioner 711c based on a simple 1/x curve, the uku levels are squeezed 712c upwards to eventually lead to a unit operator 713c 714 if (iconv .gt. 0) then 715 ratio = dble(nconv+iconv) / dble(nvar) 716 shift = uku_min / (1.0d0-ratio) 717 shift = shift + h(iconv+nroot,iconv+nroot) 718c 719c do a regular shift, which also seems to work 720c 721 do k = 1, nvar 722 uku(k) = uku_max + (uku0(k)-uku_max)*(uku_max-shift) 723 & / (uku_max-uku_min) 724 end do 725c 726c move cursor to end of storage file 727c 728 do i = 1, nconv 729 read (ivb1) (rk(k),k=1,nvar) 730 end do 731c 732c norm of residual 733c 734 do j = 1, iconv 735 rnorm = 0.0d0 736 do k = 1, nvar 737 r(k) = 0.0d0 738 rk(k) = 0.0d0 739 do i = 1, nbasis 740 r(k) = r(k)+c(i,j)*rho(k,i) 741 rk(k) = rk(k)+c(i,j)*rhok(k,i) 742 end do 743 rnorm = rnorm + (rk(k)-h(j,j)*r(k))**2 744 end do 745 rnorm = sqrt(rnorm) 746c 747c component of root in R-space 748c 749 do i = 1, 3 750 tmp1(i) = 0.0d0 751 do k = 1, nvar 752 tmp1(i) = tmp1(i) + ur(k,i)*r(k) 753 end do 754 end do 755 rcomp = 0.0d0 756 do k = 1, nvar 757 sum = 0.0d0 758 do i = 1, 3 759 sum = sum + ur(k,i)*tmp1(i) 760 end do 761 rcomp = rcomp + sum*sum 762 end do 763 rcomp = sqrt(rcomp) 764c 765c write the converged mode to formatted and binary files 766c 767 ivib = irange + ifactor*(nconv+j) 768 if ((header.or.verbose) .and. j.eq.1) then 769 header = .false. 770 write (iout,190) 771 190 format (/,' Converged Normal Modes from Iterative', 772 & ' Vibrational Analysis :') 773 write (iout,200) 774 200 format (/,4x,'Mode',7x,'Frequency',8x,'Delta',10x, 775 & 'R Norm',10x,'Orthog') 776 if (.not. verbose) then 777 write (iout,210) 778 210 format () 779 end if 780 end if 781 dfreq = freqold(j) - freq(j) 782 write (iout,220) ivib,freq(j),dfreq,rnorm,rcomp 783 220 format (i8,f15.3,3d16.4) 784 call prtvib (ivib,r) 785 write (ivb1) (r(k),k=1,nvar) 786 end do 787 rewind (unit=ivb1) 788c 789c update total number of vectors locked on disk 790c 791 nconv = nconv + iconv 792 if (freq(iconv)*factor .ge. fmax*factor) then 793 done = .true. 794 close (unit=ivb1) 795 end if 796 end if 797c 798c shift frequency arrays by iconv 799c 800 do i = 1, npair 801 freq(i) = freq(i+iconv) 802 freqold(i) = freqold(i+iconv) 803 end do 804 do k = 1, nvar 805 do j = 1, nroot+iconv 806 tmp1(j) = 0.0d0 807 tmp2(j) = 0.0d0 808 do i = 1, nbasis 809 tmp1(j) = tmp1(j) + c(i,j)*rho(k,i) 810 tmp2(j) = tmp2(j) + c(i,j)*rhok(k,i) 811 end do 812 end do 813c 814c old solution fills up 2nd block 815c 816 do j = 1, nroot 817 rho(k,j+nroot+iconv) = rho(k,j+iconv) 818 rhok(k,j+nroot+iconv) = rhok(k,j+iconv) 819 end do 820c 821c new solution fills up 1st block 822c 823 do j = 1, nroot 824 rho(k,j+iconv) = tmp1(j+iconv) 825 rhok(k,j+iconv) = tmp2(j+iconv) 826 end do 827c 828c shift index down by iconv 829c 830 do j = 1, npair 831 rho(k,j) = rho(k,j+iconv) 832 rhok(k,j) = rhok(k,j+iconv) 833 end do 834c 835c orthogonalize 2nd block to 1st + iconv roots 836c 837 do j = 1, nroot 838 do i = 1, nroot 839 rho(k,j+nroot) = rho(k,j+nroot) 840 & - c(j+iconv,i+iconv)*rho(k,i) 841 rhok(k,j+nroot) = rhok(k,j+nroot) 842 & - c(j+iconv,i+iconv)*rhok(k,i) 843 end do 844 do i = 1, iconv 845 rho(k,j+nroot) = rho(k,j+nroot) - c(j+iconv,i)*tmp1(i) 846 rhok(k,j+nroot) = rhok(k,j+nroot) - c(j+iconv,i)*tmp2(i) 847 end do 848 end do 849 end do 850c 851c orthogonalize 2nd block on itself 852c 853 sum = 0.0d0 854 do k = 1, nvar 855 sum = sum + rho(k,nroot+1)*rho(k,nroot+1) 856 end do 857 sum = sqrt(sum) 858c 859c normalize leading vector 860c 861 do k = 1, nvar 862 rho(k,nroot+1) = rho(k,nroot+1) / sum 863 rhok(k,nroot+1) = rhok(k,nroot+1) / sum 864 end do 865c 866c orthogonalize the rest one-by-one 867c 868 if (nroot .gt. 1) then 869 do i = 2, max(2,nroot) 870 do j = 1, i-1 871 sum = 0.0d0 872 do k = 1, nvar 873 sum = sum + rho(k,i+nroot)*rho(k,j+nroot) 874 end do 875 do k = 1, nvar 876 rho(k,i+nroot) = rho(k,i+nroot)-sum*rho(k,j+nroot) 877 rhok(k,i+nroot) = rhok(k,i+nroot)-sum*rhok(k,j+nroot) 878 end do 879 end do 880 sum = 0.0d0 881 do k = 1, nvar 882 sum = sum + rho(k,i+nroot)*rho(k,i+nroot) 883 end do 884 sum = sqrt(sum) 885 do k = 1, nvar 886 rho(k,i+nroot) = rho(k,i+nroot) / sum 887 rhok(k,i+nroot) = rhok(k,i+nroot) / sum 888 end do 889 end do 890 end if 891c 892c print a header for the current iteration 893c 894 if (verbose) then 895 write (iout,230) iter,iconv,nconv 896 230 format (/,' Iteration',i7,11x,'New Modes',i6,10x, 897 & ' Total Modes',i6,/) 898 write (iout,240) 899 240 format (4x,'Mode',7x,'Frequency',8x,'Delta',10x, 900 & 'R Norm',10x,'Orthog') 901 end if 902c 903c norm of residual 904c 905 do i = 1, nroot 906 rnorm = 0.0d0 907 do k = 1, nvar 908 rnorm = rnorm + (rhok(k,i)-h(i+iconv,i+iconv)*rho(k,i))**2 909 end do 910 rnorm = sqrt(rnorm) 911c 912c calculate root's component in R-space 913c 914 do j = 1, 3 915 tmp1(j) = 0.0d0 916 do k = 1, nvar 917 tmp1(j) = tmp1(j) + ur(k,j)*rho(k,i) 918 end do 919 end do 920 rcomp = 0.0d0 921 do k = 1, nvar 922 sum = 0.0d0 923 do j = 1, 3 924 sum = sum + ur(k,j)*tmp1(j) 925 end do 926 rcomp = rcomp + sum*sum 927 end do 928 rcomp = sqrt(rcomp) 929 dfreq = freqold(i) - freq(i) 930 if (verbose) then 931 write (iout,250) irange+ifactor*(i+nconv), 932 & freq(i),dfreq,rnorm,rcomp 933 250 format (i8,f15.3,3d16.4) 934 end if 935 end do 936c 937c save vectors needed to restart a calculation 938c 939 if (mod(iter,isave) .eq. 0) then 940 rewind (unit=ivb2) 941 do i = 1, npair 942 write (ivb2) (rho(k,i),k=1,nvar) 943 write (ivb2) (rhok(k,i),k=1,nvar) 944 end do 945 end if 946c 947c prepare restart if finished or iterations exhausted 948c 949 if (done .or. iter.eq.maxiter) then 950 write (iout,260) nconv-nlock 951 260 format (/,' Number of Converged Normal Modes :',6x,i12) 952 rewind (ivb2) 953 do i = 1, npair 954 write (ivb2) (rho(k,i),k=1,nvar) 955 write (ivb2) (rhok(k,i),k=1,nvar) 956 end do 957 close (unit=ivb2) 958 goto 270 959 end if 960c 961c as above, make sure no prior roots are mixed in the basis 962c 963 do i = 1, npair 964 do k = 1, nvar 965 r(k) = rho(k,i) 966 end do 967 call qonvec (nvar,nr,u,r,rk) 968 do k = 1, nvar 969 rho(k,i) = rk(k) 970 end do 971 do k = 1, nvar 972 r(k) = rhok(k,i) 973 end do 974 call qonvec (nvar,nr,u,r,rk) 975 do k = 1, nvar 976 rhok(k,i) = rk(k) 977 end do 978 end do 979c 980c project out locked roots from components of rho 981c 982 call project (nvar,nconv,ivb1,npair,0) 983 call projectk (nvar,nconv,ivb1,npair,0) 984c 985c setup next iteration; solution residue, Davidson weight 986c 987 do i = 1, nroot 988 do k = 1, nvar 989 rk(k) = rhok(k,i) - h(i+iconv,i+iconv)*rho(k,i) 990 end do 991c 992c use Davidson preconditioner if finding low frequencies 993c 994 ii = i + iconv 995 if (factor .gt. 0.0d0) then 996 call preconblk (nvar,nblk,iblk,uku,uu,h(ii,ii),hmin(i),rk) 997 end if 998c 999c project residual onto R-space 1000c 1001 call qonvec (nvar,nr,u,rk,r) 1002 do k = 1, nvar 1003 rho(k,i+npair) = r(k) 1004 end do 1005 end do 1006c 1007c project out locked roots from components of rho 1008c 1009 call project (nvar,nconv,ivb1,nroot,npair) 1010c 1011c reload and orthogonalize to 1st + 2nd 1012c 1013 do i = 1, nroot 1014 do k = 1, nvar 1015 r(k) = rho(k,i+npair) 1016 end do 1017 call gsort (nvar,npair+i-1,r) 1018 do k = 1, nvar 1019 rho(k,i+npair) = r(k) 1020 end do 1021 end do 1022c 1023c store K on rho 1024c 1025 do i= 1, nroot 1026 do k = 1, nvar 1027 r(k) = rho(k,i+npair) 1028 end do 1029 call konvec (nvar,xm,xe,r,rk) 1030 call qonvec(nvar,nr,u,rk,r) 1031 do k = 1, nvar 1032 rhok(k,i+npair) = factor * r(k) 1033 end do 1034 end do 1035c 1036c project out locked roots from components of rhok 1037c 1038 call projectk (nvar,nconv,ivb1,nroot,npair) 1039 goto 150 1040 270 continue 1041c 1042c perform deallocation of some local arrays 1043c 1044 deallocate (iblk) 1045 deallocate (xe) 1046 deallocate (xm) 1047 deallocate (r) 1048 deallocate (rk) 1049 deallocate (hmin) 1050 deallocate (uku) 1051 deallocate (uku0) 1052 deallocate (uu) 1053 deallocate (u) 1054 deallocate (ur) 1055 deallocate (freq) 1056 deallocate (freqold) 1057 deallocate (tmp1) 1058 deallocate (tmp2) 1059 deallocate (h) 1060 deallocate (c) 1061c 1062c perform any final tasks before program exit 1063c 1064 call final 1065 end 1066c 1067c 1068c ############################################################## 1069c ## ## 1070c ## subroutine trigger -- get initial trial eigenvectors ## 1071c ## ## 1072c ############################################################## 1073c 1074c 1075c "trigger" constructs a set of initial trial vectors for 1076c use during sliding block iterative matrix diagonalization 1077c 1078c 1079 subroutine trigger (nvar,nbasis,nr,ifactor,nblk,iblk,u,uu,r) 1080 implicit none 1081 integer i,j,k,m 1082 integer k0,k1,k2 1083 integer nvar,nbasis 1084 integer nr,ifactor 1085 integer nblk,nguess 1086 integer iblk(*) 1087 real*8 w,sum 1088 real*8 random 1089 real*8 r(*) 1090 real*8 uu(*) 1091 real*8 u(nvar,*) 1092 real*8, allocatable :: tmp(:) 1093 external random 1094c 1095c 1096c set the number of random guesses 1097c 1098 nguess = 1 + int(dble(nbasis)/dble(nblk)) 1099c 1100c zero out the trial vector 1101c 1102 do k = 1, nvar 1103 r(k) = 0.0d0 1104 end do 1105c 1106c create overlap with the entire P-space 1107c 1108 k0 = 0 1109 k1 = 1 1110 do i = 1, nblk 1111 if (i .gt. 1) then 1112 k0 = k0 + 9*iblk(i-1)**2 1113 k1 = k1 + 3*iblk(i-1) 1114 end if 1115 k2 = k1 + 3*iblk(i) - 1 1116c 1117c scan over rows of the Hessian 1118c 1119 m = 0 1120 do j = 1, 3*iblk(i) 1121 if (ifactor .eq. 1) then 1122 if (j .gt. min(nguess,3*iblk(i))) then 1123 w = 0.0d0 1124 else 1125 w = random() - 0.5d0 1126 end if 1127 else 1128 if (j .lt. (3*iblk(i)-min(nguess,3*iblk(i))+1)) then 1129 w = 0.0d0 1130 else 1131 w = random() - 0.5d0 1132 end if 1133 end if 1134 do k = k1, k2 1135 m = m + 1 1136 r(k) = r(k) + w*uu(k0+m) 1137 end do 1138 end do 1139 end do 1140c 1141c perform dynamic allocation of some local arrays 1142c 1143 allocate (tmp(nvar)) 1144c 1145c project the vector onto P-space 1146c 1147 call qonvec (nvar,nr,u,r,tmp) 1148c 1149c perform a normalization 1150c 1151 sum = 0.0d0 1152 do i = 1, nvar 1153 sum = sum + tmp(i)**2 1154 end do 1155 sum = sqrt(sum) 1156 do i = 1, nvar 1157 r(i) = tmp(i) / sum 1158 end do 1159c 1160c perform deallocation of some local arrays 1161c 1162 deallocate (tmp) 1163 return 1164 end 1165c 1166c 1167c ################################################################ 1168c ## ## 1169c ## subroutine trbasis -- set translation/rotation vectors ## 1170c ## ## 1171c ################################################################ 1172c 1173c 1174c "trbasis" forms translation and rotation basis vectors used 1175c during vibrational analysis via block iterative diagonalization 1176c 1177c 1178 subroutine trbasis (nvar,nr,xe,u,ur) 1179 use atomid 1180 use atoms 1181 implicit none 1182 integer i,j,k 1183 integer nvar,nr 1184 real*8 tmass,sum 1185 real*8 ra,rha,pr 1186 real*8 cm(3) 1187 real*8 r(3) 1188 real*8 e(3,3) 1189 real*8 c(3,3) 1190 real*8 xe(*) 1191 real*8 u(nvar,*) 1192 real*8 ur(nvar,*) 1193c 1194c 1195c zero out the translation and rotation vectors 1196c 1197 do i = 1, 6 1198 do j = 1, nvar 1199 u(j,i) = 0.0d0 1200 end do 1201 end do 1202c 1203c get the total mass of the system 1204c 1205 tmass = 0.0d0 1206 do i = 1, n 1207 tmass = tmass + mass(i) 1208 end do 1209c 1210c set basis vectors for translations 1211c 1212 do i = 1, n 1213 u(3*i-2,1) = sqrt(mass(i)/tmass) 1214 u(3*i-1,1) = 0.0d0 1215 u(3*i,1) = 0.0d0 1216 u(3*i-2,2) = 0.0d0 1217 u(3*i-1,2) = sqrt(mass(i)/tmass) 1218 u(3*i,2) = 0.0d0 1219 u(3*i-2,3) = 0.0d0 1220 u(3*i-1,3) = 0.0d0 1221 u(3*i,3) = sqrt(mass(i)/tmass) 1222 end do 1223c 1224c move center of mass to origin 1225c 1226 do i = 1, 3 1227 cm(i) = 0.0d0 1228 end do 1229 do i = 1, n 1230 do j = 1, 3 1231 cm(j) = cm(j) + xe(3*(i-1)+j)*mass(i) 1232 end do 1233 end do 1234 do i = 1, n 1235 do j = 1, 3 1236 xe(3*(i-1)+j) = xe(3*(i-1)+j) - cm(j)/tmass 1237 end do 1238 end do 1239c 1240c get the moments of inertia 1241c 1242 do i = 1, 3 1243 e(i,i) = 0.0d0 1244 end do 1245 do i = 1, n 1246 e(1,1) = e(1,1) + ((xe(3*i-1)**2+xe(3*i)**2))*mass(i) 1247 e(2,2) = e(2,2) + ((xe(3*i-2)**2+xe(3*i)**2))*mass(i) 1248 e(3,3) = e(3,3) + ((xe(3*i-2)**2+xe(3*i-1)**2))*mass(i) 1249 end do 1250 do i = 1, 2 1251 do j = i+1, 3 1252 e(i,j) = 0.0d0 1253 do k = 1, n 1254 e(i,j) = e(i,j) - xe(3*(k-1)+i)*xe(3*(k-1)+j)*mass(k) 1255 end do 1256 e(j,i) = e(i,j) 1257 end do 1258 end do 1259c 1260c diagonalize to get principal axes 1261c 1262 call jacobi (3,e,cm,c) 1263c 1264c construction of principle rotations 1265c 1266 do i = 1, 3 1267 do j = 1, n 1268 ra = 0.0d0 1269 pr = 0.0d0 1270 do k = 1, 3 1271 cm(k) = xe(3*(j-1)+k) 1272 ra = ra + cm(k)**2 1273 pr = pr + cm(k)*c(k,i) 1274 end do 1275 rha = sqrt(ra-pr**2) 1276 r(1) = c(2,i)*cm(3) - c(3,i)*cm(2) 1277 r(2) = c(3,i)*cm(1) - c(1,i)*cm(3) 1278 r(3) = c(1,i)*cm(2) - c(2,i)*cm(1) 1279 sum = 0.0d0 1280 do k = 1, 3 1281 sum = sum + r(k)**2 1282 end do 1283 sum = sqrt(sum) 1284 do k = 1, 3 1285 ur(3*(j-1)+k,i) = sqrt(mass(j)) * rha*r(k)/sum 1286 end do 1287 end do 1288 sum = 0.0d0 1289 do j = 1, nvar 1290 sum = sum + ur(j,i)**2 1291 end do 1292 sum = sqrt(sum) 1293 do j = 1, nvar 1294 ur(j,i) = ur(j,i) / sum 1295 end do 1296 end do 1297c 1298c set basis vectors for rotation 1299c 1300 if (nr .eq. 6) then 1301 do i = 1, 3 1302 do j = 1, nvar 1303 u(j,i+3) = ur(j,i) 1304 end do 1305 end do 1306 end if 1307 return 1308 end 1309c 1310c 1311c ################################################################# 1312c ## ## 1313c ## subroutine preconblk -- precondition atom block Hessian ## 1314c ## ## 1315c ################################################################# 1316c 1317c 1318c "preconblk" applies a preconditioner to an atom block section 1319c of the Hessian matrix 1320c 1321c 1322 subroutine preconblk (nvar,nblk,iblk,uku,uu,h,hmin,rk) 1323 implicit none 1324 integer i,j,k,l 1325 integer nvar,nblk 1326 integer k0,k1,k2,l2 1327 integer iblk(*) 1328 real*8 h,hmin 1329 real*8 uku(*) 1330 real*8 rk(*) 1331 real*8 uu(*) 1332 real*8, allocatable :: d(:) 1333 real*8, allocatable :: work(:) 1334c 1335c 1336c find smallest element of |h-uku| 1337c 1338 hmin = abs(h-uku(1)) 1339 do k = 2, nvar 1340 if (abs(h-uku(k)) .lt. hmin) then 1341 hmin = abs(h-uku(k)) 1342 end if 1343 end do 1344c 1345c perform dynamic allocation of some local arrays 1346c 1347 allocate (d(nvar)) 1348 allocate (work(nvar)) 1349c 1350c assign values to temporary array 1351c 1352 do k = 1, nvar 1353 d(k) = h - uku(k) 1354 end do 1355c 1356c invert array via d=hmin/d, where hmin=min{|d(k)|} 1357c 1358 do k = 1, nvar 1359 d(k) = hmin / d(k) 1360 end do 1361c 1362c create overlap with the entire rk array 1363c 1364 k0 = 0 1365 k1 = 1 1366 do i = 1, nblk 1367 if (i .gt. 1) then 1368 k0 = k0 + 9*iblk(i-1)**2 1369 k1 = k1 + 3*iblk(i-1) 1370 end if 1371 k2 = k1 + 3*iblk(i) - 1 1372c 1373c scan over rows of the Hessian, first part 1374c 1375 l = 0 1376 do j = 1, 3*iblk(i) 1377 l2 = k1 + j - 1 1378 work(l2) = 0.0d0 1379 do k = k1, k2 1380 l = l + 1 1381 work(l2) = work(l2) + uu(k0+l)*rk(k) 1382 end do 1383 end do 1384c 1385c zero out the segment 1386c 1387 do k = k1, k2 1388 rk(k) = 0.0d0 1389 end do 1390c 1391c scan over rows of the Hessian, second part 1392c 1393 l = 0 1394 do j = 1, 3*iblk(i) 1395 l2 = k1 + j - 1 1396 do k = k1, k2 1397 l = l + 1 1398 rk(k) = rk(k) + uu(k0+l)*d(l2)*work(l2) 1399 end do 1400 end do 1401 end do 1402c 1403c perform deallocation of some local arrays 1404c 1405 deallocate (d) 1406 deallocate (work) 1407 return 1408 end 1409c 1410c 1411c ################################################################ 1412c ## ## 1413c ## subroutine gsort -- orthogonal vector via Gram-Schmidt ## 1414c ## ## 1415c ################################################################ 1416c 1417c 1418c "gsort" uses the Gram-Schmidt algorithm to build orthogonal 1419c vectors for sliding block interative matrix diagonalization 1420c 1421c 1422 subroutine gsort (nvar,nb,r0) 1423 use vibs 1424 implicit none 1425 integer i,j 1426 integer nvar,nb 1427 real*8 sum 1428 real*8 r0(*) 1429 real*8, allocatable :: s(:) 1430 real*8, allocatable :: proj(:) 1431c 1432c 1433c perform dynamic allocation of some local arrays 1434c 1435 allocate (s(nb)) 1436 allocate (proj(nvar)) 1437c 1438c make overlap between two basis sets 1439c 1440 do i = 1, nb 1441 s(i) = 0.0d0 1442 do j = 1, nvar 1443 s(i) = s(i) + r0(j)*rho(j,i) 1444 end do 1445 end do 1446c 1447c start the Gram-Schmidt procedure 1448c 1449 do i = 1, nvar 1450 proj(i) = 0.0d0 1451 end do 1452c 1453c construct projector 1454c 1455 do i = 1, nb 1456 do j = 1, nvar 1457 proj(j) = proj(j) + s(i)*rho(j,i) 1458 end do 1459 end do 1460c 1461c apply projector and normalize new vector 1462c 1463 sum = 0.0d0 1464 do i = 1, nvar 1465 proj(i) = r0(i) - proj(i) 1466 sum = sum + proj(i)*proj(i) 1467 end do 1468 sum = sqrt(sum) 1469 do i = 1, nvar 1470 proj(i) = proj(i) / sum 1471 end do 1472c 1473c return original array updated 1474c 1475 do i = 1, nvar 1476 r0(i) = proj(i) 1477 end do 1478c 1479c perform deallocation of some local arrays 1480c 1481 deallocate (s) 1482 deallocate (proj) 1483 return 1484 end 1485c 1486c 1487c ################################################################ 1488c ## ## 1489c ## subroutine qonvec -- block iterative vibration utility ## 1490c ## ## 1491c ################################################################ 1492c 1493c 1494c "qonvec" is a vector utility routine used during sliding 1495c block iterative matrix diagonalization 1496c 1497c 1498 subroutine qonvec (nvar,nr,u,rk,r) 1499 implicit none 1500 integer i,j,nvar,nr 1501 real*8 rku(6) 1502 real*8 rk(*) 1503 real*8 r(*) 1504 real*8 u(nvar,*) 1505c 1506c 1507c operate on vector rk with u-transpose 1508c 1509 do i = 1, nr 1510 rku(i) = 0.0d0 1511 do j = 1, nvar 1512 rku(i) = rku(i) + u(j,i)*rk(j) 1513 end do 1514 end do 1515c 1516c operate with u on the resultant 1517c 1518 do i = 1, nvar 1519 r(i) = 0.0d0 1520 do j = 1, nr 1521 r(i) = r(i) + u(i,j)*rku(j) 1522 end do 1523 end do 1524c 1525c subtract new product from r 1526c 1527 do i = 1, nvar 1528 r(i) = rk(i) - r(i) 1529 end do 1530 return 1531 end 1532c 1533c 1534c ################################################################# 1535c ## ## 1536c ## subroutine project -- remove known vectors from current ## 1537c ## ## 1538c ################################################################# 1539c 1540c 1541c "project" reads locked vectors from a binary file and projects 1542c them out of the components of the set of trial eigenvectors 1543c using the relation Y = X - U * U^T * X 1544c 1545c 1546 subroutine project (nvar,nconv,ivb1,ns,m) 1547 use vibs 1548 implicit none 1549 integer i,j,k 1550 integer nvar,nconv 1551 integer ivb1,ns,m 1552 real*8, allocatable :: temp(:) 1553 real*8, allocatable :: u(:) 1554c 1555c 1556c zero the temporary storage array 1557c 1558 do k = 1, nvar 1559 do i = 1, ns 1560 rwork(k,i+m) = 0.0d0 1561 end do 1562 end do 1563c 1564c perform dynamic allocation of some local arrays 1565c 1566 allocate (temp(ns)) 1567 allocate (u(nvar)) 1568c 1569c read and scan over the locked eigenvectors 1570c 1571 do i = 1, nconv 1572 read (ivb1) (u(k),k=1,nvar) 1573 do j = 1, ns 1574 temp(j) = 0.0d0 1575 do k = 1, nvar 1576 temp(j) = temp(j) + u(k)*rho(k,j+m) 1577 end do 1578 end do 1579 do j = 1, ns 1580 do k = 1, nvar 1581 rwork(k,j+m) = rwork(k,j+m) + u(k)*temp(j) 1582 end do 1583 end do 1584 end do 1585c 1586c perform deallocation of some local arrays 1587c 1588 deallocate (temp) 1589 deallocate (u) 1590c 1591c project locked vectors out of the current set 1592c 1593 do k = 1, nvar 1594 do i = 1, ns 1595 rho(k,i+m) = rho(k,i+m) - rwork(k,i+m) 1596 end do 1597 end do 1598 if (nconv .gt. 0) rewind (unit=ivb1) 1599 return 1600 end 1601c 1602c 1603c ################################################################## 1604c ## ## 1605c ## subroutine projectk -- remove known vectors from current ## 1606c ## ## 1607c ################################################################## 1608c 1609c 1610c "projectk" reads locked vectors from a binary file and projects 1611c them out of the components of the set of trial eigenvectors 1612c using the relation Y = X - U * U^T * X 1613c 1614c 1615 subroutine projectk (nvar,nconv,ivb1,ns,m) 1616 use vibs 1617 implicit none 1618 integer i,j,k 1619 integer nvar,nconv 1620 integer ivb1,ns,m 1621 real*8, allocatable :: temp(:) 1622 real*8, allocatable :: u(:) 1623c 1624c 1625c zero the temporary storage array 1626c 1627 do k = 1, nvar 1628 do i = 1, ns 1629 rwork(k,i+m) = 0.0d0 1630 end do 1631 end do 1632c 1633c perform dynamic allocation of some local arrays 1634c 1635 allocate (temp(ns)) 1636 allocate (u(nvar)) 1637c 1638c read and scan over the locked eigenvectors 1639c 1640 do i = 1, nconv 1641 read (ivb1) (u(k),k=1,nvar) 1642 do j = 1, ns 1643 temp(j) = 0.0d0 1644 do k = 1, nvar 1645 temp(j) = temp(j) + u(k)*rhok(k,j+m) 1646 end do 1647 end do 1648 do j = 1, ns 1649 do k = 1, nvar 1650 rwork(k,j+m) = rwork(k,j+m) + u(k)*temp(j) 1651 end do 1652 end do 1653 end do 1654c 1655c perform deallocation of some local arrays 1656c 1657 deallocate (temp) 1658 deallocate (u) 1659c 1660c project locked vectors out of the current set 1661c 1662 do k = 1, nvar 1663 do i = 1, ns 1664 rhok(k,i+m) = rhok(k,i+m) - rwork(k,i+m) 1665 end do 1666 end do 1667 if (nconv .gt. 0) rewind (unit=ivb1) 1668 return 1669 end 1670c 1671c 1672c ############################################################## 1673c ## ## 1674c ## subroutine konvec -- evaluate Hessian-vector product ## 1675c ## ## 1676c ############################################################## 1677c 1678c 1679c "konvec" finds a Hessian-vector product via finite-difference 1680c evaluation of the gradient based on atomic displacements 1681c 1682c 1683 subroutine konvec (nvar,xm,qe,uvec,kuvec) 1684 use atomid 1685 use atoms 1686 use units 1687 implicit none 1688 integer i,j,k,nvar 1689 real*8 e,term 1690 real*8 sum,eps 1691 real*8 xm(*) 1692 real*8 qe(*) 1693 real*8 uvec(*) 1694 real*8 kuvec(*) 1695 real*8, allocatable :: delta(:) 1696 real*8, allocatable :: grd1(:,:) 1697 real*8, allocatable :: grd2(:,:) 1698c 1699c 1700c estimate displacement based on total average 1701c 1702 sum = 0.0d0 1703 do i = 1, nvar 1704 sum = sum + uvec(i)*uvec(i)/xm(i) 1705 end do 1706c 1707c perform dynamic allocation of some local arrays 1708c 1709 allocate (delta(nvar)) 1710 allocate (grd1(3,n)) 1711 allocate (grd2(3,n)) 1712c 1713c store the coordinate displacements 1714c 1715 eps = 0.001d0 / sqrt(sum) 1716 do i = 1, nvar 1717 delta(i) = eps * uvec(i) / sqrt(xm(i)) 1718 end do 1719c 1720c compute the forward displacement 1721c 1722 do i = 1, n 1723 k = 3 * (i-1) 1724 x(i) = bohr * (qe(k+1)+delta(k+1)) 1725 y(i) = bohr * (qe(k+2)+delta(k+2)) 1726 z(i) = bohr * (qe(k+3)+delta(k+3)) 1727 end do 1728 call gradient (e,grd1) 1729c 1730c compute the backward displacement 1731c 1732 do i = 1, n 1733 k = 3 * (i-1) 1734 x(i) = bohr * (qe(k+1)-delta(k+1)) 1735 y(i) = bohr * (qe(k+2)-delta(k+2)) 1736 z(i) = bohr * (qe(k+3)-delta(k+3)) 1737 end do 1738 call gradient (e,grd2) 1739c 1740c update via finite differences 1741c 1742 term = 0.5d0 * bohr / (eps * hartree) 1743 do i = 1, n 1744 k = 3 * (i-1) 1745 do j = 1, 3 1746 kuvec(k+j) = term * (grd1(j,i)-grd2(j,i)) / sqrt(xm(k+j)) 1747 end do 1748 end do 1749c 1750c perform deallocation of some local arrays 1751c 1752 deallocate (delta) 1753 deallocate (grd1) 1754 deallocate (grd2) 1755 return 1756 end 1757c 1758c 1759c ################################################################# 1760c ## ## 1761c ## subroutine transform -- diagonalize trial basis vectors ## 1762c ## ## 1763c ################################################################# 1764c 1765c 1766c "transform" diagonalizes the current basis vectors to produce 1767c trial roots for sliding block iterative matrix diagonalization 1768c 1769c 1770 subroutine transform (ns,nb,h,c) 1771 implicit none 1772 integer i,j,k,ns,nb 1773 real*8 h(nb,*) 1774 real*8 c(nb,*) 1775 real*8, allocatable :: e1(:) 1776 real*8, allocatable :: h1(:) 1777 real*8, allocatable :: c1(:,:) 1778c 1779c 1780c perform dynamic allocation of some local arrays 1781c 1782 allocate (e1(ns)) 1783 allocate (h1((ns+1)*ns/2)) 1784 allocate (c1(ns,ns)) 1785c 1786c pack the upper triangle of matrix 1787c 1788 k = 0 1789 do i = 1, ns 1790 do j = i, ns 1791 k = k + 1 1792 h1(k) = h(i,j) 1793 end do 1794 end do 1795c 1796c perform the matrix diagonalization 1797c 1798 call diagq (ns,ns,h1,e1,c1) 1799c 1800c copy values into the return arrays 1801c 1802 do i = 1, ns 1803 do j = 1, ns 1804 h(i,j) = 0.0d0 1805 c(i,j) = c1(i,j) 1806 end do 1807 h(i,i) = e1(i) 1808 end do 1809c 1810c perform deallocation of some local arrays 1811c 1812 deallocate (e1) 1813 deallocate (h1) 1814 deallocate (c1) 1815 return 1816 end 1817c 1818c 1819c ############################################################# 1820c ## ## 1821c ## subroutine diagblk -- diagonalization for atom block ## 1822c ## ## 1823c ############################################################# 1824c 1825c 1826c "diagblk" performs diagonalization of the Hessian for a 1827c block of atoms within a larger system 1828c 1829c 1830 subroutine diagblk (k0,k1,n,vector,wres) 1831 implicit none 1832 integer i,j,k,m 1833 integer n,k0,k1 1834 real*8 wres(*) 1835 real*8 vector(*) 1836 real*8, allocatable :: hval(:) 1837 real*8, allocatable :: hres(:) 1838 real*8, allocatable :: hvec(:,:) 1839c 1840c 1841c perform dynamic allocation of some local arrays 1842c 1843 allocate (hval(n)) 1844 allocate (hres((n+1)*n/2)) 1845 allocate (hvec(n,n)) 1846c 1847c pack the upper triangle of matrix 1848c 1849 k = 0 1850 do i = 1, n 1851 m = k0 + (i-1)*n 1852 do j = i, n 1853 k = k + 1 1854 hres(k) = vector(m+j) 1855 end do 1856 end do 1857c 1858c perform the matrix diagonalization 1859c 1860 call diagq (n,n,hres,hval,hvec) 1861c 1862c copy values into return arrays 1863c 1864 k = 0 1865 do i = 1, n 1866 do j = 1, n 1867 k = k + 1 1868 vector(k0+k) = hvec(j,i) 1869 end do 1870 end do 1871 do i = 1, n 1872 wres(k1+i-1) = hval(i) 1873 end do 1874c 1875c perform deallocation of some local arrays 1876c 1877 deallocate (hval) 1878 deallocate (hres) 1879 deallocate (hvec) 1880 return 1881 end 1882c 1883c 1884c ######################################################### 1885c ## ## 1886c ## subroutine prtvib -- output of vibrational mode ## 1887c ## ## 1888c ######################################################### 1889c 1890c 1891c "prtvib" writes to an external disk file a series of 1892c coordinate sets representing motion along a vibrational 1893c normal mode 1894c 1895c 1896 subroutine prtvib (ivib,r) 1897 use atoms 1898 use files 1899 implicit none 1900 integer i,j,k 1901 integer ivib,ixyz 1902 integer lext,nview 1903 integer freeunit 1904 real*8 ratio 1905 real*8 r(*) 1906 real*8, allocatable :: xorig(:) 1907 real*8, allocatable :: yorig(:) 1908 real*8, allocatable :: zorig(:) 1909 character*7 ext 1910 character*240 xyzfile 1911c 1912c 1913c create a name for the vibrational displacement file 1914c 1915 lext = 3 1916 call numeral (ivib,ext,lext) 1917 xyzfile = filename(1:leng)//'.'//ext(1:lext) 1918 ixyz = freeunit () 1919 call version (xyzfile,'new') 1920 open (unit=ixyz,file=xyzfile,status='new') 1921c 1922c perform dynamic allocation of some local arrays 1923c 1924 allocate (xorig(n)) 1925 allocate (yorig(n)) 1926 allocate (zorig(n)) 1927c 1928c store the original atomic coordinates 1929c 1930 do i = 1, n 1931 xorig(i) = x(i) 1932 yorig(i) = y(i) 1933 zorig(i) = z(i) 1934 end do 1935c 1936c make file with plus and minus the current vibration 1937c 1938 nview = 3 1939 do i = -nview, nview 1940 ratio = dble(i) / dble(nview) 1941 do k = 1, n 1942 j = 3 * (k-1) 1943 x(k) = xorig(k) + ratio*r(j+1) 1944 y(k) = yorig(k) + ratio*r(j+2) 1945 z(k) = zorig(k) + ratio*r(j+3) 1946 end do 1947 call prtxyz (ixyz) 1948 end do 1949 close (unit=ixyz) 1950c 1951c restore the original atomic coordinates 1952c 1953 do i = 1, n 1954 x(i) = xorig(i) 1955 y(i) = yorig(i) 1956 z(i) = zorig(i) 1957 end do 1958c 1959c perform deallocation of some local arrays 1960c 1961 deallocate (xorig) 1962 deallocate (yorig) 1963 deallocate (zorig) 1964 return 1965 end 1966c 1967c 1968c ############################################################### 1969c ## ## 1970c ## subroutine hessblk -- Hessian elements for atom block ## 1971c ## ## 1972c ############################################################### 1973c 1974c 1975c "hessblk" calls subroutines to calculate the Hessian elements 1976c for each atom in turn with respect to Cartesian coordinates 1977c 1978c 1979 subroutine hessblk (amass,k0,i1,i2,vector) 1980 use atoms 1981 use bound 1982 use couple 1983 use hescut 1984 use hessn 1985 use inform 1986 use iounit 1987 use limits 1988 use mpole 1989 use potent 1990 use rigid 1991 use usage 1992 use vdw 1993 use vdwpot 1994 use units 1995 implicit none 1996 integer i,j,k 1997 integer ii,k0 1998 integer i1,i2 1999 real*8 ami,amik 2000 real*8 cutoff,rdn 2001 real*8 amass(*) 2002 real*8 vector(*) 2003 logical first 2004 save first 2005 data first / .true. / 2006c 2007c 2008c maintain any periodic boundary conditions 2009c 2010 if (use_bounds .and. .not.use_rigid) call bounds 2011c 2012c update the pairwise interaction neighbor lists 2013c 2014 if (use_list) call nblist 2015c 2016c many implicit solvation models require Born radii 2017c 2018 if (use_born) call born 2019c 2020c alter partial charges and multipoles for charge flux 2021c 2022 if (use_chgflx) call alterchg 2023c 2024c modify bond and torsion constants for pisystem 2025c 2026 if (use_orbit) call picalc 2027c 2028c compute the induced dipoles at polarizable atoms 2029c 2030 if (use_polar) then 2031 call chkpole 2032 call rotpole 2033 call induce 2034 end if 2035c 2036c calculate the "reduced" atomic coordinates 2037c 2038 if (use_vdw) then 2039 do i = 1, n 2040 ii = ired(i) 2041 rdn = kred(i) 2042 xred(i) = rdn*(x(i)-x(ii)) + x(ii) 2043 yred(i) = rdn*(y(i)-y(ii)) + y(ii) 2044 zred(i) = rdn*(z(i)-z(ii)) + z(ii) 2045 end do 2046 end if 2047c 2048c perform dynamic allocation of some global arrays 2049c 2050 if (first) then 2051 first = .false. 2052 if (.not. allocated(hessx)) allocate (hessx(3,n)) 2053 if (.not. allocated(hessy)) allocate (hessy(3,n)) 2054 if (.not. allocated(hessz)) allocate (hessz(3,n)) 2055 end if 2056c 2057c zero out the Hessian elements for the current atom 2058c 2059 ii = 0 2060 do i = i1, i2 2061 if (use(i)) then 2062 do k = i1, i2 2063 do j = 1, 3 2064 hessx(j,k) = 0.0d0 2065 hessy(j,k) = 0.0d0 2066 hessz(j,k) = 0.0d0 2067 end do 2068 end do 2069c 2070c remove any previous use of the replicates method 2071c 2072 cutoff = 0.0d0 2073 call replica (cutoff) 2074c 2075c call the local geometry Hessian component routines 2076c 2077 if (use_bond) call ebond2 (i) 2078 if (use_angle) call eangle2 (i) 2079 if (use_strbnd) call estrbnd2 (i) 2080 if (use_urey) call eurey2 (i) 2081 if (use_angang) call eangang2 (i) 2082 if (use_opbend) call eopbend2 (i) 2083 if (use_opdist) call eopdist2 (i) 2084 if (use_improp) call eimprop2 (i) 2085 if (use_imptor) call eimptor2 (i) 2086 if (use_tors) call etors2 (i) 2087 if (use_pitors) call epitors2 (i) 2088 if (use_strtor) call estrtor2 (i) 2089 if (use_angtor) call eangtor2 (i) 2090 if (use_tortor) call etortor2 (i) 2091c 2092c call the van der Waals Hessian component routines 2093c 2094 if (use_vdw) then 2095 if (vdwtyp .eq. 'LENNARD-JONES') call elj2 (i) 2096 if (vdwtyp .eq. 'BUCKINGHAM') call ebuck2 (i) 2097 if (vdwtyp .eq. 'MM3-HBOND') call emm3hb2 (i) 2098 if (vdwtyp .eq. 'BUFFERED-14-7') call ehal2 (i) 2099 if (vdwtyp .eq. 'GAUSSIAN') call egauss2 (i) 2100 end if 2101c 2102c call the electrostatic Hessian component routines 2103c 2104 if (use_charge) call echarge2 (i) 2105 if (use_chgdpl) call echgdpl2 (i) 2106 if (use_dipole) call edipole2 (i) 2107 if (use_mpole) call empole2 (i) 2108 if (use_polar) call epolar2 (i) 2109 if (use_rxnfld) call erxnfld2 (i) 2110c 2111c call any miscellaneous Hessian component routines 2112c 2113 if (use_solv) call esolv2 (i) 2114 if (use_metal) call emetal2 (i) 2115 if (use_geom) call egeom2 (i) 2116 if (use_extra) call extra2 (i) 2117c 2118c store Hessian for the current atom block as a vector 2119c 2120 ami = bohr**2 / (hartree*sqrt(amass(i))) 2121 do k = i1, i2 2122 amik = ami / sqrt(amass(k)) 2123 do j = 1, 3 2124 ii = ii + 1 2125 vector(k0+ii) = hessx(j,k) * amik 2126 end do 2127 end do 2128 do k = i1, i2 2129 amik = ami / sqrt(amass(k)) 2130 do j = 1, 3 2131 ii = ii + 1 2132 vector(k0+ii) = hessy(j,k) * amik 2133 end do 2134 end do 2135 do k = i1, i2 2136 amik = ami / sqrt(amass(k)) 2137 do j = 1, 3 2138 ii = ii + 1 2139 vector(k0+ii) = hessz(j,k) * amik 2140 end do 2141 end do 2142 end if 2143 end do 2144 return 2145 end 2146