program orbcpl c c rotational strengths between orbitals c implicit none integer*4 nmo,nao,noe,nq,m,ndo,na,nb,noa,nz,nd0,nmo2,io,i, 1iv,j,k,ii,jj,ia,ib,n,n2,OMIN,OMAX parameter (nd0=21) c nd0 .. number of elements c 1.. 3 el. dipole, c 4.. 6 gradient, c 7.. 9 magnetic dipole, c 10..12 magnetic "GIAO" c 13..18 quadrupole c 19..21 GIAO real*8 a,b,ds,rs,g,dx,dy,dz,rx,ry,rz,ms real*8,allocatable::ea(:),eb(:),bp(:),bpd(:),bpa(:) integer*4,allocatable::iua(:),iub(:),vind(:),oind(:), 1vspin(:),ospin(:) logical lopen,lex,lalphaonly character*1 ab(2) character*4 s4 data ab/'A','B'/ OMIN=0 OMAX=0 lalphaonly=.false. inquire(file='ORBCPL.OPT',exist=lex) if(lex)then open(9,file='ORBCPL.OPT') 1001 read(9,9000,end=1111,err=1111)s4 9000 format(a4) if(s4.eq.'OMIN')read(9,*)OMIN if(s4.eq.'OMAX')read(9,*)OMAX if(s4.eq.'ALPH')read(9,*)lalphaonly goto 1001 1111 close(9) endif call wsos(nmo,nao,noe,nq,m,ndo,na,nb,noa,nz,ea,eb,0,iua,iub) allocate(ea(nmo),eb(nmo),iua(nmo),iub(nmo)) call wsos(nmo,nao,noe,nq,m,ndo,na,nb,noa,nz,ea,eb,1,iua,iub) write(6,700)nmo,nao,noe,nq,m,ndo,na,nb,noa,nz 700 format(/,' NMO NAO NOE NQ M NDO NA NB NAT NZ',/, 110i5,/) if(nb.eq.0)then na=ndo nb=ndo write(6,*)'closed shell' lopen=.false. else write(6,*)'open shell' lopen=.true. if(lalphaonly)then write(6,*)'beta forgotten' nb=0 do 7033 i=1,nmo 7033 iub(i)=0 endif endif nmo2=nmo**2 allocate(bp(nd0*nmo2),bpd(nd0*nmo2),bpa(nd0*nmo2)) c : call rwm(bp,nmo,nd0) c : if(lopen)call rwmd(bpd,nmo) inquire(file='PAX.MO.SCR.TXT',exist=lex) c : if(lex)call rwma(bpa,nmo) n=nmo n2=n**2 c sort orbitals: write(6,*)'sorting SOs' allocate(vind(nmo*2),vspin(nmo*2),oind(nmo*2),ospin(nmo*2)) io=0 iv=0 if(lopen)then 1 ia=0 do 2 i=1,nmo if(iua(i).eq.1)then ia=i a=ea(i) goto 101 endif 2 enddo 101 ib=0 do 3 i=1,nmo if(iub(i).eq.1)then ib=i b=eb(i) goto 102 endif 3 enddo if(ia.eq.0.and.ib.eq.0)goto 4 102 if(a.lt.b.or.ib.eq.0)then iua(ia)=0 if(ia.le.na)then io=io+1 oind(io)=ia ospin(io)=1 else iv=iv+1 vind(iv)=ia vspin(iv)=1 endif else iub(ib)=0 if(ib.le.nb)then io=io+1 oind(io)=ib ospin(io)=2 else iv=iv+1 vind(iv)=ib vspin(iv)=2 endif endif goto 1 else io=ndo iv=nmo-ndo do 10 i=1,io oind(i)=i 10 ospin(i)=1 do 11 i=1,iv vind(i)=i+ndo 11 vspin(i)=1 endif 4 write(6,*)io, ' occupied:' do 8 i=1,io write(6,6002)oind(i),ab(ospin(i)) 8 if(mod(i,16).eq.0)write(6,*) 6002 format(i4,a1,$) write(6,*) write(6,*)iv, ' virtual:' do 9 i=1,iv write(6,6002)vind(i),ab(vspin(i)) 9 if(mod(i,16).eq.0)write(6,*) write(6,*) write(6,*)OMIN,' ... ',io,' occupied X 1 ... ',OMAX,' virtual' open(9,file='ORBCPL.TXT') open(19,file='DIPOLES.TXT') if(OMIN.eq.0)OMIN=1 if(OMAX.eq.0)OMAX=iv do 7 ii=omin,io do 7 jj=1,OMAX j=oind(ii) k=vind(jj) dx=0.0d0 dy=0.0d0 dz=0.0d0 rx=0.0d0 ry=0.0d0 rz=0.0d0 if(ospin(ii).eq.vspin(jj))then if(vspin(jj).eq.1)then dx=bp(n2*(1-1)+n*(k-1)+j) dy=bp(n2*(2-1)+n*(k-1)+j) dz=bp(n2*(3-1)+n*(k-1)+j) rx=bp(n2*(7-1)+n*(k-1)+j) ry=bp(n2*(8-1)+n*(k-1)+j) rz=bp(n2*(9-1)+n*(k-1)+j) else dx=bpd(n2*(1-1)+n*(k-1)+j) dy=bpd(n2*(2-1)+n*(k-1)+j) dz=bpd(n2*(3-1)+n*(k-1)+j) rx=bpd(n2*(7-1)+n*(k-1)+j) ry=bpd(n2*(8-1)+n*(k-1)+j) rz=bpd(n2*(9-1)+n*(k-1)+j) endif else if(lex)then dx=bpa(n2*(1-1)+n*(k-1)+j) dy=bpa(n2*(2-1)+n*(k-1)+j) dz=bpa(n2*(3-1)+n*(k-1)+j) rx=bpa(n2*(7-1)+n*(k-1)+j) ry=bpa(n2*(8-1)+n*(k-1)+j) rz=bpa(n2*(9-1)+n*(k-1)+j) endif endif ds=dx*dx+dy*dy+dz*dz rs=rx*dx+ry*dy+rz*dz ms=rx*rx+ry*ry+rz*rz ds=ds*2.541765d0**2 rs=rs*2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 ms=ms*(2.0d0*9.2740154d-21*1.0d-18/1.0d-36)**2 g=0.0d0 if(ms+ds.gt.0.0d0)g=rs/(ms+ds) write(9,900)ii,ab(ospin(ii)),j,jj,ab(vspin(jj)),k,ds,ms,ms+ds,rs,g 900 format(i4,1x,a1,2i4,1x,a1,i4,5(' ',g12.5)) 7 write(19,901)ii,jj,dx,dy,dz,rx,ry,rz 901 format(2i4,6f12.6) close(9) close(19) end subroutine report(s) character*(*) s write(6,*)s stop end subroutine wsos(nmo,nao,noe,nq,m,ndo,na,nb,noa,nz,ea,eb,ic, 1iua,iub) c read SOS.OUT c ic=0 .. only dimensions c 1 .. also MO energies implicit none integer*4 nao,noe,nq,m,ndo,na,nb,noa,nz,nmo,i,iua(*),iub(*),ic real*8 ea(*),eb(*) character*80 s80 open(8,file='SOS.OUT',status='old') nmo=0 nao=0 noe=0 nq=0 m=0 ndo=0 na=0 nb=0 noa=0 nz=0 1 read(8,100,end=888,err=888)s80 100 format(a80) if(s80(2:11).eq.'nmo set to')then read(s80(12:len(s80)),*)nmo endif if(s80(2:26).eq.'Number of basis functions')then read(s80(39:54),*)nao endif if(s80(2:20).eq.'Number of electrons')then read(s80(39:54),*)noe endif if(s80(2:7).eq.'Charge')then read(s80(39:54),*)nq endif if(s80(2:13).eq.'Multiplicity')then read(s80(39:54),*)m endif if(s80(2:35).eq.'Number of doubly occupied orbitals')then read(s80(39:54),*)ndo na=ndo endif if(s80(2:37).eq.'Number of singly occupied a-orbitals')then read(s80(39:54),*)na endif if(s80(2:37).eq.'Number of singly occupied b-orbitals')then read(s80(39:54),*)nb endif if(s80(2:26).eq.'Number of atoms')then read(s80(39:54),*)noa endif if(s80(2:26).eq.'Sum of atomic charges')then read(s80(39:54),*)nz if(ic.eq.0)goto 888 endif if(s80(2:6).eq.'Alpha')then write(6,*)nmo,ic read(8,*)(ea(i),i=1,nmo) do 2 i=1,nmo eb(i)=ea(i) iua(i)=1 2 iub(i)=1 endif if(s80(2:5).eq.'Beta')then read(8,*)(eb(i),i=1,nmo) do 3 i=1,nmo 3 iub(i)=1 goto 888 endif goto 1 888 close(8) return end subroutine rwm(b,n,nd0) implicit none integer*4 i,j,k,n,n2,nd00,nd0 real*8 b(*) real*8 ,allocatable::p(:,:) parameter (nd00=21) character*2 ss(nd00) character*15 fn data ss/'PX','PY','PZ','VX','VY','VZ','LX','LY','LZ', 1 'DX','DY','DZ','XX','XY','XZ','YY','YZ','ZZ', 1 'EX','EY','EZ'/ c 1 2 3 4 5 6 7 8 9 c 10 11 12 13 14 15 16 17 18 c 19 20 21 logical lex c D .. r x rj c E .. r x ri if(nd0.ne.nd00)call report('nd0 <> nd00') allocate(p(n,n)) n2=n**2 do 1 i=1,9 if(i.gt.12.and.i.le.18)then fn='Q'//ss(i)//'.MO.SCR.TXT' else fn=ss(i)//'.MO.SCR.TXT' endif inquire(file=fn,exist=lex) if(lex)then open(38,file=fn) call rmtr38(p,n,n,1) close(38) do 3 j=1,n do 3 k=1,n 3 b(n2*(i-1)+n*(k-1)+j)=p(k,j) write(6,8000)i,fn 8000 format(i4,2x,a15,'read') else write(6,*)ss(i)//'.MO.SCR.TXT not present' do 2 j=1,n do 2 k=1,n 2 b(n2*(i-1)+n*(k-1)+j)=0.0d0 stop endif 1 continue return end subroutine rmtr38(A,n0,n,ic) c ic=0 .. reads triangle of a supposedly symmetric matrix c ic=1 .. reads it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(n0,n) read(38,*) read(38,*) N1=1 1 N3=MIN(N1+4,N) read(38,*) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=N3 130 read(38,*)J,(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 read(38,*) return end subroutine rwmd(b,n) implicit none integer*4 i,j,k,n,n2 real*8 b(*) real*8 ,allocatable::p(:,:) character*3 ss(18) data ss/'PBX','PBY','PBZ','VBX','VBY','VBZ','LBX','LBY','LBZ', 1'WBX','WBY','WBZ','BXX','BXY','BXZ','BYY','BYZ','BZZ'/ logical lex allocate(p(n,n)) n2=n**2 do 1 i=1,9 if(i.gt.12)then inquire(file='Q'//ss(i)//'.MO.SCR.TXT',exist=lex) else inquire(file=ss(i)//'.MO.SCR.TXT',exist=lex) endif if(lex)then if(i.gt.12)then open(38,file='Q'//ss(i)//'.MO.SCR.TXT') else open(38,file=ss(i)//'.MO.SCR.TXT') write(6,*)ss(i)//'.MO.SCR.TXT read' endif call rmtr38(p,n,n,1) close(38) do 3 j=1,n do 3 k=1,n 3 b(n2*(i-1)+n*(k-1)+j)=p(k,j) else write(6,*)ss(i)//'.MO.SCR.TXT not present' do 2 j=1,n do 2 k=1,n 2 b(n2*(i-1)+n*(k-1)+j)=0.0d0 endif 1 continue return end subroutine rwma(b,n) implicit none integer*4 i,j,k,n,n2 real*8 b(*) real*8 ,allocatable::p(:,:) character*3 ss(18) data ss/'PAX','PAY','PAZ','VAX','VAY','VAZ','LAX','LAY','LAZ', 1'WAX','WAY','WAZ','AXX','AXY','AXZ','AYY','AYZ','AZZ'/ logical lex allocate(p(n,n)) n2=n**2 do 1 i=1,9 if(i.gt.12)then inquire(file='Q'//ss(i)//'.MO.SCR.TXT',exist=lex) else inquire(file=ss(i)//'.MO.SCR.TXT',exist=lex) endif if(lex)then if(i.gt.12)then open(38,file='Q'//ss(i)//'.MO.SCR.TXT') else open(38,file=ss(i)//'.MO.SCR.TXT') write(6,*)ss(i)//'.MO.SCR.TXT read' endif call rmtr38(p,n,n,1) close(38) do 3 j=1,n do 3 k=1,n 3 b(n2*(i-1)+n*(k-1)+j)=p(k,j) else write(6,*)ss(i)//'.MO.SCR.TXT not present' do 2 j=1,n do 2 k=1,n 2 b(n2*(i-1)+n*(k-1)+j)=0.0d0 endif 1 continue return end