1C 2c 11 jul 01: header fixed; send to Gijs 3c 7 jun 01: determine runtype, scftype, dft, mp2 from output 4c 6 jun 01: uhf fixed 5c 4 jun 01: prepares MOLDEN Format file from CADPAC output for: 6c { rhf | uhf } 7c { ene | opt | hess } 8c { dft | mp2 } 9 10c 11c Mariusz Klobukowski 12c cad2mol is a minor modification of Chris Lovallo's hon2mol 13c 14C This program will take a CADPAC output file and convert it to 15C Molden Format for use with Molden (as filename.mdn). The program 16C should work for SP energy, geometry optimizations, and Hessians for 17C any single-determinant wavefunction 18 19C For the orbitals to be read correctly, NOPRINT OCCVECTORS must not 20c be used 21 22 program cad2mol 23 implicit double precision(a-e,g-h,o-z) 24 implicit logical(f) 25 character filenm*70 26 parameter(inp=1,iout=2,itmp=3,maxbf=2000) 27 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 28 common /flags2/ flgmc 29 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 30 dimension occ(maxbf) 31 32C Open necessary files (input, output, and temp) 33 34 call open_files(inp,iout,itmp,filenm) 35 36C Print title onto files 37 38 call write_files(iout,itmp,'[Molden Format]') 39 40C Read and setup parameters of run 41 42 call find_flags(inp,occ) 43 rewind(inp) 44 45C Find and write coordinates and basis set 46 47 call find_coord(inp,iout,itmp) 48 49C Find and write eigenvectors and eigenvalues 50 51 call find_eigen(inp,iout,occ) 52 53C Next (and last for SP energy) is SCF energy convergence 54 55 rewind(inp) 56 call find_conv1(inp,iout) 57 58C End if the run is SP energy, skip to frequencies if run is Hessian 59 60 if (irun .eq. 0) goto 9999 61 if (irun .eq. 2) goto 10 62 63C Find and write out geometry optimization data 64 65 call find_gopt(inp,iout) 66 if (irun .eq. 3) goto 10 67 goto 9999 68 69C Find and write frequencies and normal modes (for Hessian runs) 70 71 10 call find_freq(inp,iout,itmp) 72 73C Close files 74 75 9999 call close_files(inp,iout,itmp) 76 77C That's it! It's over! Bye now! 78 79 end 80C----------------------------------------------------------------------- 81 subroutine open_files(inp,iout,itmp,filenm) 82 character filenm*70 83 84C Open input file (CADPAC output file) 85 86 write(6,*) 'Opening necessary files...' 87 call getarg(1,filenm) 88 if(filenm.eq.' ') then 89 write(6,fmt='('' Enter the input file name (with .out): '',$)') 90 read (5,fmt='(a70)') filenm 91 endif 92 open(unit=inp,file=filenm,access='SEQUENTIAL', 93 1 status='OLD',form='FORMATTED') 94 95C Open output and temp files 96 97 do i=1,67 98 if (filenm(i:i+3) .eq. '.out') filenm(i:i+3)='.mdn' 99 enddo 100 open(unit=iout,file=filenm,access='SEQUENTIAL', 101 1 status='UNKNOWN',form='FORMATTED') 102 open(unit=itmp,access='SEQUENTIAL', 103 1 status='UNKNOWN',name='TEMP',form='FORMATTED') 104c 1 status='SCRATCH',form='FORMATTED') 105 return 106 end 107C----------------------------------------------------------------------- 108 subroutine find_flags(inp,occ) 109 implicit double precision(a-e,g-h,o-z) 110 implicit logical(f) 111 character line*130 112c character kyword*26 113 parameter(maxbf=2000) 114 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 115 common /flags2/ flgmc 116 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 117 dimension occ(*) 118 119C This subroutine searches the CADPAC file for various data pertaining 120C to the run (type of wavefunction, number of basis functions and 121C occupied orbitals, etc.) 122 123 write(6,*) 'Beginning initialization...' 124 125C Initialize flags 126 127 flgrhf=.false. 128 flguhf=.false. 129 flgopn=.false. 130 flgdft=.false. 131 flggvb=.false. 132 flgmp2=.false. 133 flgmp4=.false. 134 flgmc= .false. 135 136C Search for type of run and properties - decode from command line 137C Determine type of wavefunction - decode from command line 138c cad2mol job.out 139c { rhf | uhf | rohf | gvb | mp2 | mp4 } 140c { ene | opt | hess } 141c { dft } 142 143* inparg=iargc() 144* write(6,'('' inparg:'',i3)') inparg 145 146 irun = -1 147 iwfn = -1 148* do narg=2,inparg 149* call getarg(narg,kyword) 150* write(6,'('' narg:'',i3,'' kyword:'',a)') narg,kyword 151* call up_case(kyword,26) 152* call last_char(kyword,26,loc) 153* if (kyword(1:loc).eq.'RHF') then 154* flgrhf=.true. 155* iwfn=0 156* else if(kyword(1:loc).eq.'ROHF') then 157* flguhf=.true. 158* flgopn=.true. 159* iwfn=0 160* else if(kyword(1:loc).eq.'GVB') then 161* flguhf=.true. 162* flgopn=.true. 163* iwfn=0 164* else if(kyword(1:loc).eq.'UHF') then 165* flguhf=.true. 166* flgopn=.true. 167* iwfn=0 168* else if(kyword(1:loc).eq.'MP2') then 169* iwfn=5 170* flgmp2=.true. 171* else if(kyword(1:loc).eq.'MP4') then 172* iwfn=6 173* flgmp4=.true. 174* else if(kyword(1:loc).eq.'CAS') then 175* iwfn=1 176* else if(kyword(1:loc).eq.'DFT') then 177* flgdft=.true. 178* else if(kyword(1:loc).eq.'ENE') then 179* irun=0 180* else if(kyword(1:loc).eq.'OPT') then 181* irun=1 182* else if(kyword(1:loc).eq.'HESS') then 183* irun=2 184* else if(kyword(1:loc).eq.'OPHS') then 185* irun=3 186* else 187* write(6,*) 'wrong parameter on command line' 188* call error 189* endif 190* enddo 191 192 rewind(inp) 193 ifound=0 194 call search(inp,line,'RUN TYPE = ENERGY',ifound,0) 195 if(ifound.eq.1) then 196 write(6,*) 'RUN TYPE = ENERGY' 197 irun=0 198 endif 199 rewind(inp) 200 ifound=0 201 call search(inp,line,'RUN TYPE = OPTIMIZE',ifound,0) 202 if(ifound.eq.1) then 203 write(6,*) 'RUN TYPE = OPTIMIZE' 204 irun=1 205 endif 206 rewind(inp) 207 ifound=0 208 call search(inp,line,'RUN TYPE = SECDER',ifound,0) 209 if(ifound.eq.1.and.irun.eq.-1) then 210 irun=2 211 write(6,*) 'RUN TYPE = SECDER -- only' 212 endif 213 if(ifound.eq.1.and.irun.eq.1) then 214 irun=3 215 write(6,*) 'RUN TYPE = SECDER -- after OPTIMIZE' 216 endif 217 rewind(inp) 218 ifound=0 219 call search(inp,line,'RUN TYPE = FORCE',ifound,0) 220 if(ifound.eq.1.and.irun.eq.-1) then 221 irun=2 222 write(6,*) 'RUN TYPE = FORCE -- only' 223 endif 224 if(ifound.eq.1.and.irun.eq.1) then 225 irun=3 226 write(6,*) 'RUN TYPE = FORCE -- after OPTIMIZE' 227 endif 228 rewind(inp) 229 ifound=0 230 call search(inp,line,'SCF TYPE = CLOSED',ifound,0) 231 if(ifound.eq.1) then 232 flgrhf=.true. 233 iwfn=0 234 write(6,*) 'SCF TYPE = CLOSED' 235 endif 236 rewind(inp) 237 ifound=0 238 call search(inp,line,'UHF Calculation',ifound,0) 239 if(ifound.eq.1) then 240 iwfn=0 241 flguhf=.true. 242 flgopn=.true. 243 write(6,*) 'UHF Calculation' 244 endif 245 rewind(inp) 246 ifound=0 247 call search(inp,line,'CLOSED-SHELL KOHN-SHAM',ifound,0) 248 if(ifound.eq.1) then 249 flgdft=.true. 250 write(6,*) 'CLOSED-SHELL KOHN-SHAM' 251 endif 252 rewind(inp) 253 ifound=0 254 call search(inp,line,'Unrestricted KOHN-SHAM',ifound,0) 255 if(ifound.eq.1) then 256 iwfn=0 257 flgdft=.true. 258 flguhf=.true. 259 flgopn=.true. 260 write(6,*) 'Unrestricted KOHN-SHAM' 261 endif 262 rewind(inp) 263 ifound=0 264 call search(inp,line,'RHF MP2 Calculation',ifound,0) 265 if(ifound.eq.1) then 266 iwfn=5 267 flgmp2=.true. 268 write(6,*) 'RHF MP2 Calculation' 269 endif 270 rewind(inp) 271 ifound=0 272 call search(inp,line,'UHF MP2 Calculation',ifound,0) 273 if(ifound.eq.1) then 274 iwfn=5 275 flgmp2=.true. 276 flguhf=.true. 277 flgopn=.true. 278 write(6,*) 'UHF MP2 Calculation' 279 endif 280 281 write(6,fmt='('' irun ='',i2)') irun 282 write(6,fmt='('' iwfn ='',i2)') iwfn 283 write(6,fmt='('' flgrhf ='',L2)') flgrhf 284 write(6,fmt='('' flguhf ='',L2)') flguhf 285 write(6,fmt='('' flgopn ='',L2)') flgopn 286 write(6,fmt='('' flgdft ='',L2)') flgdft 287 write(6,fmt='('' flgmp2 ='',L2)') flgmp2 288 289 if (irun.ne.0. and. 290 1 irun.ne.1 .and. 291 1 irun.ne.2 .and. 292 1 irun.ne.3) then 293 write(6,*) 'Only RUNTYPs=0-3 (SP energy, geom opt, and Hessian) 294 1 are currently implemented!' 295 call error 296 endif 297 298 if (iwfn .ne. 0 .and. iwfn .ne. 5 .and. iwfn .ne. 1) then 299 write(6,*) 'Only single-determinant and MCSCF wavefunctions are 300 1 currently implemented!' 301 call error 302 endif 303 304 rewind (inp) 305 306C Find number of basis functions 307 308 call search(inp,line,'Number of Basis Functions',ifound,-1) 309 read(unit=line,fmt='(46x,i5)') nbf 310 write(*,*) 'nbf =',nbf 311 if (nbf .ge. maxbf) then 312 write(6,*) 'Error code 5! Number of basis functions is greater 313 1 than parameter MAXBF!' 314 write(6,*) 'Check parameter MAXBF and run program again.' 315 call error 316 endif 317 318 319C Determine type of run and orbital occupancies 320 321 if (iwfn .eq. 1 .or. iwfn .eq. 10 .or. iwfn .eq. 12) flgmc=.true. 322 if (flgmc .or. iwfn .eq. 2) then 323 call search(inp,line,'NUMBER OF ORBITALS',ifound,-1) 324 read(unit=line,fmt='(27x,i5)') norb 325 call search(inp,line,'# OF CORE ORBITALS',ifound,-1) 326 read(unit=line,fmt='(26x,i5)') ndoc 327 goto 20 328 endif 329 330C The program does not search for MP2/4 lines, it sets the 331C flags according to iwfn 332 333 10 if (iwfn .eq. 5) flgmp2=.true. 334 if (iwfn .eq. 6) flgmp4=.true. 335 if(flguhf) then 336 call search(inp,line,'occupied alpha orbitals',ifound,-1) 337 read(unit=line,fmt='(46x,i5)') ndoc 338 else 339 call search(inp,line, 340 1 'Number of doubly occupied orbitals',ifound,-1) 341 read(unit=line,fmt='(46x,i5)') ndoc 342 endif 343 write(*,*) 'ndoc =',ndoc 344 if (.not. flgopn) goto 20 345 if (flguhf) then 346 call search(inp,line,'occupied beta orbitals',ifound,-1) 347 read(unit=line,fmt='(46x,i5)') nsoc 348 write(*,*) 'nsoc=',nsoc 349 goto 20 350 endif 351 352 20 if (.not. flguhf) then 353 do i=1,ndoc 354 occ(i)=2.0 355 enddo 356 if (nsoc .gt. 0) then 357 do i=ndoc+1,ndoc+1+nsoc 358 occ(i)=1.0 359 enddo 360 do i=ndoc+1+nsoc,nbf 361 occ(i)=0.0 362 enddo 363 else 364 do i=ndoc+1,nbf 365 occ(i)=0.0 366 enddo 367 endif 368 goto 30 369 else 370 do i=1,ndoc 371 occ(i)=1.0 372 enddo 373 do i=ndoc+1,nbf 374 occ(i)=0.0 375 enddo 376 endif 377 378C Determine maximum number of iterations (SCF or CI) 379 380 30 if (iwfn .eq. 2) then 381 call search(inp,line,'MAX. NUMB. OF ITERATIONS',ifound,-1) 382 read(unit=line,fmt='(29x,i6)') maxit 383 else 384 if(flguhf) then 385 call search(inp,line,'Maximum Number of Iterations',ifound,-1) 386 else 387 call search(inp,line,'Maximum number of iterations',ifound,-1) 388 endif 389 read(unit=line,fmt='(32x,i6)') maxit 390 write(*,*) 'maxit=',maxit 391 endif 392 maxit=maxit+1 393 394 write(6,*) 'Initialization done!' 395 return 396 end 397C----------------------------------------------------------------------- 398 subroutine find_coord(inp,iout,itmp) 399 implicit double precision(a-e,g-h,o-z) 400 implicit logical(f) 401 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 402 character line*130,lintmp*130,atmnam*8 403 parameter(maxat=100) 404 dimension atmnam(maxat) 405 406C This subroutine finds and writes out the coordinates of the atoms 407C It searches for CHARGE, and then skips the --- and blank lines 408C It then calls find_bas to write out the basis set 409 410 write(6,*) 'Reading coordinates...' 411 call write_files(iout,itmp,'[Atoms] AU') 412 call search(inp,line,'Atom Nuclear',ifound,-1) 413 call skip(inp,2) 414 i=1 415 10 call read_line(inp,line) 416 read(unit=line,fmt='(4x,a8)') atmnam(i) 417 if (atmnam(i) .eq. '--------') goto 20 418 read(unit=line,fmt='(4x,a8,2x,f5.1,3f15.8)') atmnam(i),chrg 419 1 ,x,y,z 420 ichrg=int(chrg) 421 write(*,fmt='(1x,a8,i3,i3,3(2x,f12.8))') atmnam(i),i,ichrg,x,y,z 422 write(unit=lintmp,fmt='(1x,a8,i3,i3,3(2x,f12.7))') atmnam(i),i 423 1 ,ichrg,x,y,z 424 call write_files(iout,itmp,lintmp) 425 i=i+1 426 goto 10 427 20 nat=i-1 428 write(*,*) 'nat=',nat 429 write(6,*) 'Done writing coordinates!' 430 call find_bas(inp,iout,itmp,atmnam) 431 return 432 end 433C----------------------------------------------------------------------- 434 subroutine find_bas(inp,iout,itmp,atmnam) 435 implicit double precision(a-e,g-h,o-z) 436 implicit logical(f) 437 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 438 character line*130,lintmp*130,temp8*8,temp4*4,atmnam*8,shltyp*4 439 parameter(mxprim=30) 440 dimension atmnam(*),coef1(mxprim),coef2(mxprim),exp(mxprim) 441 442C This subroutine reads and writes out the basis set, including looping 443C over the non-symmetry unique atoms 444 445 write(6,*) 'Reading basis set...' 446 call skip(inp,6) 447 call write_files(iout,itmp,'[GTO]') 448 do j=1,nat 449* write(*,*) j,atmnam(j) 450c ... reset cgtf count for each new atom 451 ikount=0 452 write(unit=lintmp,fmt='(i3,a5)') j, ' 0 ' 453 call write_files(iout,itmp,lintmp) 454 call read_line(inp,line) 455* write(*,*) '>'//line(1:76) 456 read(unit=line,fmt='(1x,a8)') temp8 457* write(*,*) 'temp8'//temp8 458 if (temp8 .ne. atmnam(j)) then 459C Either an error has occurred, or this is not a symmetry unique atom. 460C Checking for an atom with the same name and copying basis set... 461 backspace(inp) 462 temp8=atmnam(j) 463 do k=1,j-1 464 if (temp8 .eq. atmnam(k)) then 465 rewind(itmp) 466 call search(itmp,line,atmnam(k),ifound,-1) 467 read(unit=line,fmt='(10x,i3)') n 468 write(unit=temp8,fmt='(i3,a5)') n, ' 0 ' 469 call search(itmp,line,temp8,ifound,-1) 470 write(unit=temp8,fmt='(i3,a5)') n+1, ' 0 ' 471 10 call read_line(itmp,line) 472* write(*,*) '>'//line(1:76) 473 if (line(1:8) .ne. temp8) then 474 call write_files(iout,0,line) 475 goto 10 476 endif 477 call ffwrd(itmp) 478 goto 40 479 endif 480 enddo 481 482C Nope! There was an error! Stopping... 483 484 write(6,*) 'Error code 1 reading basis set!' 485 call error 486 endif 487 488 l=1 489 20 call read_line(inp,line) 490* write(*,*) '>'//line(1:76) 491 if(line(1:6).eq.' Model') call read_line(inp,line) 492 read(unit=line,fmt='(12x,i2,5x,a4,6x,f16.6,f15.6)') 493 1 jkount,temp4,exp(l),coef1(l) 494* write(*,*) 'l ikount jkount',l,ikount,jkount,exp(l) 495 if(jkount.eq.0) then 496* backspace(inp) 497 ikount=jkount 498 goto 30 499 endif 500 if(ikount.gt.0 .and. jkount.ne.ikount) then 501 ikount=jkount 502 backspace(inp) 503 goto 30 504 else 505 ikount=jkount 506 endif 507 shltyp=temp4 508 if (shltyp(1:1) .eq. 'L') then 509 read(unit=line,fmt='(19x,a4,6x,f16.6,2f15.6)') 510 1 shltyp,exp(l),coef1(l),coef2(l) 511 endif 512 l=l+1 513 goto 20 514 30 nprim=l-1 515* write(*,*) 'nprim exp',nprim,exp(nprim) 516 517C Write basis set 518 519 if (shltyp(1:1) .ne. 'L') then 520 write(unit=lintmp,fmt='(a2,i5,a5)') shltyp(1:1),nprim, ' 1.00' 521 call write_files(iout,itmp,lintmp) 522 do l=1,nprim 523 write(unit=lintmp,fmt='(2(1x,d17.10))') exp(l),coef1(l) 524 call write_files(iout,itmp,lintmp) 525 enddo 526 elseif (shltyp(1:1) .eq. 'L') then 527 write(unit=lintmp,fmt='(1x,a2,i4,a5)') 'SP',nprim, 528 1 ' 1.00' 529 call write_files(iout,itmp,lintmp) 530 do l=1,nprim 531 write(unit=lintmp,fmt='(3(1x,d17.10))') exp(l),coef1(l) 532 1 ,coef2(l) 533 call write_files(iout,itmp,lintmp) 534 enddo 535 else 536 write(6,*) 'Error code 3 writing basis set!' 537 call error 538 endif 539 if(jkount.eq.0) goto 40 540 l=1 541 goto 20 542 40 call blank_line(iout) 543 enddo 544 call blank_line(iout) 545 write(6,*) 'Done writing basis set!' 546 return 547 end 548C----------------------------------------------------------------------- 549 subroutine find_eigen(inp,iout,occ) 550 implicit double precision(a-e,g-h,o-z) 551 implicit logical(f) 552 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 553 common /flags2/ flgmc 554 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 555 character line*130,lintmp*130 556 parameter(maxbf=2000) 557 parameter (zero=0.0d+00, one=1.0d+00) 558 dimension occ(*),orbs(maxbf,maxbf),eval(maxbf) 559 560C This subroutine finds and writes out the molecular orbitals and 561C orbital energies, along with their occupancies 562 563 write(6,*) 'Reading orbitals and eigenvalues...' 564 call write_files(iout,0,'[MO]') 565c ... 1st MO output 566 first_SCF=.true. 567 if(flguhf) then 568 call search(inp,line,'Alpha-spin Orbitals',ifound,-1) 569 else 570 call search(inp,line,'Molecular Orbitals',ifound,-1) 571 endif 572c ... 2nd MO output - used if opt+hess in the same run 573 if (irun .eq. 1 .or. irun.eq.2 .or. irun.eq.3) then 574 if(flguhf) then 575 call search(inp,line,'Alpha-spin orbitals',ifound,0) 576 else 577 call search(inp,line,'Molecular orbitals',ifound,0) 578 endif 579 endif 580 first_SCF=.false. 581 if(ifound.eq.0) then 582 write(*,*) '... only hessian' 583 first_SCF=.true. 584 rewind(inp) 585 if(flguhf) then 586 call search(inp,line,'Alpha-spin Orbitals',ifound,-1) 587 else 588 call search(inp,line,'Molecular Orbitals',ifound,-1) 589 endif 590 endif 591c ... locate the 1st eigenvalue 592 call search(inp,line,' 1 ',ifound,-1) 593 backspace(inp) 594c ... read the number of MOs printed 595 found=.false. 596 nmos=0 597 do while (.not.found) 598 call read_line(inp,line) 599 if(line(1:20).ne.' ' .and. 600 1 line(1:13).ne.' Eigenvectors') then 601 nmos=nmos+1 602 else 603 found=.true. 604 backspace(inp) 605 endif 606 enddo 607 write(*,*) 'nmos=',nmos 608 write(*,*) 'first_SCF ',first_SCF 609 if(.not.flguhf) then 610 if(irun.eq.0.or.first_SCF) call skip(inp,2) 611 endif 612 if(flguhf.and.first_SCF) call skip(inp,1) 613 if (flggvb) call skip(inp,1) 614 j=1 615 10 k=j+4 616 call skip(inp,2) 617 call read_line(inp,line) 618 read(unit=line,fmt='(12x,5(f13.6))') (eval(i), i=j,k) 619 write(*,fmt='(5(i4,f13.6))') (i,eval(i), i=j,k) 620 call skip(inp,1) 621 do l=1,nbf 622 call read_line(inp,line) 623 read(unit=line,fmt='(12x,5(f13.6))') (orbs(i,l), i=j,k) 624 enddo 625 if (k .ge. nmos) then 626 do i=1,nmos 627 write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i) 628 call write_files(iout,0,lintmp) 629 call write_files(iout,0,' Spin= Alpha') 630 if(flguhf.and.i.gt.ndoc) occ(i)=zero 631 write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i) 632 call write_files(iout,0,lintmp) 633 do j=1,nbf 634 write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j) 635 call write_files(iout,0,lintmp) 636 enddo 637 enddo 638 if (.not. flguhf) goto 20 639 640C Now read beta orbitals (for UHF runs only) 641 642 write(6,*) 'Done writing alpha orbitals...' 643 call find_beta(inp,iout,occ) 644 write(6,*) 'Done writing beta orbitals...' 645 goto 20 646 endif 647 j=j+5 648 goto 10 649 20 write(6,*) 'Done writing orbitals and eigenvalues!' 650 return 651 end 652C----------------------------------------------------------------------- 653 subroutine find_beta(inp,iout,occ) 654 655C This subroutine finds and write the beta molecular orbitals and 656C orbital energies for a UHF wavefunction 657 658 implicit double precision(a-e,g-h,o-z) 659 implicit logical(f) 660 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 661 common /flags2/ flgmc 662 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 663 character line*130,lintmp*130 664 parameter(maxbf=2000) 665 dimension occ(*),orbs(maxbf,maxbf),eval(maxbf) 666 667c ... locate start of current beta's 668 call search(inp,line,'Beta-spin',ifound,-1) 669c ... locate 1st beta MO 670 call search(inp,line,' 1 ',ifound,-1) 671 backspace(inp) 672c ... read the number of MOs printed 673 found=.false. 674 nmos=0 675 do while (.not.found) 676 call read_line(inp,line) 677 if(line(1:20).ne.' '.and. 678 1 line(1:13).ne.' Eigenvectors') then 679 nmos=nmos+1 680 else 681 found=.true. 682 backspace(inp) 683 endif 684 enddo 685 write(*,*) 'nmos=',nmos 686 call read_line(inp,line) 687 if(line(1:13).ne.' Eigenvectors') backspace(inp) 688 j=1 689 10 k=j+4 690 call skip(inp,2) 691 call read_line(inp,line) 692 read(unit=line,fmt='(12x,5(f13.6))') (eval(i), i=j,k) 693 write(*,fmt='(5(i4,f13.6))') (i,eval(i), i=j,k) 694 call skip(inp,1) 695 do l=1,nbf 696 call read_line(inp,line) 697 read(unit=line,fmt='(12x,5(f13.6))') (orbs(i,l), i=j,k) 698 enddo 699 if (k .ge. nmos) then 700 do i=1,nmos 701 write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i) 702 call write_files(iout,0,lintmp) 703 call write_files(iout,0,' Spin= Beta') 704 if(flguhf.and.i.gt.nsoc) occ(i)=zero 705 write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i) 706 call write_files(iout,0,lintmp) 707 do j=1,nbf 708 write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j) 709 call write_files(iout,0,lintmp) 710 enddo 711 enddo 712 goto 20 713 endif 714 j=j+5 715 goto 10 716 20 return 717 end 718C---------------------------------------------------------------------- 719 subroutine find_eigmc(inp,iout,occ) 720 implicit double precision(a-e,g-h,o-z) 721 implicit logical(f) 722 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 723 common /flags2/ flgmc 724 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 725 character line*130,lintmp*130 726 parameter(maxbf=2000) 727 dimension occ(*),orbs(maxbf,maxbf),eval(maxbf) 728 729C This subroutine finds and writes out the natural orbitals and their 730C occupancies (energies are all set to 0 if the Fock energies are not 731C printed) 732 733 do i=1,nbf 734 eval(i)=0.0 735 enddo 736 rewind(inp) 737 call search(inp,line,'EIGENVALUES OF -MC-',ifound,0) 738 if (ifound .eq. 1) then 739 call skip(inp,2) 740 j=1 741 10 k=j+5 742 call read_line(inp,line) 743 read(unit=line,fmt='(6(10x,f11.4))') (eval(i), i=j,k) 744 j=j+6 745 if (k .lt. nbf) goto 10 746 endif 747 rewind(inp) 748 call write_files(iout,0,'[MO]') 749 if (irun .ne. 1) goto 30 750 n=0 751 20 call search(inp,line,'CAN.CORE + ) NAT.VAL',ifound,0) 752 if (ifound .eq. 1) then 753 n=n+1 754 goto 20 755 endif 756 rewind(inp) 757 do i=1,n 758 call search(inp,line,'CAN.CORE + ) NAT.VAL',ifound,-1) 759 enddo 760 goto 40 761 30 call search(inp,line,'CAN.CORE + ) NAT.VAL',ifound,-1) 762 40 call skip(inp,9) 763 call read_line(inp,line) 764 j=1 765 50 k=j+6 766 read(unit=line,fmt='(15x,7(f15.10))') (occ(i), i=j,k) 767 call skip(inp,2) 768 do l=1,nbf 769 call read_line(inp,line) 770 read(unit=line,fmt='(15x,7(f15.10))') (orbs(i,l), i=j,k) 771 enddo 772 if (k .ge. norb) then 773 do i=1,norb 774 write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i) 775 call write_files(iout,0,lintmp) 776 call write_files(iout,0,' Spin= Alpha') 777 write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i) 778 call write_files(iout,0,lintmp) 779 do j=1,nbf 780 write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j) 781 call write_files(iout,0,lintmp) 782 enddo 783 enddo 784 goto 60 785 endif 786 j=j+7 787 call skip(inp,8) 788 call read_line(inp,line) 789 goto 50 790 60 return 791 end 792C----------------------------------------------------------------------- 793 subroutine find_eigci(inp,iout,occ) 794 implicit double precision(a-e,g-h,o-z) 795 implicit logical(f) 796 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 797 common /flags2/ flgmc 798 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 799 character line*130,lintmp*130 800 parameter(maxbf=2000) 801 dimension occ(*),orbs(maxbf,maxbf),eval(maxbf) 802 803C This subroutine finds and writes out the natural orbitals and their 804C occupancies (energies are all set to 0) 805 806 do i=1,nbf 807 eval(i)=0.0 808 enddo 809 rewind(inp) 810 call write_files(iout,0,'[MO]') 811 if (irun .ne. 1) goto 30 812 n=0 813 20 call search(inp,line,'NATURAL ORBITALS IN ATOMIC',ifound,0) 814 if (ifound .eq. 1) then 815 n=n+1 816 goto 20 817 endif 818 rewind(inp) 819 do i=1,n 820 call search(inp,line,'NATURAL ORBITALS IN ATOMIC',ifound,-1) 821 enddo 822 goto 40 823 30 call search(inp,line,'NATURAL ORBITALS IN ATOMIC',ifound,-1) 824 40 call skip(inp,6) 825 call read_line(inp,line) 826 j=1 827 50 k=j+6 828 read(unit=line,fmt='(15x,7(f15.10))') (occ(i), i=j,k) 829 call skip(inp,2) 830 do l=1,nbf 831 call read_line(inp,line) 832 read(unit=line,fmt='(15x,7(f15.10))') (orbs(i,l), i=j,k) 833 enddo 834 if (k .ge. norb) then 835 do i=1,norb 836 write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i) 837 call write_files(iout,0,lintmp) 838 call write_files(iout,0,' Spin= Alpha') 839 write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i) 840 call write_files(iout,0,lintmp) 841 do j=1,nbf 842 write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j) 843 call write_files(iout,0,lintmp) 844 enddo 845 enddo 846 goto 60 847 endif 848 j=j+7 849 call skip(inp,5) 850 call read_line(inp,line) 851 goto 50 852 60 return 853 end 854C----------------------------------------------------------------------- 855 subroutine find_conv1(inp,iout) 856 implicit double precision(a-e,g-h,o-z) 857 implicit logical(f) 858 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 859 common /flags2/ flgmc 860 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 861 character line*130,lintmp*130,temp4*4 862 parameter(maxits=101) 863 dimension etot(maxits) 864 865C This subroutine finds and writes out the SCF energy convergence 866C (the first one) 867 868 write(6,*) 'Reading energy convergence...' 869 if(flgdft) then 870 write(6,*) '... DFT calculations - no energy convergence' 871 return 872 endif 873 if (maxit .ge. maxits) then 874 write(6,*) 'Error code 5! Number of iterations is greater than 875 1 parameter MAXITS!' 876 write(6,*) 'Check parameter MAXITS and run program again.' 877 call error 878 endif 879 if (flgmc) then 880 call search(inp,line,' CYCLE TOTAL ENERGY',ifound,-1) 881 else if (iwfn .eq. 2) then 882 call search(inp,line,' ITER. MAX.DEV.',ifound,-1) 883 call skip(inp,1) 884 else 885* call search(inp,line,'Cycle Total energy',ifound,-1) 886 call search(inp,line,'Cycle ',ifound,-1) 887 if(.not.flguhf) call skip(inp,1) 888 endif 889 i=0 890 finish=.false. 891 do while(.not.finish) 892 call read_line(inp,line) 893 if (iwfn .eq. 2) then 894 read(unit=line,fmt='(a4,23x,f16.9)') temp4,etot(i) 895 else 896 if(line(1:5).ne.' '.and.line(1:4).ne.' Den' 897 1 .and.line(1:5).ne.' Last' 898 1 .and.line(1:5).ne.' HOMO') then 899 i=i+1 900 read(unit=line,fmt='(5x,f18.9)') etot(i) 901* write(*,*) 'i e',i,etot(i) 902 endif 903 endif 904 if (line(1:4) .eq. ' Den') then 905 finish=.true. 906 endif 907 enddo 908 ncyc=i 909 10 call write_files(iout,0,'[SCFCONV]') 910 write(unit=lintmp,fmt='(a22,i5)') 'scf-first 1 THROUGH',ncyc 911 call write_files(iout,0,lintmp) 912 do i=1,ncyc 913 write(unit=lintmp,fmt='(f14.6)') etot(i) 914 call write_files(iout,0,lintmp) 915 enddo 916 write(6,*) 'Finished writing energy convergence!' 917 return 918 end 919C----------------------------------------------------------------------- 920 subroutine find_gopt(inp,iout) 921 922C This subroutine finds and writes the data pertaining to the geometry 923C optimization process. 924 925 implicit double precision(a-e,g-h,o-z) 926 implicit logical(f) 927 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 928 common /flags2/ flgmc 929 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 930 character line*130,lintmp*130,temp4*4,temp8*8,atmnam*8 931 parameter(maxstp=200,bohr2a=5.29177249D-01,maxits=101,maxat=100) 932 dimension estep(maxstp),xmxfor(maxstp),rmsfor(maxstp) 933 dimension ecyc(maxits),atmnam(maxat),geom(3,maxat,maxstp) 934* dimension grdint(5*maxat) 935 936 write(6,*) 'Reading geometry optimization information...' 937 938c ... DFT calculations do not show SCF convergence 939c ... irrelevant for MP2 940 if(flgdft .or. flgmp2) goto 25 941 rewind(inp) 942 nstep=0 943 10 if (flgmc) then 944 call search(inp,line,' CYCLE TOTAL ENERGY',ifound,0) 945 else if (iwfn .eq. 0) then 946 if(flguhf) then 947 call search(inp,line,' UHF Calculation',ifound,0) 948 else 949 call search(inp,line,' CLOSED-SHELL SCF CALCULATION',ifound,0) 950 endif 951 else 952 call search(inp,line,' Cycle ',ifound,0) 953 endif 954 if (ifound .eq. 1) then 955 nstep=nstep+1 956 goto 10 957 endif 958 write(*,*) 'nstep=',nstep 959 960c ... position to the last SCF 961 rewind(inp) 962 do i=1,nstep 963 if (flgmc) then 964 call search(inp,line,' CYCLE TOTAL ENERGY',ifound,-1) 965 else if (iwfn .eq. 0) then 966 if(flguhf) then 967 call search(inp,line,' UHF Calculation',ifound,0) 968 else 969 call search(inp,line,' CLOSED-SHELL SCF CALCULATION',ifound,0) 970 endif 971 call skip(inp,1) 972 else 973 call search(inp,line,' CYCLE TOTAL ENERGY',ifound,-1) 974 call skip(inp,1) 975 endif 976 enddo 977c ... located last SCF; find 1st scf cycle 978 call search(inp,line,' 1 ',ifound,-1) 979 backspace(inp) 980 i=0 981 finish=.false. 982 do while(.not.finish) 983 call read_line(inp,line) 984 if (iwfn .eq. 2) then 985 read(unit=line,fmt='(a4,23x,f16.9)') temp4,ecyc(i) 986 else 987 if(line(1:5).ne.' '.and.line(1:4).ne.' Den' 988 1 .and.line(1:5).ne.' Last') then 989 i=i+1 990 read(unit=line,fmt='(5x,f18.9)') ecyc(i) 991 write(*,*) 'i e',i,ecyc(i) 992 endif 993 endif 994 if (line(1:4) .eq. ' Den') then 995 finish=.true. 996 endif 997 enddo 998 ncyc=i 999 write(*,*) 'ncyc=',ncyc 1000c ... write scf convergence - last SCF 1001 20 write(unit=lintmp,fmt='(a21,i5)') 'scf-last 1 THROUGH',ncyc 1002 call write_files(iout,0,lintmp) 1003 do i=1,ncyc 1004 write(unit=lintmp,fmt='(f14.6)') ecyc(i) 1005 call write_files(iout,0,lintmp) 1006 enddo 1007 1008 25 continue 1009 1010C The program only uses those steps that have a progress line at the end 1011C of them (one with NSERCH, etc.), so nstep above may not be the nstep 1012C used later (see line 30) 1013 1014 rewind(inp) 1015 nstep=0 1016 11 continue 1017 call search(inp,line,' RMS of gradient ',ifound,0) 1018 if (ifound .eq. 1) then 1019 nstep=nstep+1 1020 goto 11 1021 endif 1022 write(*,*) 'n RMS=',nstep 1023 1024 rewind(inp) 1025 do i=1,nstep 1026 call search(inp,line,' RMS of gradient',ifound,0) 1027 if (ifound .eq. 0) goto 30 1028 call backs(inp,line,' Molecular Geometry',ifound,0) 1029 call skip(inp,2) 1030 do j=1,nat 1031 call read_line(inp,line) 1032 read(unit=line,fmt='(6x,a8,2x,3f20.10)') atmnam(j), 1033 1 geom(1,j,i),geom(2,j,i), geom(3,j,i) 1034 write(*,*) j,geom(1,j,i),geom(2,j,i), geom(3,j,i) 1035 do k=1,3 1036 geom(k,j,i) = geom(k,j,i)*bohr2a 1037 enddo 1038 enddo 1039 call search(inp,line,' Current energy ',ifound,-1) 1040 read(unit=line,fmt='(32x,f19.10)') estep(i) 1041 write(*,*) 'estep(i)',i,estep(i) 1042 call skip(inp,2) 1043 call read_line(inp,line) 1044 read(unit=line,fmt='(32x,f13.6)') rmsfor(i) 1045 write(*,*) 'rmsfor(i)',rmsfor(i) 1046 call read_line(inp,line) 1047 read(unit=line,fmt='(32x,f13.6)') xmxfor(i) 1048 write(*,*) 'xmxfor(i)',xmxfor(i) 1049 enddo 1050 nstep=i 1051 write(*,*) 'nstep',nstep 1052 1053C Now write everything out 1054 1055 30 nstep=i-1 1056 call write_files(iout,0,'[GEOCONV]') 1057 call write_files(iout,0,'energy') 1058 do i=1,nstep 1059 write(unit=lintmp,fmt='(1x,f13.6)') estep(i) 1060 call write_files(iout,0,lintmp) 1061 enddo 1062 call write_files(iout,0,'max-force') 1063 do i=1,nstep 1064 write(unit=lintmp,fmt='(4x,f10.6)') xmxfor(i) 1065 call write_files(iout,0,lintmp) 1066 enddo 1067 call write_files(iout,0,'rms-force') 1068 do i=1,nstep 1069 write(unit=lintmp,fmt='(4x,f10.6)') rmsfor(i) 1070 call write_files(iout,0,lintmp) 1071 enddo 1072 1073C Note that the program prints out only the Cartesian geometries 1074C (no matter what coordinates were used in the optimization) 1075 1076 call write_files(iout,0,'[GEOMETRIES] XYZ') 1077 write(unit=temp8,fmt='(i5)') nat 1078 do i=1,nstep 1079 write(unit=lintmp,fmt='(a9,f13.6)') 'scf done:', estep(i) 1080 call write_files(iout,0,temp8) 1081 call write_files(iout,0,lintmp) 1082 do j=1,nat 1083 write(unit=lintmp,fmt='(1x,a8,3(f11.6,2x))') atmnam(j), 1084 1 geom(1,j,i),geom(2,j,i), geom(3,j,i) 1085 call write_files(iout,0,lintmp) 1086 enddo 1087 enddo 1088 write(6,*) 'Finished writing geometry optimization information!' 1089 return 1090 end 1091C----------------------------------------------------------------------- 1092 subroutine find_freq(inp,iout,itmp) 1093 1094C This subroutine finds and writes the frequencies and normal modes of 1095C vibration. 1096 1097 implicit double precision(a-e,g-h,o-z) 1098 implicit logical(f) 1099 common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4 1100 common /flags2/ flgmc 1101 common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb 1102 character line*130,lintmp*130 1103 parameter(maxat=100) 1104 dimension efreq(3*maxat),vibcor(3*maxat,3*maxat) 1105 1106 write(6,*) 'Reading vibrational frequencies and displacements...' 1107 call write_files(iout,0,'[FREQ]') 1108 rewind(inp) 1109 call search(inp,line,'Transformation matrix from',ifound,-1) 1110 j=1 1111 10 k=j+5 1112 call skip(inp,3) 1113 call read_line(inp,line) 1114 read(unit=line,fmt='(17x,6f10.3)') (efreq(i), i=j,k) 1115 write(*,fmt='(6(i3,f9.3))') (i,efreq(i), i=j,k) 1116 call skip(inp,2) 1117 do l=1,3*nat 1118 call read_line(inp,line) 1119 read(unit=line,fmt='(17x,6f10.5)') (vibcor(i,l),i=j,k) 1120* write(*,fmt='(6f10.5)') (vibcor(i,l),i=j,k) 1121 enddo 1122 if (k .ge. 3*nat) then 1123 do i=1,3*nat 1124 write(unit=lintmp,fmt='(f10.4)') efreq(i) 1125 call write_files(iout,0,lintmp) 1126 enddo 1127 call write_files(iout,0,'[FR-COORD]') 1128 rewind(itmp) 1129 call skip(itmp,2) 1130 do i=1,nat 1131 call read_line(itmp,line) 1132 read(unit=line,fmt='(1x,a8,6x,3(2x,f12.7))') atmnam,x,y,z 1133 write(unit=lintmp,fmt='(a8,3(1x,f12.6))') atmnam,x,y,z 1134 call write_files(iout,0,lintmp) 1135 enddo 1136 call write_files(iout,0,'[FR-NORM-COORD]') 1137 do i=1,3*nat 1138 write(unit=lintmp,fmt='(a9,i5)') 'vibration',i 1139 call write_files(iout,0,lintmp) 1140 do j=1,3*nat,3 1141 write(unit=lintmp,fmt='(3(f12.6,1x))') vibcor(i,j), 1142 1 vibcor(i,j+1),vibcor(i,j+2) 1143 call write_files(iout,0,lintmp) 1144 enddo 1145 enddo 1146 goto 20 1147 endif 1148 j=j+6 1149 goto 10 1150 20 write(6,*) 'Finished writing vibrational frequencies and 1151 1 displacements!' 1152 return 1153 end 1154C---------------------------------------------------------------------- 1155 subroutine close_files(inp,iout,itmp) 1156 1157C Closes all open files (temp file is deleted) 1158 1159 write(6,*) 'Closing files...' 1160 close(unit=inp, err=10, status='KEEP') 1161 close(unit=iout, err=20, status='KEEP') 1162 close(unit=itmp, err=30, status='DELETE') 1163 goto 40 1164 10 write(6,*) 'An error occurred while closing the input file!' 1165 call error 1166 20 write(6,*) 'An error occurred while closing the output file!' 1167 call error 1168 30 write(6,*) 'An error occurred while closing the temp file!' 1169 call error 1170 40 write(6,*) 'All done! Ending program...' 1171 return 1172 end 1173C----------------------------------------------------------------------- 1174 subroutine blank_line(iout) 1175 1176C Writes a blank line to unit=iout. 1177 1178 write(iout,fmt='(a)') ' ' 1179 return 1180 end 1181C---------------------------------------------------------------------- 1182 subroutine error 1183 write(6,*) 'An error has occurred while processing the file.' 1184 write(6,*) 'Aborting run...' 1185 stop 1186 end 1187C----------------------------------------------------------------------- 1188 subroutine ffwrd(iunit) 1189 1190C Fast-forward unit=iunit to the end of the file 1191 1192 10 read(iunit,fmt='(a130)',end=20) line 1193 goto 10 1194 20 return 1195 end 1196C---------------------------------------------------------------------- 1197 subroutine last_char(line,length,loc) 1198 character*(*) line 1199 1200C Given string LINE of length LENGTH locate and return in LOC 1201C the position of the last character. 1202 1203 loc=length 1204 do i=length,1,-1 1205 if(line(i:i).ne.' ') then 1206 loc=i 1207 return 1208 endif 1209 enddo 1210 return 1211 end 1212C----------------------------------------------------------------------- 1213 subroutine read_line(iunit,line) 1214 character line*130 1215 read(iunit,fmt='(a130)') line 1216 return 1217 end 1218C----------------------------------------------------------------------- 1219 subroutine search(iunit,line,text,ifound,idbug) 1220 character*(*) text 1221 character line*130 1222 1223C Scan unit IUNIT for the next occurence of TEXT. 1224C If idbug > 0 write line to unit idbug (probably 6). 1225C If idbug = -1 then stop program if TEXT not found 1226 1227 ifound=0 1228 l=len(text) 1229 10 read(iunit,fmt='(a130)',end=20) line 1230 if (idbug .gt. 0) then 1231 call write_files(idbug,0,line) 1232 endif 1233 do i=1,130-l+1 1234 if (line(i:i+l-1) .eq. text) then 1235 ifound=1 1236 return 1237 endif 1238 enddo 1239 goto 10 1240 20 if (idbug .eq. 1) then 1241 write(6,*) 'Cannot find ', text(1:l+1) 1242 call error 1243 endif 1244 return 1245 end 1246C----------------------------------------------------------------------- 1247 subroutine skip(iunit,nlines) 1248 character char*1 1249 1250C Skip next NLINES on unit IUNIT 1251 1252 do i=1,nlines 1253 read(iunit,fmt='(1a1)') char 1254 enddo 1255 return 1256 end 1257C----------------------------------------------------------------------- 1258 subroutine write_files(iout1,iout2,text) 1259 character*(*) text 1260 1261C Writes TEXT to files iout1 and iout2 (usually output and temp files) 1262C or just output file (iout2=0) 1263 1264 loc1=len(text) 1265 call last_char(text,loc1,loc) 1266 write(iout1,'(a)') text(1:loc) 1267 if (iout2 .eq. 0) goto 10 1268 write(iout2,'(a)') text(1:loc) 1269 1270 10 return 1271 end 1272C---------------------------------------------------------------------- 1273 1274 subroutine up_case(string,lenstr) 1275c convert all lower case characters to upper case 1276C 1277 CHARACTER*(*) STRING 1278 CHARACTER*26 UCASE,LCASE 1279C 1280 DATA UCASE /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ 1281 DATA LCASE /'abcdefghijklmnopqrstuvwxyz'/ 1282C 1283C CONVERTS LOWER CASE TO UPPER CASE IN THE GIVEN STRING. 1284C 1285 DO 100 I=1,LENSTR 1286 IC = INDEX(LCASE,STRING(I:I)) 1287 IF (IC.GT.0) STRING(I:I) = UCASE(IC:IC) 1288 100 CONTINUE 1289 RETURN 1290 END 1291C----------------------------------------------------------------------- 1292 subroutine backs(iunit,line,text,ifound,idbug) 1293 character*(*) text 1294 character line*130 1295 1296C Scan unit IUNIT backwards for the occurence of TEXT. 1297C If idbug > 0 write line to unit idbug (probably 6). 1298C If idbug = -1 then stop program if TEXT not found 1299 1300 ifound=0 1301 l=len(text) 1302 10 read(iunit,fmt='(a130)',end=20) line 1303 if (idbug .gt. 0) then 1304 call write_files(idbug,0,line) 1305 endif 1306 do i=1,130-l+1 1307 if (line(i:i+l-1) .eq. text) then 1308 ifound=1 1309 return 1310 endif 1311 enddo 1312 backspace(iunit) 1313 backspace(iunit) 1314 goto 10 1315 20 if (idbug .eq. 1) then 1316 write(6,*) 'Cannot find ', text(1:l+1) 1317 call error 1318 endif 1319 return 1320 end 1321