1 subroutine rdgdud(idebug,ibefo,istatio,irtype,istats, 2 & focc,focb,nocc,nocb,ncols,ncolb,coo,ianz) 3 4c THIS IS REALLY rdgaus 5 6 implicit double precision ( a-h,o-z) 7 parameter (numat1=20000) 8 parameter (mxonh=100) 9 common /athlp/ iatoms, mxnat 10 common /qmchar/ qch(numat1),ihasesp 11 parameter (numatm=2000) 12 parameter (mxel=100) 13 character*137 line,lint 14 character*2 els 15 common /coord / xyz(3,numatm) 16 common /pseudo /ipseud,ivale(numatm) 17 common /moldat/ natoms, norbs, nelecs,nat(numatm) 18 common /orbhlp/ mxorb,iuhf,ispd 19 common /gauori/ nzm,nso,nio,nzo,ioropt,ifor, 20 & ixyz98,iopr,isymm,irc,imp2,icntp,itd 21 common /gauver/ ivers 22 common /rdwr/ iun1,iun2,iun3,iun4,iun5 23 character*2 elemnt 24 common /elem/elemnt(mxel) 25 common /oniomh/ ioni,nion,ionih(mxonh),natonh(mxonh), 26 & iomap(numatm),icntat(mxonh),fct(mxonh), 27 & xyzi(3,mxonh) 28 common /nmr/ shlnuc(numatm),ihsnmr 29 common /uvspec/ ihasex 30 logical rohf 31 dimension focc(*),focb(*),v1(3) 32 dimension coo(3,*),ianz(*) 33 34 istats = 1 35 irtype = 0 36 istatio = 0 37 isymm = 1 38 rohf = .false. 39 ig94 = 0 40 ixyz98 = 0 41 iopr = 0 42 irc = 0 43 imp2 = 0 44 itd = 0 45 icntp = 0 46 ihasex = 0 47 48 call search(line,' ***********************',istat) 49 call nxtlin(line,jstat) 50 if (line(11:12).eq.'DV') then 51 ivers = 2003 52 else 53 read(line,'(10x,i2)') ivers 54 if (ivers.eq.3) ivers = 2003 55 if (ivers.eq.9) ivers = 2009 56 if (ivers.eq.16) ivers = 2016 57 if (ivers.ge.94) ig94 = 1 58 endif 59 if (ivers.ge.98) then 60 ixyz98 = 1 61 call search(line,'will use Cartesian coordinates',istat) 62 if (istat.ne.0) ixyz98 = 2 63 call rewfil 64 endif 65 66 if (ig94.eq.0) then 67 call search(line,'Restricted open shell SCF:',istat) 68 if (istat.ne.0) rohf = .true. 69 call rewfil 70 endif 71 72 call searchd(line,'Frequencies --', 73 & 'Excitation energies and oscillator strengths:', 74 & istat) 75 if (istat.ne.0) then 76 if (index(line,'Frequencies --').ne.0) irtype = 4 77 if (index(line,'Excitation').ne.0) then 78 ihasex = 1 79 call gttrns(ist) 80 endif 81 endif 82 call rewfil 83 84 call search(line,'Symmetry turned off',istat) 85 if (istat.ne.0) isymm = 0 86 call rewfil 87 88 call search(line,'IRC-IRC-IRC',istat) 89 if (istat.ne.0) then 90 irc = 1 91 call search(line,'Integration scheme',istat) 92 if (istat.ne.0) then 93 if (icdex(line,'HPC').ne.0) irc = 2 94 endif 95 endif 96 call rewfil 97 98 call searchd(line,'Computing MP2 derivatives', 99 & 'Computing MP2/KS-MP2 derivatives',istat) 100 if (istat.ne.0) imp2 = 1 101 call rewfil 102 103 call search(line,'Computing CIS/TD-HF/TD-KS derivatives',istat) 104 if (istat.ne.0) itd = 1 105 call rewfil 106 107 108 call search(line,'Counterpoise: corrected',istat) 109 if (istat.ne.0) icntp = 1 110 call rewfil 111 112 call searchd(line,'ONIOM: saving','ONIOM: Cut',istat) 113 if (istat.ne.0) then 114 ioni = 1 115 nion = 0 116 do while (index(line,'Cut').ne.0) 117 if (nion.lt.mxonh) then 118 nion = nion + 1 119 id = index(line,'/') + 1 120 els = line(id:id+1) 121 if (els(2:2).eq.' ') then 122 els = ' '//els(1:1) 123 endif 124 do i=1,99 125 if (els.eq.elemnt(i)) natonh(nion) = i 126 end do 127 if (line(id+8:id+8).eq.' ') then 128 read(line(id+2:id+38),'(i6,7x,i6,8x,f10.6)') 129 & ionih(nion),icntat(nion),fct(nion) 130 else if (line(id+7:id+7).eq.' ') then 131 read(line(id+2:id+35),'(i5,7x,i5,8x,f9.6)') 132 & ionih(nion),icntat(nion),fct(nion) 133 else 134 print*,'error reading ONIOM: Cut atoms' 135 endif 136 endif 137 call redel(line,1) 138 end do 139 if (idebug.eq.1) then 140 print*,'Oniom boundary atoms:' 141 do i=1,nion 142 print*,natonh(i),' ',ionih(i) 143 end do 144 endif 145 else 146 ioni = 0 147 endif 148 149 if (ioni.eq.1) then 150 call searchd(line,'high level on model system', 151 & 'generating new system at layer 1',istat) 152 if (istat.eq.0) then 153 call inferr('no high level on model system found!',1) 154 goto 1000 155 endif 156 endif 157c 158c look for number of orbitals and electrons 159c 160 call search(line,'primitive gaussians',istat) 161 if (istat.eq.0) then 162 call inferr('no primitive gaussian found!',1) 163 goto 1000 164 endif 165 if (index(line,'primitive gaussians').eq.30) then 166 read(line,'(1x,i3)',err=1000,end=1000) norbs 167 else 168 read(line,'(1x,i5)',err=1000,end=1000) norbs 169 endif 170 if (norbs.gt.mxorb) 171 & call inferr('Exceeding MaxNum of Orbitals!',1) 172 call search(line,'alpha electrons',istat) 173 if (istat.eq.0) then 174 call inferr('no alpha electrons found!',1) 175 goto 1000 176 endif 177 if (index(line,'alpha electrons').eq.6) then 178 read(line,'(1x,i3,21x,i4)',err=1000,end=1000) neleca,nelecb 179 else 180 read(line,'(1x,i5,20x,i5)',err=1000,end=1000) neleca,nelecb 181 endif 182 nelecs = neleca + nelecb 183 iorbs = 0 184 call search(line,'One-electron integrals computed using',istat) 185 if (istat.eq.0) then 186 call rewfil 187 call search(line,'alpha electrons',istat) 188 endif 189 do i=1,30 190 call nxtlin(line,jstat) 191 if (jstat.eq.1) goto 999 192 if (index(line,'RelInt:').eq.0) then 193 if (index(line,'NBsUse=').ne.0) then 194 read(line,'(9x,i5)',err=1000,end=1000) iorbs 195 endif 196 endif 197 end do 198999 continue 199 if (ig94.eq.1) then 200 ncols = norbs 201 ncolb = norbs 202 if (iorbs.ne.0) then 203 ncols = iorbs 204 ncolb = iorbs 205 endif 206 else 207 ncols = max0(neleca+5,nelecb+5) 208 ncols = min0(ncols,norbs) 209 if (rohf) ncols = norbs 210 ncolb = ncols 211 endif 212c 213 do i=1,min0(mxorb,norbs) 214 focc(i) = 0.0d0 215 end do 216 do i=1,min0(mxorb,neleca) 217 focc(i) = 1.0d0 218 end do 219 220 call rewfil 221 if (ioni.eq.1) then 222 call searchd(line,'high level on model system', 223 & 'generating new system at layer 1',istat) 224 if (istat.eq.0) then 225 call inferr('no high level on model system found!',1) 226 goto 1000 227 endif 228 endif 229c call searchd(line,'UHF open shell SCF','E(UHF)',istat) 230 call search(line,'Beta Molecular Orbital Coefficients',istat) 231 if (istat.ne.0) then 232 iuhf=1 233 nocc=neleca 234 nocb=nelecb 235 do i=1,min0(mxorb,norbs) 236 focb(i) = 0.0d0 237 end do 238 do i=1,min0(mxorb,nelecb) 239 focb(i) = 1.0d0 240 end do 241 else 242 nocc=max0(neleca,nelecb) 243 do i=1,min0(mxorb,nelecb) 244 focc(i) = focc(i) + 1.0d0 245 end do 246 endif 247c 248 call search(line,'-- Stationary point found',istat) 249 call rewfil 250 ihssta = istat 251 252c if (ivers.ge.2009) then 253c ibefo = 1 254c write(iun3,*) 255c & 'Warning G09: using orbitals input geometry!' 256c write(iun3,*) ' ' 257c endif 258 259 if (istat.eq.0.or.ibefo.eq.1.or.(istat.eq.1.and.irtype.eq.4)) 260 &then 261 write(iun3,*)'First standard orientation encountered used' 262 write(iun3,*)'for density calculations' 263 call search(line,'Z-MATRIX (ANGSTROMS AND DEGREES)',istat) 264 if (istat.eq.0) then 265 call haszm(.false.) 266 else 267 call redel(line,2) 268 if (irtype.eq.4) then 269 call convzmat(coo,ianz,iatoms,1,0,1) 270 else 271 call convzmat(coo,ianz,iatoms,1,1,1) 272 endif 273 endif 274 call rdcor(idebug,istat) 275 if (istat.eq.0) goto 1000 276 if (idebug.eq.1) 277 & call inferr('Succesfully read coordinates',0) 278 else 279 istatio = 1 280 281 write(iun3,*)'Coordinates of stationary point used' 282 write(iun3,*)'for density calculations' 283 28410 call rdcor(idebug,istat) 285 if (istat.eq.0) goto 20 286 if (idebug.eq.1) 287 & call inferr('found coordinates',0) 288 goto 10 28920 continue 290 endif 291c 292 call rewfil 293 call search(line,'Charges from ESP fit',istat) 294 if (ihssta.eq.1) call search(line,'Charges from ESP fit',istat) 295 if (istat.eq.1) then 296 call redel(line,2) 297 do i=1,natoms 298 call nxtlin(line,jstat) 299 if (jstat.eq.1.or.jstat.eq.2) goto 1000 300 if (linlen(line).eq.18) then 301 read(line,'(7x,f11.6)',err=1000,end=1000) qch(i) 302 else if (linlen(line).eq.21) then 303 read(line,'(11x,f11.6)',err=1000,end=1000) qch(i) 304 endif 305 end do 306 ihasesp = 1 307 endif 308 309 call rewfil 310 call fndor(idebug) 311 312 if (ixyz98.gt.0.and.natoms.gt.50) then 313 314c check for IOP(2/11=1) 315 316 call rewfil 317 call search(line,' 2/',istat) 318 if (istat.eq.1) then 319 if (ichar(line(1:1)).eq.32.and. 320 & ichar(line(2:2)).eq.50.and. 321 & ichar(line(3:3)).eq.47) then 322 if (index(line,'11=1').ne.0) then 323 iopr = 1 324 endif 325 endif 326 endif 327 endif 328 329 if (ioni.eq.1) then 330 331c coordinates are read, calculate the coordinates of the oniom link 332c atoms 333 334 do j=1,nion 335 do k=1,3 336 v1(k) = xyz(k,ionih(j))-xyz(k,icntat(j)) 337 end do 338 do k=1,3 339 xyzi(k,j) = xyz(k,icntat(j)) + 340 & fct(j)*v1(k) 341 end do 342 end do 343 call xyzcoo(1,0,0) 344 345 endif 346 347 call rewfil 348 if (ioni.eq.1) then 349 call searchd(line,'high level on model system', 350 & 'generating new system at layer 1',istat) 351 if (istat.eq.0) then 352 call inferr('no high level on model system found!',1) 353 goto 1000 354 endif 355 endif 356 357c if ioni=1, rdbas transforms the xyz array to those of H layer only 358 359 call rdbas(idebug,0,istat) 360 call norml 361 362 if (istat.eq.0) goto 1000 363 if (idebug.eq.1) 364 & call inferr('Succesfully read basis-set',0) 365 if (idebug.eq.1) call basprt(iun3,.false.,.false.) 366 if (norbs.gt.mxorb) goto 1000 367 368 call search(line,'Pseudopotential Parameters',istat) 369 if (istat.eq.1) then 370 ipseud = 1 371 call redel(line,4) 372 do i=1,natoms 373 call nxtlin(line,jstat) 374 call nxtlin(lint,jstat) 375 if (icdex(lint,'No pseudo').ne.0) then 376 read(line,'(14x,i3)') ivale(i) 377 else 378 read(line,'(26x,i3)') ivale(i) 379 do while(.true.) 380 call redel(line,1) 381 if (line(5:5).ne.' ') then 382 call bckfil 383 goto 100 384 endif 385 end do 386 endif 387100 continue 388 end do 389 endif 390 391 natbck = natoms 392 if (ioni.eq.1) then 393 394c rdbas sets the array iomap, when ioni=1 395 396 natmp = 0 397 do i=1,numatm 398 if (iomap(i).ne.0) natmp = natmp + 1 399 end do 400 natoms = natmp 401 endif 402 403 call rewfil 404 if (ioni.eq.1.and.istatio.eq.1) then 405 call search(line,'-- Stationary point found',istat) 406 if (istat.eq.0) then 407 call inferr('no Stationary point found!',1) 408 goto 1000 409 endif 410 endif 411 412 call nmrshl 413 if (ihsnmr.eq.2) call nmrcpl(idebug) 414 call rewfil 415 416 if (ioni.eq.1) then 417 call searchd(line,'high level on model system', 418 & 'generating new system at layer 1',istat) 419 if (istat.eq.0) then 420 call inferr('error reading vectors: '// 421 & 'no high level on model system found for Stationary point!' 422 & ,1) 423 call inferr( 424 & 'no vectors for high level on model system '// 425 & 'of stationary point: try molden -b',1) 426 goto 1000 427 endif 428 endif 429 if (istatio.eq.1) then 430 call rdvect(idebug,ig94,istat) 431 istato = istat 432 433 if (ioni.ne.1) then 434 435 if (ivers.lt.2009) then 436 437 call rdvect(idebug,ig94,istat) 438 if (istat.eq.0.and.istato.eq.0) goto 1000 439 440 else 441 call search(line,'Population analysis using',istat) 442 if (istat.eq.1) then 443 call rdvect(idebug,ig94,istat) 444 if (istat.eq.0.and.istato.eq.0) goto 1000 445 endif 446 endif 447 448 else if (ioni.eq.1.and.istato.eq.0) then 449 450 call inferr( 451 & 'no vectors for high level on model system '// 452 & 'of stationary point: use molden -b',1) 453 goto 1000 454 455 endif 456 457 else 458 459 call rdvect(idebug,ig94,istat) 460 if (istat.eq.0) goto 1000 461 462 endif 463 464 if (idebug.eq.1) call inferr('Succesfully read vectors',0) 465 466 467 if (idebug.eq.1) 468 & call inferr('Succesfully read Gaussian output',0) 469 470 return 471 4721000 if (idebug.eq.1) then 473 call inferr('Error reading Gaussian output!',1) 474 print*,line 475 endif 476 477 istats = 0 478 return 479 480 end 481 482 subroutine cubtst(iun,ijag) 483 implicit double precision ( a-h,o-z) 484 character*137 line,str 485 common /curlin/ line 486 common /rdwr/ iun1,iun2,iun3,iun4,iun5 487 common /grdhlp/ mx3d,mx3d2 488 integer currec,bigend,buflen 489 common /rdrec/ iund,currec 490 common /isend/ bigend,nlen 491 character keywrd*320, keyori*320 492 common /keywrd/ keywrd,keyori 493 equivalence (bufft(1),nfast) 494 equivalence (bufft(5),nmed) 495 equivalence (bufft(9),nslow) 496 integer getlin 497 character*4 buff(40) 498 integer*2 buf2(256) 499 character*1 bufft(1024) 500 integer*2 idum 501 502 if (ijag.eq.1) then 503c Jaguar cube 504 iuntmp = iun2 505 iun2 = iun 506 507 do while(getlin(1).eq.1) 508 ktype = nxtwrd(str,nstr,itype,rtype) 509 if (ktype.eq.1) then 510 if (nstr.ge.4) then 511 if (icdex(str,'npts').ne.0) then 512 ktype = nxtwrd(str,nstr,itype,rtype) 513 if (ktype.eq.2) npts1 = itype 514 ktype = nxtwrd(str,nstr,itype,rtype) 515 if (ktype.eq.2) npts2 = itype 516 ktype = nxtwrd(str,nstr,itype,rtype) 517 if (ktype.eq.2) npts3 = itype 518 endif 519 endif 520 endif 521 end do 522 523 elseif (ijag.eq.2) then 524 525 iund = 10 526 bigend = 1 527 nlen = 1 528 currec = 0 529 buflen = 40 530 close(10) 531 53210 open(unit=iund,file=keywrd,access="direct", 533 & recl=nlen) 534 535 call getrec(buff,buflen,1,ierr) 536 if (ierr.eq.1) then 537 bigend = 0 538 currec = 0 539 rewind(iund) 540 call getrec(buff,buflen,1,ierr) 541 if (ierr.eq.1.and.nlen.eq.1) then 542 currec = 0 543 bigend = 1 544 nlen = 4 545 close(iund) 546 goto 10 547 endif 548 endif 549 550 if (ierr.eq.1) then 551 ijag = -1 552 return 553 endif 554 555 call byter(buff(26),npts1) 556 call byter(buff(27),npts2) 557 call byter(buff(28),npts3) 558 559 rewind(iund) 560 currec = 0 561 562 elseif (ijag.eq.3) then 563 564c O map (DSN6) 565 566 iund = 10 567 bigend = 1 568 nlen = 1 569 currec = -1 570 57112 open(unit=iund,file=keywrd(1:linlen(keywrd)),access="direct", 572 & recl=nlen) 573 574 call getrc2(buf2,ierr) 575 if (ierr.eq.1) then 576 577 bigend = 0 578 currec = -1 579 call getrc2(buf2,ierr) 580 if (ierr.eq.1.and.nlen.eq.1) then 581 582 currec = -1 583 bigend = 1 584 nlen = 2 585 close(iund) 586 goto 12 587 588 endif 589 590 endif 591 592 if (ierr.eq.1) then 593 ijag = -1 594 return 595 endif 596 597 npts1 = int(buf2(4)) 598 npts2 = int(buf2(5)) 599 npts3 = int(buf2(6)) 600 601c print*,'npts ',npts1,npts2,npts3 602 currec = -1 603 604 elseif (ijag.eq.4) then 605 606c ccp4 map 607 608 open(unit=iund,file=keywrd(1:linlen(keywrd)),access="stream") 609 610 read(iund,pos=1) bufft 611 612 if (ierr.eq.1) then 613 ijag = -1 614 return 615 endif 616 617 npts1 = nfast 618 npts2 = nmed 619 npts3 = nslow 620 621 currec = -1 622 623 else 624 625 read(iun,'(a)') line 626 read(iun,'(a)') line 627 read(iun,'(a)') line 628 629 read(iun,'(i5)') npts1 630 read(iun,'(i5)') npts2 631 read(iun,'(i5)') npts3 632 633 endif 634 635 nptsmx = npts1 636 if (npts2.gt.nptsmx) nptsmx = npts2 637 if (npts3.gt.nptsmx) nptsmx = npts3 638 if (ijag.eq.3) then 639 iptsmx = nptsmx/8 640 if (nptsmx.gt.iptsmx*8) nptsmx = (iptsmx+1)*8 641 endif 642 643 if (nptsmx.gt.mx3d) then 644 call allgrd(nptsmx) 645 endif 646 647 if (ijag.eq.1) then 648 rewind iun2 649 iun2 = iuntmp 650 else 651 rewind iun 652 endif 653 654 return 655 656100 print*,'error reading grid file' 657 return 658 end 659 660 subroutine rdcubd(npts1,npts2,npts3,iposng,ipsi,istat,iun, 661 & idebug,denn,pmnn) 662 implicit double precision ( a-h,o-z) 663 664c THIS IS REALLY rdcube 665 666 parameter (numatm=2000) 667 parameter (mxel=100) 668 common /coord / xyz(3,numatm) 669 common /moldat/ natoms, norbs, nelecs,nat(numatm) 670 common /grdhlp/ mx3d,mx3d2 671 common /plane/ px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat 672 common /eulx/ ca,cb,sa,sb,cc,sc 673 logical oclose,osshad,ofill,ofrst 674 common /spacom/ grey,iscol,oclose,osshad,ofill,ofrst 675 common /zproj/ zval(numatm),nindx(numatm),nconn(numatm) 676 common /elmcom/ vdwr(mxel),vrad(mxel),icol(mxel) 677 character*137 line,lstr,str 678 common /curlin/ line 679 common /rdwr/ iun1,iun2,iun3,iun4,iun5 680 logical gnreal,oampac,ojag 681 integer getlin 682 dimension orig(3),tmp(3),dtmp(6) 683 dimension denn(*),pmnn(*) 684 685 ojag = .false. 686 if (istat.eq.1) ojag = .true. 687 688 istat = 1 689 itype = 0 690 toang = 0.52917706d0 691 rneg = 0.0d0 692 693 if (ojag) then 694 695c Jaguar cube file 696 697 iuntmp = iun2 698 iun2 = iun 699 do while(.true.) 700 if (getlin(1).eq.1) then 701 if (icdex(line,'&end').ne.0) goto 10 702 ktype = nxtwrd(str,nstr,itype,rtype) 703 if (ktype.eq.1) then 704 if (nstr.eq.4) then 705 if (icdex(str,'npts').ne.0) then 706 ktype = nxtwrd(str,nstr,itype,rtype) 707 if (ktype.eq.2) npts1 = itype 708 ktype = nxtwrd(str,nstr,itype,rtype) 709 if (ktype.eq.2) npts2 = itype 710 ktype = nxtwrd(str,nstr,itype,rtype) 711 if (ktype.eq.2) npts3 = itype 712 endif 713 elseif (nstr.eq.6) then 714 if (icdex(str,'orig').ne.0) then 715 if (.not.gnreal(orig,3,.false.)) goto 100 716 endif 717 elseif (nstr.eq.7) then 718 if (icdex(str,'extent').ne.0) then 719 if (str(7:7).eq.'x'.or.str(7:7).eq.'X') then 720 if (.not.gnreal(v1,3,.false.)) goto 100 721 endif 722 if (str(7:7).eq.'y'.or.str(7:7).eq.'Y') then 723 if (.not.gnreal(v2,3,.false.)) goto 100 724 endif 725 if (str(7:7).eq.'z'.or.str(7:7).eq.'Z') then 726 if (.not.gnreal(tmp,3,.false.)) goto 100 727 endif 728 endif 729 endif 730 endif 731 endif 732 end do 73310 natoms = 0 734 nat(1) = 99 735 736 else 737 738c Gaussian Cube 739 740 read(iun,'(a)') line 741 read(iun,'(a)') line 742 if (index(line,'Density').ne.0) itype = 0 743 if (index(line,'MO coefficients').ne.0) itype = 1 744 ipsi = itype 745 746 read(iun,'(a)') line 747 oampac = (linlen(line).eq.44) 748 if (oampac) then 749 read(line,'(i5,3f13.6)') natoms,(orig(i),i=1,3) 750 else 751 read(line,'(i5,3f12.6)') natoms,(orig(i),i=1,3) 752 endif 753 if (natoms.lt.0) itype = 1 754 natoms = iabs(natoms) 755 756 if (oampac) then 757 read(iun,'(i5,3f13.6)') npts1,(v1(i),i=1,3) 758 read(iun,'(i5,3f13.6)') npts2,(v2(i),i=1,3) 759 read(iun,'(i5,3f13.6)') npts3,(tmp(i),i=1,3) 760 else 761 read(iun,'(i5,3f12.6)') npts1,(v1(i),i=1,3) 762 read(iun,'(i5,3f12.6)') npts2,(v2(i),i=1,3) 763 read(iun,'(i5,3f12.6)') npts3,(tmp(i),i=1,3) 764 endif 765 766 endif 767 768 if (npts1.gt.mx3d.or.npts2.gt.mx3d.or.npts3.gt.mx3d) then 769 lstr = 'dimension greater than maximum !' 770 print*,'npts1 ',npts1,' npts2 ',npts2,' npts3 ',npts3, 771 & ' maxdim ',mx3d 772 goto 100 773 endif 774 775 if (ojag) then 776 777 r(1) = vlen(v1) 778 r(2) = vlen(v2) 779 r(3) = vlen(tmp) 780 781 px = orig(1) + (v1(1) + v2(1) + tmp(1))/2.0d0 782 py = orig(2) + (v1(2) + v2(2) + tmp(2))/2.0d0 783 pz = orig(3) + (v1(3) + v2(3) + tmp(3))/2.0d0 784 785 else 786 787 r(1) = vlen(v1)*(npts1-1) 788 r(2) = vlen(v2)*(npts2-1) 789 r(3) = vlen(tmp)*(npts3-1) 790 791 px = orig(1) + ((npts1-1)*v1(1) + (npts2-1)*v2(1) + 792 & (npts3-1)*tmp(1))/2.0d0 793 py = orig(2) + ((npts1-1)*v1(2) + (npts2-1)*v2(2) + 794 & (npts3-1)*tmp(2))/2.0d0 795 pz = orig(3) + ((npts1-1)*v1(3) + (npts2-1)*v2(3) + 796 & (npts3-1)*tmp(3))/2.0d0 797 798 endif 799 800 call vsc1(v1,1.0d0,1.0d-4) 801 call vsc1(v2,1.0d0,1.0d-4) 802 call vsc1(tmp,1.0d0,1.0d-4) 803 804 805 cx = tmp(1) 806 cy = tmp(2) 807 cz = tmp(3) 808 809 if (.not.ojag) then 810 do i=1,natoms 811 if (oampac) then 812 read(iun,'(i5,4f13.6)') nat(i),tdum,(xyz(j,i),j=1,3) 813 else 814 read(iun,'(i5,4f12.6)') nat(i),tdum,(xyz(j,i),j=1,3) 815 endif 816 end do 817 endif 818 819 if (idebug.eq.1) then 820 print*,'r=',(r(i),i=1,3) 821 print*,'px py pz ',px,py,pz 822 print*,'cx cy cz ',cx,cy,cz 823 print*,'v1=',(v1(i),i=1,3) 824 print*,'v2=',(v2(i),i=1,3) 825 if (ojag) then 826 do i=1,natoms 827 print*, nat(i),(xyz(j,i),j=1,3) 828 end do 829 endif 830 endif 831 832 call crprod(v1,v2,dtmp) 833 call impsc(dtmp,tmp,cosb) 834 if (1.0d0-dabs(cosb).gt.0.001d0) then 835 lstr = 'Only rectangular grids supported !' 836 print*,'ERROR: Only rectangular grids supported ',cosb 837 goto 100 838 endif 839 840 if (.not.ojag.and.itype.eq.1) then 841 read(iun,'(i5)') no 842 if (no.gt.1) then 843 lstr = 'ONLY single orbital cube files are supported !' 844 goto 100 845 endif 846 nlines = (no + 1) / 10 847 if (no+1 - nlines*10.gt.0) nlines = nlines + 1 848 do i=1,nlines-1 849 read(iun,'(a)') line 850 end do 851 endif 852 853 ij = 0 854 iposng = 0 855 856 if (ojag) then 857 858c Jaguar cube file 859 860 do i=1,npts1 861 do j=1,npts2 862 ij = ij + 1 863 do k=1,npts3 864 865 if (getlin(0).eq.1) then 866 ktype = nxtwrd(str,nstr,itype,rtype) 867 if (ktype.eq.3) then 868 if (rtype.lt.0.0d0) iposng = 1 869 denn((k-1)*mx3d2 + ij) = rtype 870 endif 871 else 872 lstr = 'Not enough lines' 873 goto 100 874 endif 875 876 end do 877 end do 878 end do 879 880 else 881 882c Gaussian Cube 883 884 do i=1,npts1 885 do j=1,npts2 886 ij = ij + 1 887 kk = 0 888 do while (kk.lt.npts3) 889 n = 6 890 if (npts3-kk.lt.n) n = npts3-kk 891 read(iun,'(6E13.5)') (dtmp(k),k=1,n) 892 do k=1,n 893 denn((npts3-(kk+k))*mx3d2 + ij) = dtmp(k) 894 if (dtmp(k).lt.rneg) rneg = dtmp(k) 895 end do 896 kk = kk + 6 897 end do 898 end do 899 end do 900 endif 901 902c allow tolerance for density, since gaussian density cubes often 903c have small negative values 904 905 if (dabs(rneg).gt.1.0d08) iposng = 1 906 907 iplat = 0 908 909 ca = v2(2) 910 sa = -v2(1) 911 sb = -v1(3) 912 if (dabs(ca).gt.0.001d0) then 913 cb = v1(1) / ca 914 else 915 cb = v1(2) / sa 916 endif 917 918 cc = 1.0d0 919 sc = 0.0d0 920 921 922 call parrat 923 call proato 924 925 do i=1,npts3 926 pmnn(i) = 100000.d0 927 end do 928 929 do i=1,natoms 930 nconn(i) = 0 931 do j=1,natoms 932 if (i.ne.j) then 933 dmaxsq = (vdwr(nat(i)) + vdwr(nat(j)))**2 934 dijsq = ((xyz(1,i)-xyz(1,j))*toang)**2 935 & +((xyz(2,i)-xyz(2,j))*toang)**2 936 & +((xyz(3,i)-xyz(3,j))*toang)**2 937 if (dijsq.lt.dmaxsq) then 938 nconn(i) = nconn(i) + 1 939 endif 940 endif 941 end do 942 end do 943 944 call xyzcoo(1,0,0) 945 946 ofrst = .false. 947 948 if (ojag) iun2 = iuntmp 949 950 lstr = 'found cube file' 951 call inferr(lstr,0) 952 return 953 954100 continue 955110 istat = 0 956 if (ojag) iun2 = iuntmp 957 call inferr(lstr,0) 958 959 return 960 end 961 962 subroutine rdvasd(npts1,npts2,npts3,iposng,istat,lenf, 963 & idocub,idebug,denn,pmnn,bucket, 964 & coo,ianz,iatclr,iconn, 965 & nnat,norg,icent,inorm,ncon,nspg,ichx, 966 & nopr,ir,it, 967 & xa,ya,yb,za,zb,zc,a,b,c,alpha,beta,gamma) 968 implicit double precision ( a-h,o-z) 969 970c THIS IS REALLY rdvasp 971 972 parameter (numatm=2000) 973 parameter (mxel=100) 974 parameter (mxcon=10) 975 parameter (lnbuck=10) 976 common /coord / xyz(3,numatm) 977 common /moldat/ natoms, norbs, nelecs,nat(numatm) 978 common /grdhlp/ mx3d,mx3d2 979 common /plane/ px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat 980 common /eulx/ ca,cb,sa,sb,cc,sc 981 integer*2 ir,it 982 logical oclose,osshad,ofill,ofrst 983 common /spacom/ grey,iscol,oclose,osshad,ofill,ofrst 984 common /zproj/ zval(numatm),nindx(numatm),nconn(numatm) 985 common /elmcom/ vdwr(mxel),vrad(mxel),icol(mxel) 986 common /athlp/ iatoms, mxnat 987 character*2 elemnt 988 common /elem/elemnt(mxel) 989 990 common /rdwr/ iun1,iun2,iun3,iun4,iun5 991 character*137 line,lstr 992 logical opfil,gnreal 993 integer getlin 994 character*137 str 995 character*2 tocapf 996 character keywrd*320, keyori*320 997 common /keywrd/ keywrd,keyori 998 dimension tmp(3),dtmp(10),xyzt(3) 999 dimension ieltyp(20),ielnum(20) 1000 dimension denn(*),pmnn(*),bucket(*) 1001 dimension coo(3,*),ianz(*),iatclr(*),iconn(mxcon+1,*), 1002 & ir(3,3,192),it(3,192) 1003 1004 logical felem 1005 character*137 strtmp 1006 common /vasp/nheadx,natx 1007 1008 1009 istat = 1 1010 itype = 0 1011 toang = 0.52917706d0 1012 todeg = 45.0d0 / datan(1.0d0) 1013 lstr = ' ' 1014 1015 if (idocub.eq.1) then 1016 iuntmp = iun2 1017 iun2 = 21 1018 1019 if (.not.opfil(21,keywrd(1:lenf),lenf,1,1,0)) then 1020 lstr = keywrd(1:lenf)//' non-existent' 1021 goto 110 1022 endif 1023 endif 1024 1025 1026c title 1027 1028 natt = 0 1029 1030 if (idocub.eq.1) then 1031 call nxtlin(line,jstat) 1032 else 1033 1034 felem = .false. 1035 1036 if (getlin(0).eq.1) then 1037 do while(.true.) 1038 ktype = nxtwrd(str,nstr,itype,rtype) 1039 if (ktype.eq.1.and.(nstr.eq.1.or.nstr.eq.2)) then 1040 if (nstr.eq.1) then 1041 str(2:2) = str(1:1) 1042 str(1:1) = ' ' 1043 endif 1044 ity = 0 1045 do i=1,99 1046 if (tocapf(str).eq.tocapf(elemnt(i))) ity = i 1047 end do 1048 if (ity.ne.0) then 1049 felem = .true. 1050 natt = natt + 1 1051 ieltyp(natt) = ity 1052 endif 1053 else 1054 goto 10 1055 endif 1056 end do 1057 endif 1058 endif 1059 106010 continue 1061 1062c scale factor 1063 1064 if (getlin(0).eq.1) then 1065 ktype = nxtwrd(str,nstr,itype,rtype) 1066 if (ktype.eq.3) then 1067 scl = rtype 1068 endif 1069 endif 1070 1071c unit vectors 1072 1073 if (getlin(0).eq.1) then 1074 if (.not.gnreal(v1,3,.false.)) goto 100 1075 else 1076 goto 100 1077 endif 1078 1079 if (getlin(0).eq.1) then 1080 if (.not.gnreal(v2,3,.false.)) goto 100 1081 else 1082 goto 100 1083 endif 1084 1085 if (getlin(0).eq.1) then 1086 if (.not.gnreal(tmp,3,.false.)) goto 100 1087 else 1088 goto 100 1089 endif 1090 1091 r(1) = vlen(v1)*scl 1092 r(2) = vlen(v2)*scl 1093 r(3) = vlen(tmp)*scl 1094 if (r(1).le.0.0d0.or.r(2).le.0.0d0.or.r(3).le.0.0d0) goto 100 1095 1096 a = r(1) 1097 b = r(2) 1098 c = r(3) 1099 call impsc(v2,tmp,csa) 1100 call impsc(v1,tmp,csb) 1101 call impsc(v1,v2,csc) 1102 alpha = dacos(csa)*todeg 1103 beta = dacos(csb)*todeg 1104 gamma = dacos(csc)*todeg 1105 nspg = 1 1106 1107 call prcell(nspg,a,b,c,alpha,beta,gamma) 1108 call setop(xa,ya,yb,za,zb,zc,a,b,c,alpha,beta,gamma,1) 1109 call cprot(nspg,nopr,icent,ir,it,.false.) 1110 if (idebug.eq.1) call prop(nopr,ir,it) 1111 1112c In new POSCARs there might be a line here to give the elements names 1113 1114 if (getlin(0).eq.1) then 1115 1116 call gtplin(strtmp) 1117 ktype = nxtwrd(str,nstr,itype,rtype) 1118 1119 if (ktype.eq.1) then 1120 1121 if (.not.felem) then 1122 call setlin(strtmp,0) 1123 natt = 0 1124 do while(.true.) 1125 ktype = nxtwrd(str,nstr,itype,rtype) 1126 if (ktype.eq.1.and.(nstr.eq.1.or.nstr.eq.2)) then 1127 if (nstr.eq.1) then 1128 str(2:2) = str(1:1) 1129 str(1:1) = ' ' 1130 endif 1131 ity = 0 1132 do i=1,99 1133 if (tocapf(str).eq.tocapf(elemnt(i))) ity = i 1134 end do 1135 if (ity.ne.0) then 1136 felem = .true. 1137 natt = natt + 1 1138 ieltyp(natt) = ity 1139 endif 1140 else 1141 goto 11 1142 endif 1143 end do 1144 11 continue 1145 endif 1146 1147c we now read the number of atoms of each sort 1148 1149 natoms = 0 1150 ity = 0 1151 if (getlin(0).eq.1) then 1152 do i=1,50 1153 ktype = nxtwrd(str,nstr,itype,rtype) 1154 if (ktype.eq.2) then 1155 natoms = natoms + itype 1156 ity = ity + 1 1157 ielnum(ity) = itype 1158 else 1159 goto 20 1160 endif 1161 end do 1162 else 1163 goto 100 1164 endif 1165 1166 20 continue 1167 1168 elseif (ktype.eq.2) THEN 1169 1170c in fact, it is already the number of atoms of each sort 1171c we read them from the stored line 1172 1173 call setlin(strtmp,0) 1174 natoms = 0 1175 ity = 0 1176 do i=1,50 1177 ktype = nxtwrd(str,nstr,itype,rtype) 1178 if (ktype.eq.2) then 1179 natoms = natoms + itype 1180 ity = ity + 1 1181 ielnum(ity) = itype 1182 else 1183 goto 21 1184 endif 1185 end do 1186 1187 21 continue 1188 endif 1189 endif 1190 1191 natx = natoms 1192 1193c if (natt.eq.ity) then 1194 if (natt.ge.ity) then 1195 natt=ity 1196 k = 0 1197 do i=1,natt 1198 do j=1,ielnum(i) 1199 k = k + 1 1200 ianz(k) = ieltyp(i) 1201 nat(k) = ieltyp(i) 1202 end do 1203 end do 1204 endif 1205 1206c Check for Direct (=fractional coordinates) 1207 1208 if (getlin(0).eq.1) then 1209 ktype = nxtwrd(str,nstr,itype,rtype) 1210 if (ktype.eq.1) then 1211 call tocap(str,nstr) 1212 if (str(1:1).eq.'S') then 1213 if (getlin(0).eq.1) then 1214 ktype = nxtwrd(str,nstr,itype,rtype) 1215 if (ktype.ne.1) goto 100 1216 endif 1217 endif 1218 1219 if (str(1:1).ne.'C'.and.str(1:1).ne.'c'.and. 1220 & str(1:1).ne.'K'.and.str(1:1).ne.'k') then 1221 str = 'DIRECT' 1222 nstr = 6 1223 endif 1224 else 1225 goto 100 1226 endif 1227 end if 1228 1229c read in fractional coordinates 1230 1231 do i=1,natoms 1232 if (natt.ne.ity) then 1233 nat(i) = 99 1234 ianz(i) = 99 1235 endif 1236 if (getlin(0).eq.1) then 1237 if (.not.gnreal(xyzt,3,.false.)) goto 100 1238 do j=1,3 1239 xyz(j,i) = v1(j)*scl*xyzt(1) + v2(j)*scl*xyzt(2) + 1240 & tmp(j)*scl*xyzt(3) 1241 coo(j,i) = xyzt(j) 1242 end do 1243 else 1244 goto 100 1245 endif 1246 end do 1247 1248c if (idocub.eq.0) then 1249 1250 ncon = 1 1251 inorm = 0 1252 iatoms = natoms 1253 nnat = natoms 1254 norg = natoms 1255 nstor = mxnat-iatoms 1256 do i=1,iatoms 1257 1258 do j=1,3 1259 coo(j,nstor+i) = coo(j,i) 1260 end do 1261 1262 call fr2crt(coo(1,i),xa,ya,yb,za,zb,zc) 1263 do j=1,3 1264 coo(j,i) = coo(j,i)/toang 1265 end do 1266 ianz(nstor+i) = ianz(i) 1267 iatclr(nstor+i) = 1 1268 1269 end do 1270 1271 call doconn 1272 call dohcon(0) 1273 1274 do i=1,iatoms 1275 do j=1,iconn(1,i)+1 1276 iconn(j,nstor+i) = iconn(j,i) 1277 end do 1278 end do 1279 1280 1281c endif 1282 1283 ichx = 1 1284 1285 if (idocub.eq.0) then 1286 istat = 2 1287 return 1288 endif 1289 1290 call nxtlin(line,jstat) 1291 1292c x,y,z dimension of grid 1293 1294 if (getlin(0).eq.1) then 1295 ktype = nxtwrd(str,nstr,itype,rtype) 1296 if (ktype.eq.2) then 1297 npts1 = itype 1298 else 1299 goto 100 1300 endif 1301 ktype = nxtwrd(str,nstr,itype,rtype) 1302 if (ktype.eq.2) then 1303 npts2 = itype 1304 else 1305 goto 100 1306 endif 1307 ktype = nxtwrd(str,nstr,itype,rtype) 1308 if (ktype.eq.2) then 1309 npts3 = itype 1310 else 1311 goto 100 1312 endif 1313 else 1314 goto 100 1315 endif 1316 1317 if (npts1.gt.mx3d.or.npts2.gt.mx3d.or.npts3.gt.mx3d) then 1318 lstr = 'dimension greater than maximum !' 1319 print*,'npts1 ',npts1,' npts2 ',npts2,' npts3 ',npts3, 1320 & ' maxdim ',mx3d 1321 goto 100 1322 endif 1323 1324 px = (scl*v1(1) + scl*v2(1) + scl*tmp(1))/2.0d0 1325 py = (scl*v1(2) + scl*v2(2) + scl*tmp(2))/2.0d0 1326 pz = (scl*v1(3) + scl*v2(3) + scl*tmp(3))/2.0d0 1327 1328c do i=1,natoms 1329c xyz(1,i) = xyz(1,i) + px 1330c xyz(2,i) = xyz(2,i) + py 1331c xyz(3,i) = xyz(3,i) + pz 1332c end do 1333 1334 call vsc1(v1,1.0d0,1.0d-4) 1335 call vsc1(v2,1.0d0,1.0d-4) 1336 call vsc1(tmp,1.0d0,1.0d-4) 1337 1338 1339 cx = tmp(1) 1340 cy = tmp(2) 1341 cz = tmp(3) 1342 1343 if (idebug.eq.1) then 1344 print*,'r=',(r(i),i=1,3) 1345 print*,'px py pz ',px,py,pz 1346 print*,'cx cy cz ',cx,cy,cz 1347 print*,'v1=',(v1(i),i=1,3) 1348 print*,'v2=',(v2(i),i=1,3) 1349 do i=1,natoms 1350 print*, nat(i),(xyz(j,i),j=1,3) 1351 end do 1352 endif 1353 1354 call crprod(v1,v2,dtmp) 1355 call impsc(dtmp,tmp,cosb) 1356 if (1.0d0-dabs(cosb).gt.0.001d0) then 1357 lstr = 'Only rectangular grids supported !' 1358 print*,'ERROR: Only rectangular grids supported ',cosb 1359 goto 100 1360 endif 1361 1362 j = 1 1363 k = 1 1364 ijkt = 0 1365 iposng = 0 1366 nall = npts1*npts2*npts3 1367 ib = 0 1368 nb = (npts1/lnbuck)*lnbuck 1369 if (npts1.gt.nb) nb = nb + lnbuck 1370 1371 do while (ijkt.lt.nall) 1372c fill bucket 1373c ijkt = ijk 1374 do while (ib.lt.nb.and.ijkt.lt.nall) 1375 n = lnbuck 1376 if (nall-ijkt.lt.n) n = nall-ijkt 1377c print*,'n=',n,' bucket',bucket(ib),' ',nall,' ',ijkt 1378 read(21,'(10f8.3)') (bucket(ib+l),l=1,n) 1379 ib = ib + n 1380 ijkt = ijkt + n 1381 end do 1382c full bucket 1383 ib = ib - npts1 1384c fill density array 1385 do i=1,npts1 1386 denn((k-1)*mx3d2 + (i-1)*npts2+j) = bucket(i) 1387 end do 1388c print*,j,' ',k 1389c print*,(bucket(l),l=1,npts1) 1390c move access bucket to beginning of bucket 1391 do l=1,ib 1392 bucket(l) = bucket(npts1+l) 1393 end do 1394 j = j + 1 1395 if (j.gt.npts2) then 1396 j = 1 1397 k = k + 1 1398 endif 1399c ijk = ijk + npts1 1400 end do 1401 1402 iplat = 0 1403 1404 ca = v2(2) 1405 sa = -v2(1) 1406 sb = -v1(3) 1407 if (dabs(ca).gt.0.001d0) then 1408 cb = v1(1) / ca 1409 else 1410 cb = v1(2) / sa 1411 endif 1412 1413 cc = 1.0d0 1414 sc = 0.0d0 1415 1416 1417 call parrat 1418 call proato 1419 1420 do i=1,npts3 1421 pmnn(i) = 100000.d0 1422 end do 1423 1424 do i=1,natoms 1425 nconn(i) = 0 1426 do j=1,natoms 1427 if (i.ne.j) then 1428 dmaxsq = (vdwr(nat(i)) + vdwr(nat(j)))**2 1429 dijsq = ((xyz(1,i)-xyz(1,j))*toang)**2 1430 & +((xyz(2,i)-xyz(2,j))*toang)**2 1431 & +((xyz(3,i)-xyz(3,j))*toang)**2 1432 if (dijsq.lt.dmaxsq) then 1433 nconn(i) = nconn(i) + 1 1434 endif 1435 endif 1436 end do 1437 end do 1438 1439 if (idocub.eq.1) then 1440 close(21) 1441 call xyzcoo(1,0,0) 1442 endif 1443 1444 ofrst = .false. 1445 1446 if (idocub.eq.1) iun2 = iuntmp 1447 1448 lstr = 'found cube file' 1449 call inferr(lstr,0) 1450 return 1451 1452100 if (idocub.eq.1) close(21) 1453110 if (idocub.eq.1) iun2 = iuntmp 1454 istat = 0 1455 call inferr(lstr,0) 1456 1457 return 1458 end 1459 1460 subroutine wrvasd(iun,coo,ianz, 1461 & xa,ya,yb,za,zb,zc,a,b,c,alpha,beta,gamma) 1462 implicit double precision ( a-h,o-z) 1463 parameter (mxel=100) 1464 parameter (mxetyp=20) 1465 common /athlp/ iatoms, mxnat 1466 character*2 elemnt 1467 common /elem/elemnt(mxel) 1468 logical fndel 1469 dimension tmp(3),tr(3),rr(3,3),ieltyp(mxetyp),ielnum(mxetyp) 1470 dimension coo(3,*),ianz(*) 1471 1472 toang = 0.52917706d0 1473 1474 natt = 0 1475 do i=1,mxetyp 1476 ielnum(i) = 0 1477 end do 1478 1479 do i=1,iatoms 1480 1481 if (ianz(i).gt.0.and.ianz(i).lt.99) then 1482 1483 fndel = .false. 1484 do j=1,natt 1485 if (ianz(i).eq.ieltyp(j)) then 1486 ielnum(j) = ielnum(j) + 1 1487 fndel = .true. 1488 endif 1489 end do 1490 1491 if (.not.fndel) then 1492 natt = natt + 1 1493 ieltyp(natt) = ianz(i) 1494 ielnum(natt) = 1 1495 endif 1496 1497 endif 1498 1499 end do 1500 1501 write(iun,'(20(a2,1x))') (elemnt(ieltyp(i)),i=1,natt) 1502 1503 scl = 1.0d0 1504 write(iun,'(f12.6)') scl 1505 1506 call setrr(alpha,beta,gamma,a,b,c,rr) 1507 1508 tmp(1) = 1.0d0 1509 tmp(2) = 0.0d0 1510 tmp(3) = 0.0d0 1511c call fr2crt(tmp,xa,ya,yb,za,zb,zc) 1512 do k=1,3 1513 tr(k) = trc(tmp,rr,k) 1514 end do 1515 1516c write(iun,'(3(f12.6,1x))') (tmp(i),i=1,3) 1517 write(iun,'(3(f12.6,1x))') (tr(i),i=1,3) 1518 1519 tmp(1) = 0.0d0 1520 tmp(2) = 1.0d0 1521 tmp(3) = 0.0d0 1522c call fr2crt(tmp,xa,ya,yb,za,zb,zc) 1523 do k=1,3 1524 tr(k) = trc(tmp,rr,k) 1525 end do 1526c write(iun,'(3(f12.6,1x))') (tmp(i),i=1,3) 1527 write(iun,'(3(f12.6,1x))') (tr(i),i=1,3) 1528 1529 tmp(1) = 0.0d0 1530 tmp(2) = 0.0d0 1531 tmp(3) = 1.0d0 1532c call fr2crt(tmp,xa,ya,yb,za,zb,zc) 1533 do k=1,3 1534 tr(k) = trc(tmp,rr,k) 1535 end do 1536c write(iun,'(3(f12.6,1x))') (tmp(i),i=1,3) 1537 write(iun,'(3(f12.6,1x))') (tr(i),i=1,3) 1538 1539 write(iun,'(20(i3,1x))') (ielnum(i),i=1,natt) 1540 1541 write(iun,'(a)') 'Direct' 1542 1543 do k=1,natt 1544 1545 do i=1,iatoms 1546 if (ianz(i).eq.ieltyp(k)) then 1547 do j=1,3 1548 tmp(j) = coo(j,i) * toang 1549 end do 1550 call crt2fr(tmp,tr,xa,ya,yb,za,zb,zc) 1551 write(iun,'(3(f12.6,1x))') (tr(j),j=1,3) 1552 endif 1553 end do 1554 1555 end do 1556 1557 call inferr('Wrote file: POSCAR',0) 1558 1559 return 1560 end 1561 1562 subroutine wrmopd(iun,coo,ianz, 1563 & xa,ya,yb,za,zb,zc,a,b,c,alpha,beta,gamma) 1564 implicit double precision ( a-h,o-z) 1565 parameter (mxel=100) 1566 parameter (mxetyp=20) 1567 common /athlp/ iatoms, mxnat 1568 character*2 elemnt 1569 common /elem/elemnt(mxel) 1570 dimension tmp(3) 1571c dimension tr(3),rr(3,3) 1572 dimension coo(3,*),ianz(*) 1573 1574 toang = 0.52917706d0 1575 1576 write(iun,'(a)') 'AM1 NOINTER XYZ' 1577 write(iun,'(a)') ' ' 1578 write(iun,'(a)') ' ' 1579 1580c call setrr(alpha,beta,gamma,a,b,c,rr) 1581 1582 1583 1584 do i=1,iatoms 1585 do j=1,3 1586 tmp(j) = coo(j,i) * toang 1587 end do 1588 if (ianz(i).gt.0.and.ianz(i).lt.99) then 1589 if (i.eq.1) then 1590 write(iun,'(a2,1x,3(f12.6,a3))') elemnt(ianz(i)), 1591 & tmp(1),' 0 ',tmp(2),' 0 ',tmp(3),' 0 ' 1592 else 1593 write(iun,'(a2,1x,3(f12.6,a3))') elemnt(ianz(i)), 1594 & tmp(1),' 1 ',tmp(2),' 1 ',tmp(3),' 1 ' 1595 endif 1596 endif 1597 end do 1598 1599 1600 1601 tmp(1) = 1.0d0 1602 tmp(2) = 0.0d0 1603 tmp(3) = 0.0d0 1604 call fr2crt(tmp,xa,ya,yb,za,zb,zc) 1605c do k=1,3 1606c tr(k) = trc(tmp,rr,k) 1607c end do 1608 1609 write(iun,'(a3,3(f12.6,a3))') 1610 & 'Tv ',tmp(1),' 1 ',tmp(2),' 0 ',tmp(3),' 0 ' 1611c write(iun,'(3(f12.6,1x))') (tr(i),i=1,3) 1612 1613 tmp(1) = 0.0d0 1614 tmp(2) = 1.0d0 1615 tmp(3) = 0.0d0 1616 call fr2crt(tmp,xa,ya,yb,za,zb,zc) 1617c do k=1,3 1618c tr(k) = trc(tmp,rr,k) 1619c end do 1620c write(iun,'(3(f12.6,1x))') (tmp(i),i=1,3) 1621 write(iun,'(a3,3(f12.6,a3))') 1622 & 'Tv ',tmp(1),' 0 ',tmp(2),' 1 ',tmp(3),' 0 ' 1623c write(iun,'(3(f12.6,1x))') (tr(i),i=1,3) 1624 1625 tmp(1) = 0.0d0 1626 tmp(2) = 0.0d0 1627 tmp(3) = 1.0d0 1628 1629 call fr2crt(tmp,xa,ya,yb,za,zb,zc) 1630 1631c do k=1,3 1632c tr(k) = trc(tmp,rr,k) 1633c end do 1634c write(iun,'(3(f12.6,1x))') (tmp(i),i=1,3) 1635 1636 write(iun,'(a3,3(f12.6,a3))') 1637 & 'Tv ',tmp(1),' 0 ',tmp(2),' 0 ',tmp(3),' 1 ' 1638c write(iun,'(3(f12.6,1x))') (tr(i),i=1,3) 1639 1640 call inferr('Wrote file: mopac.dat',0) 1641 1642 return 1643 end 1644 1645 subroutine rdgrdd(npts1,npts2,npts3,iun,istat, 1646 & denn,dens,pmnn,buff) 1647 implicit double precision (a-h,o-z) 1648 parameter (mxdir=500) 1649 logical oclose,osshad,ofill,ofrst 1650 common /spacom/ grey,iscol,oclose,osshad,ofill,ofrst 1651 common /plane/ px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat 1652 common /grdhlp/ mx3d,mx3d2 1653 common /comsrf/ vo(3), rr(3),vv1(3),vv2(3),vv3(3), wo(3), 1654 & sl(3),isl 1655 common /cnthlp/ rng1,rng2,vlcnt,vlcnt2 1656 real ra,orig,den 1657 character*137 lstr 1658 character*4 buff(*) 1659 character*72 header 1660 integer currec,bigend,buflen 1661 common /rdrec/ iund,currec 1662 common /isend/ bigend,nlen 1663 dimension orig(3) 1664 dimension denn(*),pmnn(*),dens(*) 1665 1666c nnx,nny,nnz the number of points in x, y and z direction 1667c ra real step along x,y and z direction per point 1668c rx,ry,rz origin of the grid 1669 1670 toang = 0.52917706d0 1671 istat = 1 1672 lstr = ' ' 1673 1674 iplat = 0 1675 1676 buflen = mx3d2 1677 rewind(iund) 1678 call getrec(buff,buflen,0,ierr) 1679 1680 do i =1,18 1681 header(4*(i-1)+1:4*i) = buff(i)(1:4) 1682 end do 1683 1684 print*,'header = ',header 1685 1686 call byter(buff(26),npts2) 1687 call byter(buff(27),npts1) 1688 call byter(buff(28),npts3) 1689 1690 print*,'npts1=',npts1,' npts2=',npts2,' npts3=',npts3 1691 1692 call byter(buff(29),ra) 1693 1694 ra = ra / toang 1695 1696 print*,'ra=',ra 1697 1698 call byter(buff(30),orig(2)) 1699 call byter(buff(31),orig(1)) 1700 call byter(buff(32),orig(3)) 1701 1702 print*,'orgx,orgy,orgz ',orig(1),' ',orig(2),' ',orig(3) 1703 1704 v1(1) = 0.0d0 1705 v1(2) = 1.0d0 1706 v1(3) = 0.0d0 1707 1708 v2(1) = 1.0d0 1709 v2(2) = 0.0d0 1710 v2(3) = 0.0d0 1711 1712 r(1) = dble(ra)*dble(npts1-1) 1713 r(2) = dble(ra)*dble(npts2-1) 1714 r(3) = dble(ra)*dble(npts3-1) 1715 1716 do i=1,3 1717 orig(i) = orig(i) / toang 1718 rr(i) = r(i) 1719 vv1(i) = v1(i) 1720 vv2(i) = v2(i) 1721 vv3(i) = 0.0d0 1722 sl(i) = 1.0d0 1723 end do 1724 1725 vv3(3) = 1.0d0 1726 1727 isl = 1 1728 sl(1) = -1.0d0 1729 sl(2) = -r(2)/r(1) 1730 sl(3) = -r(3)/r(1) 1731 1732 call vnrm(sl) 1733 1734 px = dble(orig(2)) 1735 py = dble(orig(1)) 1736 pz = dble(orig(3)) 1737 1738 1739 vo(1) = px 1740 vo(2) = py 1741 vo(3) = pz 1742 1743 cx = 0.0d0 1744 cy = 0.0d0 1745 cz = 1.0d0 1746 1747 pmax = -1000000.d0 1748 pmin = 1000000.d0 1749 1750 do k=1,npts3 1751 1752c get plate 1753 1754 call getrec(buff,buflen,0,ierr) 1755 if (ierr.eq.1) goto 100 1756 1757 call byter(buff(1),nz) 1758 call byter(buff(2),nx) 1759 call byter(buff(3),ny) 1760 1761c print*,'nz,nx,ny ',nz,' ',nx,' ',ny 1762 1763 call getrec(buff,buflen,0,ierr) 1764 1765 ij = 0 1766 do i=1,ny 1767 do j=1,nx 1768 ij = ij + 1 1769 call byter(buff(ij),den) 1770 dens((j-1)*ny+i) = dble(den) 1771 end do 1772 end do 1773 do i=1,nx*ny 1774 denn((nz-1)*mx3d2 + i) = dens(i) 1775 if (dens(i).gt.pmax) pmax = dens(i) 1776 if (dens(i).lt.pmin) pmin = dens(i) 1777 end do 1778 1779 end do 1780 1781 rng1 = pmin 1782 rng2 = pmax 1783 1784 call messg(22) 1785 1786 ofrst = .false. 1787 1788 lstr = 'found grid file' 1789 call inferr(lstr,0) 1790 return 1791 1792100 istat = 0 1793 lstr = 'error reading GRID file' 1794 call inferr(lstr,0) 1795 return 1796 end 1797 1798 subroutine rdccpd(npts1,npts2,npts3,iun,istat, 1799 & denn,dens,pmnn,ichx,a,b,c) 1800 implicit double precision (a-h,o-z) 1801 common /comsrf/ vo(3), r(3),v1(3),v2(3),v3(3), wo(3), 1802 & sl(3),isl 1803 common /grdhlp/ mx3d,mx3d2 1804 common /cnthlp/ rng1,rng2,vlcnt,vlcnt2 1805 integer currec,bigend 1806 common /rdrec/ iund,currec 1807 common /rdwr/ iun1,iun2,iun3,iun4,iun5 1808 character*137 lstr 1809 dimension denn(*),pmnn(*),dens(*) 1810 1811 integer*4 nfast,nmed,nslow,mode,nsfast,nsmed,nsslow 1812 integer*4 mx,my,mz,mapc,mapr,maps,ispg,nsymbt 1813 real*4 abc(3),alpha,beta,gamma,dmin,dmax,dmean 1814 real*4 ori(3),dtemp,r1,r2,r3 1815 1816 equivalence (buff(1),nfast) 1817 equivalence (buff(5),nmed) 1818 equivalence (buff(9),nslow) 1819 equivalence (buff(13),mode) 1820 equivalence (buff(17),nsfast) 1821 equivalence (buff(21),nsmed) 1822 equivalence (buff(25),nsslow) 1823 equivalence (buff(29),mx) 1824 equivalence (buff(33),my) 1825 equivalence (buff(37),mz) 1826 equivalence (buff(41),abc(1)) 1827 equivalence (buff(53),alpha) 1828 equivalence (buff(57),beta) 1829 equivalence (buff(61),gamma) 1830 equivalence (buff(65),mapc) 1831 equivalence (buff(69),mapr) 1832 equivalence (buff(73),maps) 1833 equivalence (buff(77),dmin) 1834 equivalence (buff(81),dmax) 1835 equivalence (buff(85),dmean) 1836 equivalence (buff(89),ispg) 1837 equivalence (buff(93),nsymbt) 1838 equivalence (buff(197),ori(1)) 1839 character*1 buff(1024) 1840 character*1 extbuf(1024) 1841 dimension n(3),ns(3),dtemp(10000) 1842 1843 1844 todeg = 45.0d0 / datan(1.0d0) 1845 toang = 0.52917706d0 1846 istat = 1 1847 lstr = ' ' 1848 1849 read(iund,pos=1) buff 1850 1851 a = dble(abc(1)) 1852 b = dble(abc(2)) 1853 c = dble(abc(3)) 1854 1855 write(iun3,*) " " 1856 write(iun3,'(3(a,f7.3))') ,"a=",abc(1)," b=",abc(2)," c= ",abc(3) 1857 write(iun3,*) " " 1858 1859 abc(1) = abc(1) / toang 1860 abc(2) = abc(2) / toang 1861 abc(3) = abc(3) / toang 1862 1863 write(iun3,*) " " 1864 write(iun3,'(3(a,f7.3))') ,"alpha ",alpha," beta ",beta, 1865 & " gamma ",gamma 1866 write(iun3,*) " " 1867 1868 alpha = alpha / todeg 1869 beta = beta / todeg 1870 gamma = gamma / todeg 1871 1872 abc(1) = abc(1) / float(mx) 1873 abc(2) = abc(2) / float(my) 1874 abc(3) = abc(3) / float(mz) 1875 1876 v1(1) = dble(abc(1)) 1877 v1(2) = 0.0d0 1878 v1(3) = 0.0d0 1879 1880 v2(1) = dble(cos(gamma) * abc(2)) 1881 v2(2) = dble(sin(gamma) * abc(2)) 1882 v2(3) = 0.0d0 1883 1884 z1 = dble(cos(beta)) 1885 z2 = dble((cos(alpha) - cos(beta)*cos(gamma)) / sin(gamma)) 1886 z3 = dble(sqrt(1.0 - z1*z1 - z2*z2)); 1887 1888 v3(1) = z1 * dble(abc(3)) 1889 v3(2) = z2 * dble(abc(3)) 1890 v3(3) = z3 * dble(abc(3)) 1891 1892c Convert the origin from grid space to cartesian coordinates 1893c nsfast,nsmed,nsslow convert to vo 1894 1895 n(1) = nfast 1896 n(2) = nmed 1897 n(3) = nslow 1898 1899 ns(1) = nsfast 1900 ns(2) = nsmed 1901 ns(3) = nsslow 1902 1903c print*,"ns(1-3) ",ns(1),ns(2),ns(3) 1904 1905 vo(1) = v1(1) * dble(ns(mapc)) + 1906 & v2(1) * dble(ns(mapr)) + v3(1) * dble(ns(maps)) 1907 vo(2) = v2(2) * dble(ns(mapr)) + v3(2) * dble(ns(maps)) 1908 vo(3) = v3(3) * dble(ns(maps)) 1909 1910c normalise basis vectors 1911 1912 call vsc1(v1,1.0d0,1.0d-4) 1913 call vsc1(v2,1.0d0,1.0d-4) 1914 call vsc1(v3,1.0d0,1.0d-4) 1915 1916 isl = 1 1917 sl(1) = 1.0d0 1918 sl(2) = b/a 1919 sl(3) = c/a 1920 call vnrm(sl) 1921 1922 1923 r(1) = abc(mapc)*dble(n(mapc)) 1924 r(2) = abc(mapr)*dble(n(mapr)) 1925 r(3) = abc(maps)*dble(n(maps)) 1926 1927 r1 = real(r(mapc)) 1928 r2 = real(r(mapr)) 1929 r3 = real(r(maps)) 1930 1931c nfast,nmed,nslow corresponds to npts1,npts2,npts3 1932c nsfast,nsmed,nsslow corresponds to 1933 1934c print*," nfast=",nfast," nmed=",nmed," nz=",nslow 1935c print*,"mode=",mode 1936c print*,"nsfast=",nsfast," nsmed=",nsmed, 1937c + " nsslow=",nsslow 1938c print*," mx=",mx," my=",my," mz=",mz 1939c print*,"a=",a," b=",b," c=",c 1940c print*,"alpha=",alpha," beta=",beta," gamma=",gamma 1941c print*,"mapc=",mapc," mapr=",mapr," maps=",maps 1942 write(iun3,'(3(a,f7.3))') "dmin= ",dmin," dmax= ",dmax, 1943 & " dmean= ",dmean 1944 write(iun3,*) " " 1945c print*,"ispg=",ispg 1946c print*,"nsymbt=",nsymbt 1947c print*,"ORI=",(ori(i),i=1,3) 1948 1949c inquire(iund,pos=ipos) 1950 read(iund,pos=1024)(extbuf(i),i=1,nsymbt) 1951 1952c print*,extbuf(1:nsymbt) 1953 1954 pmax = -1000000.d0 1955 pmin = 1000000.d0 1956 1957c print*,"npts1,npts2,npts3 ",n(mapc),n(mapr),n(maps) 1958 npts1 = n(mapc) 1959 npts2 = n(mapr) 1960 npts3 = n(maps) 1961 mxind = mx3d*mx3d2 1962 1963 no = n(mapc)*n(mapr) 1964 1965 write(iun3,'(a,3(x,i3))') "slow,medium,fast ", 1966 & n(maps),n(mapc),n(mapr) 1967 write(iun3,*) " " 1968 1969 do k=1,n(maps) 1970 ipos = 1024+nsymbt+(k-1)*no*4 + 1 1971 read(iund,pos=ipos) (dtemp(ii),ii=1,no) 1972 1973 do i=1,n(mapc) 1974 do j=1,n(mapr) 1975 ind = i + (j-1)*n(mapc) + (k-1)*mx3d2 1976 1977 if (ind.le.mxind) then 1978 ii = j + n(mapr)*(i-1) 1979 denn(ind) = dtemp(ii) 1980 1981 if (dtemp(ii).gt.pmax) pmax = dtemp(ii) 1982 if (dtemp(ii).lt.pmin) pmin = dtemp(ii) 1983 1984 else 1985 print*,"ind ",ind," mx3d2 "," error d=",dtemp(ii) 1986 endif 1987 end do 1988 end do 1989 end do 1990 1991 close(unit=iund) 1992 1993 vlcnt = 1.00d0 1994 lstr = 'found ccp4 file' 1995 1996 call inferr(lstr,0) 1997 1998 rng1 = pmin 1999 rng2 = pmax 2000 2001 ichx = 1 2002 2003 call messg(16) 2004 2005 return 2006 2007100 istat = 0 2008 lstr = 'error reading ccp4 file' 2009 call inferr(lstr,0) 2010 return 2011 2012 end 2013 2014 subroutine rdomad(npts1,npts2,npts3,iun,istat, 2015 & denn,dens,pmnn,ichx) 2016 implicit double precision (a-h,o-z) 2017 parameter (mxdir=500) 2018 logical dolap, lapdbl 2019 common /vropt/ ivtwo,ihand,ivadd 2020 common /comsrf/ vo(3), r(3),v1(3),v2(3),v3(3), wo(3), 2021 & sl(3),isl 2022 common /grdhlp/ mx3d,mx3d2 2023 common /cnthlp/ rng1,rng2,vlcnt,vlcnt2 2024 character*137 lstr 2025 integer*2 buff(256),i1to2 2026 integer*1 bff(512) 2027 integer currec,bigend 2028 common /rdrec/ iund,currec 2029 common /isend/ bigend,nlen 2030 equivalence (buff,bff) 2031 dimension denn(*),pmnn(*),dens(*), fdum(3) 2032 2033c nnx,nny,nnz the number of points in x, y and z direction 2034c ra real step along x,y and z direction per point 2035c rx,ry,rz origin of the grid 2036 2037 istat = 1 2038 lstr = ' ' 2039 2040 toang = 0.52917706d0 2041 iplat = 0 2042 2043 currec = -1 2044 call getrc2(buff,ierr) 2045 2046c 7-9 number of grid points in x,y,z 2047 2048 npts1 = int(buff(7)) 2049 npts2 = int(buff(8)) 2050 npts3 = int(buff(9)) 2051 2052 print*,' ' 2053 print*,'O/DSN6 map file:' 2054 print*,' ' 2055 print*,'npts1 npts2 npts3 ',npts1,npts2,npts3 2056 2057c 1-3 lower limits for x,y,z 2058 2059 ior1 = int(buff(1)) 2060 ior2 = int(buff(2)) 2061 ior3 = int(buff(3)) 2062 2063 print*,'orgx,orgy,orgz ',ior1,ior2,ior3 2064 2065c 4-6 limit range for x,y,z 2066 2067 iex1 = int(buff(4)) 2068 iex2 = int(buff(5)) 2069 iex3 = int(buff(6)) 2070 2071 print*,'iex1,iex2,iex3 ',iex1,iex2,iex3 2072 2073c the grid does not cover the whole unit cell, but just the molecule 2074c the grid therefor does not run from 1 to npts in each direction but 2075c from ior1 to ior1 + iex1 -1 (in the x direction) 2076c we need npts1,npts2,npts3 however to calculate the spacing between 2077c two consequetive grid points in all three directions 2078c (together with a,b,c) 2079c 2080c the problem with existing molden maps is they have equally many 2081c points to the right and left of the origin (vo) 2082c in these new maps the origin is not the center of the map, but 2083c the start of the map in all three directions (npts1=npts2=npts3=1) 2084 2085 2086c 18 F1 2087 2088 sc = 1.0d0 / dble(buff(18)) 2089 sca = 1.0d0 / (dble(buff(18)) * toang) 2090 2091c the spacing between two consequetive grid points in all three directions 2092 2093c 10-12 cell dimensions (times F1) in x,y,z 2094c sca gets F1 (buff(18) out again 2095 2096 a = sca * dble(buff(10)) / dble(npts1) 2097 b = sca * dble(buff(11)) / dble(npts2) 2098 c = sca * dble(buff(12)) / dble(npts3) 2099 2100 print*,'a b c',sc*dble(buff(10)),sc*dble(buff(11)), 2101 * sc*dble(buff(12)) 2102 2103c 13-15 cell angles (times F1) in x,y,z 2104c sc gets F1 (buff(18) out again 2105 2106 alpha = sc * dble(buff(13)) * 3.14159265d0 / 180.0d0 2107 beta = sc * dble(buff(14)) * 3.14159265d0 / 180.0d0 2108 gamma = sc * dble(buff(15)) * 3.14159265d0 / 180.0d0 2109 2110c 16 multiplicative term (times F2) for density value to map to [0,255] 2111c 17 additive term (times F2) for density value to map to [0,255] 2112c 19 F2 2113 2114 prod = dble(buff(16)) / dble(buff(19)) 2115 plus = dble(buff(17)) / dble(buff(19)) 2116 2117 print*,'alpha beta gamma ',sc*dble(buff(13)),sc*dble(buff(14)), 2118 & sc*dble(buff(15)) 2119 2120 v1(1) = a 2121 v1(2) = 0.0d0 2122 v1(3) = 0.0d0 2123 2124 v2(1) = dcos(gamma) * b 2125 v2(2) = dsin(gamma) * b 2126 v2(3) = 0.0d0 2127 2128 z1 = dcos(beta) 2129 z2 = (dcos(alpha) - dcos(beta)*dcos(gamma)) / dsin(gamma) 2130 z3 = dsqrt(1.0 - z1*z1 - z2*z2); 2131 2132 v3(1) = z1 * c 2133 v3(2) = z2 * c 2134 v3(3) = z3 * c 2135 2136c Convert the origin from grid space to cartesian coordinates 2137 2138 vo(1) = v1(1) * dble(ior1) + 2139 & v2(1) * dble(ior2) + v3(1) * dble(ior3) 2140 vo(2) = v2(2) * dble(ior2) + v3(2) * dble(ior3) 2141 vo(3) = v3(3) * dble(ior3) 2142 2143c normalise basis vectors 2144 2145 call vsc1(v1,1.0d0,1.0d-4) 2146 call vsc1(v2,1.0d0,1.0d-4) 2147 call vsc1(v3,1.0d0,1.0d-4) 2148 2149 isl = 1 2150 sl(1) = 1.0d0 2151 sl(2) = b/a 2152 sl(3) = c/a 2153 call vnrm(sl) 2154 2155 r(1) = a*dble(iex1-1) 2156 r(2) = b*dble(iex2-1) 2157 r(3) = c*dble(iex3-1) 2158 2159 print*,'r ',(r(i),i=1,3) 2160 npts1 = iex1 2161 npts2 = iex2 2162 npts3 = iex3 2163 2164c spacing between consequetive points in the three directions 2165c rf1 = r(1) / dble(iex1-1) 2166c rf2 = r(2) / dble(iex2-1) 2167c rf3 = r(3) / dble(iex3-1) 2168c 2169c do kc=1,npts3 do ic=1,npts2 do jc=1,npts3 2170c coordinates of each grid point (jc,ic,kc) rt(1,2,3) 2171c rt(1) = v1(1)*(jc-1)*rf1 + v2(1)*(ic-1)*rf2 + v3(1)*(kc-1)*rf3 + vo(1) 2172c rt(2) = v1(2)*(jc-1)*rf1 + v2(2)*(ic-1)*rf2 + v3(2)*(kc-1)*rf3 + vo(2) 2173c rt(3) = v1(3)*(jc-1)*rf1 + v2(3)*(ic-1)*rf2 + v3(3)*(kc-1)*rf3 + vo(3) 2174c 2175c we have to change rotbck to: 2176c do i=1,3 2177c wn(i) = v1(i)*(w1)*r(1) + 2178c & v2(i)*(w2)*r(2) + 2179c & v3(i)*(w3)*r(3) + vo(i) 2180c end do 2181c 2182c alternatively we could make it a general purpose routine by specifying 2183c in common comsrf what is the location of the origin in grid coordinates 2184c 2185 2186 pmax = -1000000.d0 2187 pmin = 1000000.d0 2188 2189 ibrik = npts1/8 2190 if (npts1.gt.ibrik*8) ibrik = ibrik + 1 2191 jbrik = npts2/8 2192 if (npts2.gt.jbrik*8) jbrik = jbrik + 1 2193 kbrik = npts3/8 2194 if (npts3.gt.kbrik*8) kbrik = kbrik + 1 2195 2196 mxind = mx3d*mx3d2 2197 2198 do k=1,kbrik 2199 do j=1,jbrik 2200 do i=1,ibrik 2201 2202 call getrc2(buff,ierr) 2203 2204 do iz=1,8 2205 if ((iz + (k-1)*8).le.npts3) then 2206 do iy=1,8 2207 if ((iy + (j-1)*8).le.npts2) then 2208 do ix=1,8 2209 2210 if ((ix + (i-1)*8).le.npts1) then 2211 2212 ipnt = ix + 2213 & (iy-1)*8 + 2214 & (iz-1)*64 2215 den = dble(i1to2(bff(ipnt))) - plus 2216 den = den / prod 2217 2218 iind = (ix+(i-1)*8) + 2219 & ((iy+(j-1)*8)-1)*npts1 + 2220 & ((iz+(k-1)*8)-1)*mx3d2 2221 2222 if (iind.le.mxind) then 2223 denn(iind) = den 2224 2225 if (den.gt.pmax) pmax = den 2226 if (den.lt.pmin) pmin = den 2227 endif 2228 2229 endif 2230 2231 end do 2232 endif 2233 end do 2234 endif 2235 end do 2236 2237 end do 2238 end do 2239 end do 2240 2241 do i=1,npts3 2242 pmnn(i) = 100000.d0 2243 end do 2244 2245 print*,'Dens Max=',pmax,' Dens Min=',pmin 2246 rng1 = pmin 2247 rng2 = pmax 2248 2249 vlcnt = 1.00d0 2250 lstr = 'found omap file' 2251 call inferr(lstr,0) 2252 2253 ichx = 1 2254 2255 close(iun) 2256 call messg(16) 2257 2258 return 2259 2260100 istat = 0 2261 lstr = 'error reading OMAP/DSN6 file' 2262 call inferr(lstr,0) 2263 return 2264 end 2265 2266 subroutine dpomad(iopt,denn) 2267 implicit double precision (a-h,o-z) 2268 logical dolap, lapdbl 2269 common /comsrf/ vo(3), r(3),v1(3),v2(3),v3(3), wo(3), 2270 & sl(3),isl 2271 common /cnthlp/ rng1,rng2,vlcnt,vlcnt2 2272 integer srfmap, srfloc 2273 common /hlpsrf/ npts1,npts2,npts3,srfmap,srfloc,ifogl,itsrf 2274 common /vropt/ ivtwo,ihand,ivadd 2275 dimension origin(3),rt(3) 2276 dimension denn(*),fdum(3) 2277 2278 call curs(1) 2279 wo(1) = 0.0d0 2280 wo(2) = 0.0d0 2281 wo(3) = 0.0d0 2282 2283 do i=1,3 2284 origin(i) = 0.0d0 2285 rt(i) = 1.0d0 2286 end do 2287 2288 ipsi = 0 2289 dolap = .false. 2290 lapdbl = .false. 2291 if (iopt.eq.1) then 2292 itrans = 0 2293 else 2294 itrans = 1 2295 endif 2296 mapit = 0 2297 2298 ivtwot = ivtwo 2299 ivtwo = 4 2300 2301 call mcubes(npts1,npts2,npts3,denn,fdum,mapit,origin,vlcnt, 2302 & vlcnt2,rt,ipsi,dolap,lapdbl,itrans,iun) 2303 2304 ivtwo = ivtwot 2305 2306 call curs(0) 2307 2308 return 2309 end 2310 2311 subroutine byter(intin,intout) 2312 implicit none 2313 integer intin, intout,bigend,nlen 2314 common /isend/ bigend,nlen 2315 2316 if (bigend.eq.1) then 2317 intout = intin 2318 else 2319 call mvbits( intin, 24, 8, intout, 0 ) 2320 call mvbits( intin, 16, 8, intout, 8 ) 2321 call mvbits( intin, 8, 8, intout, 16 ) 2322 call mvbits( intin, 0, 8, intout, 24 ) 2323 endif 2324 2325 return 2326 end 2327 2328 subroutine getrec(buff,buflen,isil,ierr) 2329 implicit real (a-h,o-z) 2330 character*4 buff(*) 2331 integer recln,currec,ierr,buflen 2332 common /rdrec/ iund,currec 2333 2334 ierr = 0 2335 currec = currec + 1 2336 2337 read(iund,rec=currec,err=100) intin 2338 2339 call byter(intin,intout) 2340 2341 recln = intout/4 2342 2343 if (recln.gt.buflen) goto 100 2344 2345 do i =1,recln 2346 read(iund,rec=currec+i,err=100) buff(i) 2347 end do 2348 2349 currec = currec + recln 2350 currec = currec + 1 2351 2352 read(iund,rec=currec,err=100) intin 2353 2354 2355 call byter(intin,intchk) 2356 2357 if (intout.ne.intchk) then 2358 ierr = 1 2359 if (isil.eq.0) print*,'getrec: error reading file' 2360 endif 2361 2362 return 2363 2364100 ierr = 1 2365 return 2366 end 2367 2368 subroutine bytr2(intin) 2369 implicit none 2370 integer bigend,nlen,i 2371 integer*2 intin(256),intout 2372 common /isend/ bigend,nlen 2373 2374 do i=1,256 2375 if (bigend.ne.1) then 2376 call mvbits( intin(i), 8, 8, intout, 0 ) 2377 call mvbits( intin(i), 0, 8, intout, 8 ) 2378 intin(i) = intout 2379 endif 2380 end do 2381 2382 return 2383 end 2384 2385 integer*2 function i1to2(intin) 2386 integer*1 intin 2387 2388 i1to2 = intin 2389 if (i1to2.lt.0) i1to2 = abs(intin) + 128 2390 return 2391 end 2392 2393 subroutine getrc2(buff,ierr) 2394 implicit real (a-h,o-z) 2395 integer*2 buff(256) 2396 integer currec,ierr 2397 integer*2 into 2398 common /rdrec/ iund,currec 2399 2400 ierr = 0 2401 currec = currec + 1 2402 2403 do i =1,256 2404 read(iund,rec=currec+i,err=100) buff(i) 2405 end do 2406 2407 call bytr2(buff) 2408 if (currec.eq.0) then 2409 if (buff(19).ne.100) goto 100 2410 endif 2411 2412 currec = currec + 255 2413 2414 return 2415 2416100 ierr = 1 2417 return 2418 end 2419 2420 subroutine fndor(idebug) 2421 implicit double precision ( a-h,o-z) 2422 character*137 line 2423 common /gauori/ nzm,nso,nio,nzo,ioropt,ifor, 2424 & ixyz98,iopr,isymm,irc,imp2,icntp,itd 2425 2426 if (idebug.eq.1) print*,'start find orientations' 2427 ifor = 0 2428 nzm = 0 2429 nso = 0 2430 nio = 0 2431 nzo = 0 2432 istat = 1 2433 do while(istat.eq.1) 2434 call searchq(line,'Z-MATRIX (ANGSTROMS', 2435 & 'Standard orientation:', 2436 & 'Input orientation:', 2437 & 'Z-Matrix orientation:',istat) 2438 if (icdex(line,'Z-MATRIX (ANGSTROMS').ne.0) nzm = nzm + 1 2439 if (icdex(line,'Standard orientation:').ne.0) nso = nso + 1 2440 if (icdex(line,'Input orientation:').ne.0) nio = nio + 1 2441 if (icdex(line,'Z-Matrix orientation:').ne.0) nzo = nzo + 1 2442 end do 2443 2444 if (nso.eq.nio) ioropt = 1 2445 if (nso.eq.0.and.nio.eq.0) ioropt = 0 2446 if (nso.gt.nio) ioropt = 2 2447 if (nso.lt.nio) ioropt = 3 2448 2449 if (idebug.eq.1) then 2450 print*,'nzm=',nzm,' nso=',nso,' nio=',nio,' nzo=',nzo 2451 endif 2452 2453 return 2454 end 2455 2456 subroutine wrcubd(npts1,npts2,npts3,ipsi,denn) 2457 2458c THIS IS REALLY wrcube 2459 2460 implicit double precision (a-h,o-z) 2461 parameter (numatm=2000) 2462 common /plane/ px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat 2463 common /moldat/ natoms, norbs, nelecs,nat(numatm) 2464 common /coord / xyz(3,numatm) 2465 common /grdhlp/ mx3d,mx3d2 2466 dimension v3(3),o(3),p(3),dtmp(6),denn(*) 2467 2468 iun = 21 2469 2470c normalize the perpendicular vector 2471 2472 rn = dsqrt(cx*cx+cy*cy+cz*cz) 2473 2474 v3(1) = cx / rn 2475 v3(2) = cy / rn 2476 v3(3) = cz / rn 2477 2478 p(1) = px 2479 p(2) = py 2480 p(3) = pz 2481 2482c calculate the origin for the cube 2483 2484 do i=1,3 2485 o(i) = p(i) - 0.5d0*(v1(i)*r(1) + v2(i)*r(2) + v3(i)*r(3)) 2486 end do 2487 2488 write(iun,'(a)') 'Molden generated cube file' 2489 2490 if (ipsi.eq.0) then 2491 write(iun,*) 'Density' 2492 else 2493 if (ipsi.lt.0) then 2494 write(iun,'(''Beta Orbital '',i4)') iabs(ipsi) 2495 else 2496 write(iun,'(''Orbital '',i4)') ipsi 2497 endif 2498 endif 2499 2500c number of atoms and origin of the cube 2501 2502 write(iun,'(i5,4(f12.6))') natoms,(o(i),i=1,3) 2503 2504c number of points in each dimension, and voxel size 2505 2506 write(iun,'(i5,4(f12.6))') npts1, 2507 & (v1(i)*r(1)/dble(npts1-1),i=1,3) 2508 write(iun,'(i5,4(f12.6))') npts2, 2509 & (v2(i)*r(2)/dble(npts2-1),i=1,3) 2510 write(iun,'(i5,4(f12.6))') npts3, 2511 & (v3(i)*r(3)/dble(npts3-1),i=1,3) 2512 2513c atomic numbers and coordinates of the atoms (in bohrs) 2514 2515 do i=1,natoms 2516 write(iun,'(i5,4(f12.6))') nat(i),0.0,(xyz(j,i),j=1,3) 2517 end do 2518 2519c density data 2520 2521 ij = 0 2522 2523 do i=1,npts1 2524 do j=1,npts2 2525 ij = ij + 1 2526 kk = 0 2527 do while (kk.lt.npts3) 2528 n = 6 2529 if (npts3-kk.lt.n) n = npts3-kk 2530 do k=1,n 2531 dtmp(k) = denn((npts3-(kk+k))*mx3d2 + ij) 2532 if (dtmp(k).lt.0.0d0) then 2533 if (dtmp(k).gt.-1.0d-99) dtmp(k) = -1.0d-99 2534 else 2535 if (dtmp(k).lt.1.0d-99) dtmp(k) = 1.0d-99 2536 endif 2537 end do 2538 write(iun,'(6E13.5)') (dtmp(k),k=1,n) 2539 kk = kk + 6 2540 end do 2541 end do 2542 end do 2543 2544 return 2545 end 2546 2547 subroutine nmrshl 2548 implicit double precision (a-h,o-z) 2549 parameter (numatm=2000) 2550 character*137 line 2551 common /nmr/ shlnuc(numatm),ihsnmr 2552 common /moldat/ natoms, norbs, nelecs,nat(numatm) 2553 common /rdwr/ iun1,iun2,iun3,iun4,iun5 2554 2555 call searchd(line,'GIAO Magnetic shielding', 2556 & 'Magnetic shielding (ppm)',istat) 2557 if (istat.ne.0) then 2558 do i=1,natoms 2559 call redel(line,1) 2560 i1 = icdex(line,'Isotropic =') 2561 if (i1.ne.0) then 2562 read(line(i1+11:i1+22),'(f11.4)') shlnuc(i) 2563 endif 2564 call redel(line,4) 2565 end do 2566 ihsnmr = 1 2567 2568 call search(line,'Total nuclear spin-spin coupling J',istat) 2569 if (istat.ne.0) ihsnmr = 2 2570 else 2571 call rewfil 2572 call searchd(line,'Total nuclear spin-spin coupling J', 2573 & 'Fermi Contact (FC) contribution to J', 2574 & istat) 2575 if (istat.ne.0) ihsnmr = 2 2576 endif 2577 2578 return 2579 end 2580 2581 subroutine nmrcpd(idebug,couplj) 2582 implicit double precision (a-h,o-z) 2583 parameter (numatm=2000) 2584 character*137 line 2585 common /nmr/ shlnuc(numatm),ihsnmr 2586 common /moldat/ natoms, norbs, nelecs,nat(numatm) 2587 common /rdwr/ iun1,iun2,iun3,iun4,iun5 2588 dimension couplj(*) 2589 2590 ntimes = natoms/5 2591 nt = mod(natoms,5) 2592 if (nt.ne.0) ntimes = ntimes + 1 2593 2594 ndone = 0 2595 do l=1,ntimes 2596 call redel(line,1) 2597 nitems = 5 2598 ndo = ndone + nitems 2599 if (ndo.gt.natoms) nitems = nitems - (ndo - natoms) 2600 do i=(l-1)*5+1,natoms 2601 call redel(line,1) 2602 read(line,'(7x,5f14.10)') 2603 & (couplj((i-1)*natoms+(l-1)*5+j),j=1,nitems) 2604 end do 2605 ndone = ndone + nitems 2606 end do 2607 2608 if (idebug.eq.1) then 2609 do i=1,natoms 2610 print*,'atom ',i,' ',(couplj((i-1)*natoms+j),j=1,i) 2611 end do 2612 endif 2613 2614 return 2615 end 2616 2617 subroutine qmtot 2618 implicit double precision (a-h,o-z) 2619 parameter (numatm=2000) 2620 common /moldat/ natoms, norbs, nelecs,nat(numatm) 2621 common /totchg/ itot 2622 2623 cor = 0.0d0 2624 do i=1,natoms 2625 cor = cor + nat(i) 2626 end do 2627 itot = cor - nelecs 2628 2629 return 2630 end 2631 2632