1 subroutine rdexbas(basfile) 2 implicit double precision (a-h,o-z) 3 parameter (maxbas=15) 4 parameter (mxel=100) 5 parameter (mxaorb=40) 6 common /basext/ ibasin(mxel),eborbs(maxbas), 7 & pbas(mxaorb,mxaorb,maxbas) 8 common /basopt/ extbas 9 logical extbas 10 character*75 basfile,line 11 character*8 aname 12 13 extbas = .false. 14 15 open (unit=99,file=basfile,status='old',err=90) 16 num = 1 17 do while(.true.) 18 19 read (99,'(A8)',end=10,err=100) aname 20 read (99,*,err=100) line 21 write(6,*) aname 22 read (99,'(i4,i10)',err=100) icharge,norbs 23 ibasin(icharge) = num 24 eborbs(num) = norbs 25 write(6,*) 'norbs=', eborbs(num),'icharge=', icharge 26 do i=1,norbs 27 read (99,'(6f15.10)',err=100) (pbas(i,j,num),j=1,norbs) 28 end do 29 num = num + 1 30 read(99,*,err=100) line 31 end do 32 close(99) 33 3410 continue 35 if (num.gt.1) extbas = .true. 36 return 37 3890 continue 39 print*,'error opening file' 40 return 41 42100 continue 43 return 44 end 45 46 subroutine newdenmad( 47 & p,q) 48c THIS IS REALLY newdenmat 49 50c returns density matrix in p common /densty/ 51 implicit double precision (a-h,o-z) 52 parameter (numatm=2000) 53 parameter (maxbas=15) 54 parameter (mxel=100) 55 parameter (mxaorb=40) 56 logical valenc,bonds,ovrlap,atomic,doori,orient,dolap,doelf 57 common /rdwr/ iun1,iun2,iun3,iun4,iun5 58 common /moldat/ natoms, norbs, nelecs,nat(numatm) 59 common /option/ ipsi,idebug,valenc,bonds,ovrlap,atomic,doori, 60 & dolap,doelf 61 62 common /basext/ ibasin(mxel),eborbs(maxbas), 63 & pbas(mxaorb,mxaorb,maxbas) 64 common /orbhlp/ mxorb,iuhf,ispd 65 66 character*2 elemnt 67 common /elem/elemnt(mxel) 68 dimension p(*),q(*) 69 70 do i=1,norbs 71 do j=1,norbs 72 q((j-1)*mxorb + i) = 0.0d0 73 end do 74 end do 75 76 ix = 0 77 do 700 k=1,natoms 78 n = eborbs(ibasin(nat(k))) 79 write(6,*) 'n=', n, 'K=',k 80 if (ovrlap) then 81c 82c============ do overlap densities only ======================== 83c 84 do i=1,n 85 do j=1,n 86 p((ix+j-1)*mxorb + (ix+i)) = 0.0d0 87 end do 88 end do 89 else 90c 91c============ molecule minus atoms ============================= 92c 93 orient = .false. 94C do i=1,norien 95C if (ori(i).eq.k) then 96C orient = .true. 97C ibal = i 98C endif 99C end do 100 if (orient.and.doori) then 101 write (6,*) 'Orienting ' 102C 103 else 104c 105c----- for all other atoms substract spherical density ----- 106c 107 do i=1,n 108 do j=1,n 109c write (6,*) 'before:', ix,i,j,p((ix+j-1)*mxorb 110c $ + (ix+i)), pbas((ix+i),(ix+j),ibasin(nat(k))) 111 p((ix+j-1)*mxorb + (ix+i)) = 112 $ p((ix+j-1)*mxorb + (ix+i)) - pbas((ix+i),(ix+j) 113 $ ,ibasin(nat(k))) 114 q((ix+j-1)*mxorb + (ix+i)) = 115 $ pbas((ix+i),(ix+j),ibasin(nat(k))) 116C write (6,*) 'after:' ,p(ix+i,ix+j) 117 end do 118 end do 119 endif 120c 121c============== clear out the overlap region of the density matrix 122c if ATOMIC is supplied 123c 124 if (atomic) then 125 idavje = n 126 do i=ix+1,ix+idavje 127 if (k.ne.0) then 128 do j=1,ix 129 p((j-1)*mxorb + i) = 0.0d0 130 end do 131 endif 132 if (k.ne.natoms) then 133 do j=ix+idavje+1,norbs 134 p((j-1)*mxorb + i) = 0.0d0 135 end do 136 endif 137 end do 138 endif 139 endif 140 ix = ix + n 141 700 continue 142 if (bonds.and.idebug.eq.1) then 143 write(iun3,*)' ' 144 write(iun3,*)'***** Atomic Density Matrix *****' 145 write(iun3,*)' ' 146 call prev(q,norbs,norbs,mxorb) 147 write(iun3,*)' ' 148 endif 149 150c############### end bonds , overlap , atomic ################## 151 if (idebug.eq.1) then 152 write(iun3,'('' DENSITY MATRIX USED BY MAP'')') 153 call prev(p,norbs,norbs,mxorb) 154 endif 155 156 return 157 end 158