program guvcde c c polarizability from TDDFT by SOS c implicit none character*80 fo,s80,output character*4 key integer nm0,n,nd,ie,i,j,k,l,np,nmo,ix,ia,nat,grid(3),ib,ic,iep, 1oa,ob,idum,io,iy,iz,ig98,nat0,nptem,iw,nmo2,na,nb,nd0,oc,od,ip, 2nt parameter (nm0=1000,nat0=1000,nd0=9) c nm0 .. maximum number of orbitals and states c nat0 ... maximum number of atoms c nd0 .. number of elements 1..3 el. dipole, c 4..6 gradient, 7..9 magnetic dipole real*8 dl(nm0,3),dv(nm0,3),r(nm0,3),e(nm0),ev(nm0),uux,uuy,uuz, 1ucx,ucy,ucz,x(3,nat0),v(nat0),DD,dist,xp,yp,zp,px,py,pz, 1r0v(nm0),r0l(nm0),al(3,3),ax,ay,az,qua,wden,sum,ciji,tq, 1wmin,wmax,dw,ali(3,3),gamma,w,fr,ds1,rs,angle,fi,cijip,ct, 1qv(nm0,3,3),eau(nm0),bohr,db,wtem,wmintem,dwtem,djg(nd0), 1rsr1,rsv1,rsr2,rsv2,un(3),u0(nd0),dsr,dsv,dkg(nd0),djk(nd0),ekj, 1ekg,eps,rsv,rsr,debye,xau(3,nat0),ut(3,nm0),rsr3,rsv3,cback, 1degper,degper0,deglim,ecm,u0d(nd0) integer*4 qz(nat0) c x atomic coordinates c v atomic weights c iz atomic number logical lex,lden,lgeo,lzmat,lw,lrds,lpx,lmcd,lcub,ldeg,lopen real*8, allocatable::cij(:),rox(:),orb(:),roy(:),roz(:),ro(:), 1bp(:),dij(:),bpd(:) c cij - alpha orbitals c dij - beta orbitals integer*4, allocatable::ni(:),aij(:),bij(:),ind(:),jnd(:), 1daij(:),dbij(:),nid(:) c tq=dsqrt(2.0d0) bohr=0.529177d0 debye=2.541765d0 degper0=0.0001d0 lopen=.false. deglim=0.00001d0 nmo=0 write(6,*)' Reads UVCD from Gaussian output' write(6,*)' and writes it in a TAB file' write(6,*)' plus SOS polarizabilities.' write(6,*) lden=.false. lmcd=.true. lrds=.false. lgeo=.false. lw=.false. wden=532.0d0 fo='G.OUT' lzmat=.true. ldeg=.false. DD=2.0d0 open(9,file='GUVCDE.OPT') 1001 read(9,90,err=109,end=109)key 90 format(a4) c gaussian TDDFT output: if(key.eq.'FILE')read(9,9010)fo 9010 format(a80) c use z-matrix orientation: if(key.eq.'ZMAT')read(9,*)lzmat c remove degeneracy if present: if(key.eq.'LDEG')read(9,*)ldeg c frequency to produce .cub file if(key.eq.'WDEN')read(9,*)wden c produce .cub file: if(key.eq.'LDEN')read(9,*)lden c produce atomic weights: if(key.eq.'LWEI')read(9,*)lw c atomix extent for the weigths in A: if(key.eq.'EXTE')read(9,*)DD c calculate MCD: if(key.eq.'LMCD')read(9,*)lmcd c recalculate dipole strengths: if(key.eq.'LRDS')read(9,*)lrds c degeneracy limit detection: if(key.eq.'DEGL')read(9,*)deglim c degeneracy perturbation: if(key.eq.'DEGP')read(9,*)degper0 goto 1001 109 close(9) write(6,6070)fo,lzmat,wden,lden,lw,DD,ldeg 6070 format('FILENAME',/,a80,/,'ZMAT',/,l1,/,'WDEN (cm-1)',/,f10.2, 1/,'LDEN',/,l1,/,'LWEIGHT',/,l1,/,'EXTEND',/,f10.2,/,'LDEG',/,l1) open(2,file=fo) n=0 ie=0 nat=0 c 1 read(2,2000,end=2,err=2)s80 2000 format(a80) IF((lzmat.and.(s80(19:39).EQ.'Z-Matrix orientation:'.OR. 1 s80(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s80(20:37).EQ.'Input orientation:'.OR. 1 s80(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (s80(20:40).EQ.'Standard orientation:'.OR. 2 s80(26:46).EQ.'Standard orientation:')))THEN ig98=0 if(s80(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s80(27:44).EQ.'Input orientation:'.OR. 1 s80(26:46).EQ.'Standard orientation:')ig98=1 write(6,2000)s80 lgeo=.true. DO 2004 I=1,4 2004 READ(2,*) l=0 2005 READ(2,2000)s80 IF(s80(2:4).NE.'---')THEN l=l+1 if(l.gt.nat0)call report('too many atoms') BACKSPACE 2 if(ig98.eq.0)then READ(2,*)qz(l),qz(l),(x(i,l),i=1,3) else READ(2,*)qz(l),qz(l),x(1,l),(x(i,l),i=1,3) endif xau(1,l)=x(1,l)/bohr xau(2,l)=x(2,l)/bohr xau(3,l)=x(3,l)/bohr IF(qz(l).EQ.-1)l=l-1 GOTO 2005 ENDIF nat=l ENDIF if(s80(8:22).eq.'alpha electrons')then read(s80( 1: 6),*)na read(s80(25:31),*)nb write(6,2000)s80 endif if(s80(31:49).eq.'primitive gaussians')then read(s80(1:6),*)nmo write(6,*)nmo,' orbitals' nmo2=nmo**2 endif if(s80(2:44).eq.'Ground to excited state transition electric')then n=0 write(6,2000)s80 read(2,*) 3 read(2,*,end=21,err=21)nd,ax,ay,az n=n+1 if(n.gt.nm0)call report('Too many transitions.') dl(n,1)=ax dl(n,2)=ay dl(n,3)=az goto 3 21 backspace 2 write(6,*)n,' dipole length transitions' endif if(s80(2:51).eq. 1'Ground to excited state transition velocity dipole')then write(6,2000)s80 read(2,*) do 9 i=1,n 9 read(2,*)nd,dv(i,1),dv(i,2),dv(i,3) write(6,*)n,' dipole velocity transitions' endif if(s80(2:55).eq. 1'Ground to excited state transition velocity quadrupole')then write(6,2000)s80 read(2,*) do 7 i=1,n read(2,*)nd,qv(i,1,1),qv(i,2,2),qv(i,3,3),qv(i,1,2),qv(i,1,3), 1 qv(i,2,3) qv(i,3,2)=qv(i,2,3) qv(i,3,1)=qv(i,1,3) 7 qv(i,2,1)=qv(i,1,2) write(6,*)n,' quadrupole velocity transitions' endif c if(s80(2:44).eq.'Ground to excited state transition magnetic')then write(6,2000)s80 read(2,*) do 10 i=1,n 10 read(2,*)nd,(r(i,j),j=1,3) write(6,*)n,' magnetic transitions' endif c if(s80(53:63).eq.'R(velocity)') then c Gaussian 98 does not write R(length) nor R(velocity) do 8 i=1,n read(2,2000)s80 read(s80,*,err=99)nd,r0v(i),r0v(i),r0v(i),r0v(i) goto 991 99 write(6,2000)s80 write(6,*)i stop 991 continue read(2,2000)s80 if(s80(2:6).eq.'Total')then do 81 k=1,4 81 read(2,2000)s80 else backspace 2 endif read(2,2000)s80 if(s80(2:3).eq.'R(')then do 82 k=1,4 82 read(2,2000)s80 else backspace 2 endif 8 continue endif if(s80(54:62).eq.'R(length)') then write(6,2000)s80 do 11 i=1,n 11 read(2,*,end=2,err=2)nd,r0l(i),r0l(i),r0l(i),r0l(i) endif c if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then write(6,2000)s80 read(s80(15:18),*)nd do 203 i=1,79 if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev(nd) 203 if(s80(i:i+1).eq.'nm')read(s80(i- 9:i-2),*)e(nd) c ecm(i)=ev(i)*8065.54476345045d0 eau(nd)=ev(nd)/27.211384205943d0 ie=ie+1 if(ie.ne.nd)call report(' Inconsistent reading') if(lden)then if(nd.eq.1)allocate(cij(nmo2*n),aij(nmo2*n),dij(nmo2*n), 1 bij(nmo2*n),ni(n),daij(nmo2*n),dbij(nmo2*n),nid(n)) ni(ie)=0 nid(ie)=0 111 read(2,2000,end=299,err=299)s80 if(s80(10:11).eq.'->')then if(ni(ie).gt.nmo2)call report('too many coefficients') if(s80(8:8).eq.'A')then c open shell, alpha: lopen=.true. ni(ie)=ni(ie)+1 read(s80( 1: 7),*)aij(nmo2*(ie-1)+ni(ie)) read(s80(12:14),*)bij(nmo2*(ie-1)+ni(ie)) read(s80(18:30),*)cij(nmo2*(ie-1)+ni(ie)) else if(s80(8:8).eq.'B')then c open shell, beta: nid(ie)=nid(ie)+1 lopen=.true. read(s80( 1: 7),*)daij(nmo2*(ie-1)+nid(ie)) read(s80(12:14),*)dbij(nmo2*(ie-1)+nid(ie)) read(s80(18:30),*) dij(nmo2*(ie-1)+nid(ie)) else c closed shell: ni(ie)=ni(ie)+1 read(s80( 1: 8),*) aij(nmo2*(ie-1)+ni(ie)) read(s80(12:15),*) bij(nmo2*(ie-1)+ni(ie)) read(s80(18:30),*) cij(nmo2*(ie-1)+ni(ie)) cij(nmo2*(ie-1)+ni(ie))=cij(nmo2*(ie-1)+ni(ie))*tq endif endif goto 111 endif if(s80(10:11).eq.'<-'.and..not.lopen)then read(s80( 1: 8),*)oa read(s80(12:15),*)ob read(s80(18:30),*)cback cback=cback*tq c just quick fix for the deexitations, closed shell only: do 315 j=1,ni(ie) if(oa.eq.aij(nmo2*(ie-1)+j).and.ob.eq.bij(nmo2*(ie-1)+j))then ct=cij(nmo2*(ie-1)+j) cij(nmo2*(ie-1)+j)=dsqrt(dabs(ct**2-cback**2)) if(ct.lt.0.0d0)cij(nmo2*(ie-1)+j)=-cij(nmo2*(ie-1)+j) goto 111 endif 315 continue goto 111 endif 299 continue write(6,*)ni(ie),' alpha coefficients' write(6,*)nid(ie),' beta coefficients' endif if(nd.eq.n) goto 2 endif goto 1 c 2 close(2) c c if desired, arbitrarily perturb degeneracy if(ldeg)then k=0 do 12 i=1,n-1 degper=degper0 do 12 j=i+1,n if(dabs(eau(i)-eau(j)).lt.deglim)then eau(j)=eau(j)+degper ev(j)=eau(j)*27.211384205943d0 ecm=ev(j)*8065.54476345045d0 e(j)=1.0d7/ecm write(6,6003)i,j,e(i),e(j) 6003 format(2i4,f8.2,' -> ',f8.2) degper=degper+degper0 k=k+1 endif 12 continue write(6,6008)k 6008 format('Degeneracy perturbed ',i3,' times.') do 13 i=1,n ev(i)=eau(i)*27.211384205943d0 ecm=ev(i)*8065.54476345045d0 13 e(i)=1.0d7/ecm endif do 14 i=1,n do 141 j=1,ni(i) if(aij(nmo2*(i-1)+j).gt.nmo.or.aij(nmo2*(i-1)+j).lt.1) 1call report('wrong orbital a') 141 if(bij(nmo2*(i-1)+j).gt.nmo.or.bij(nmo2*(i-1)+j).lt.1) 1call report('wrong orbital b') do 14 j=1,nid(i) if(daij(nmo2*(i-1)+j).gt.nmo.or.daij(nmo2*(i-1)+j).lt.1) 1call report('wrong orbital da') 14 if(dbij(nmo2*(i-1)+j).gt.nmo.or.dbij(nmo2*(i-1)+j).lt.1) 1call report('wrong orbital db') c sign, unit, and Gaussian stupidity correction of dipoles!: do 5 i=1,n do 5 j=1,3 dv(i,j)=-dv(i,j)/eau(i) r( i,j)= r( i,j)/2.0d0 do 5 k=1,3 5 qv(i,j,k)=-qv(i,j,k)/eau(i) c c record spectral tables: c tttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt write(6,*)n,' energies' call titles(fo,'.l.tab',output) write(6,*)output open(31,file=output) call titles(fo,'.v.tab',output) write(6,*)output open(32,file=output) write(31,3000) 3000 format(' UVCD spectrum in length formalism', 1/,' n wavelength (nm) dipole strength (D^2)', 1' rotatory strength (cgs/10**-36)',/,80(1h-)) write(32,3100) 3100 format(' UVCD spectrum in velocity formalism', 1/,' n wavelength (nm) dipole strength (D^2)', 1' rotatory strength (cgs/10**-36)',/,80(1h-)) c do 6 i=1,n c ds=2.12689677295232d-30*f(i)/ev(i)/8065.54476345045d0/1.0d-36 call getv(nm0,dl,ds1,r0l,angle,rs,r,i) write(31,3001)i,e(i),ds1,rs,angle 3001 format(i4,f12.2,2f20.10,f8.1) call getv(nm0,dv,ds1,r0v,angle,rs,r,i) 6 write(32,3001)i,e(i),ds1,rs,angle write(31,3002) write(32,3002) 3002 format(80(1h-)) close(31) close(32) write(6,*) c tttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt c c computation of polarizabilities: gamma=0.01d0 c bandwidth = gamma * w (gamma=0...1) wmin=50.0d0 wmax=500.0d0 np=901 c wavenumber of wavelength: inquire(file='GUVCDE.OPT',exist=lex) if(lex)then open(44,file='GUVCDE.OPT') 201 read(44,440,end=202,err=202)s80 440 format(a80) if(s80(1:5).eq.'GAMMA')read(44,*)gamma if(s80(1:4).eq.'WMIN')read(44,*)wmin if(s80(1:4).eq.'WMAX')read(44,*)wmax if(s80(1:2).eq.'NP')read(44,*)np goto 201 202 close(44) endif dw=(wmax-wmin)/dble(np-1) open(44,file='SOS.TTT') open(9,file='spa.prn') write(44,446)np,wmin,wmax 446 format(i6,2f12.2) write(44,443) 443 format(' Dynamic polarizabilities') write(44,444) 444 format(10x,' Alpha ( f ) Alpha ( g) (pure velocity formalism)', 1/,'l(nm) ux uy uz') w=wmin-dw do 102 k=1,np w=w+dw do 101 i=1,3 do 101 j=1,3 al(i,j)=0.0d0 ali(i,j)=0.0d0 do 101 ie=1,n call setv(eau(ie),w,fr,fi,gamma,0) c a_ab = 2 sum wjn (1/(wjn**2-w**2)) Re ua ub al(i,j) =al(i,j) +dv(ie,i)*dv(ie,j)*fr 101 ali(i,j)=ali(i,j)+dv(ie,i)*dv(ie,j)*fi write(9,69)w,ali(1,1)+ali(2,2)+ali(3,3),al(1,1)+al(2,2)+al(3,3) 69 format(f7.1,2f10.2) write(44,441)w,(al(1,j),j=1,3),(ali(1,j),j=1,3) write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 102 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) 441 format(f8.2,6f10.2) 442 format(8x ,6f10.2) close(9) write(6,*)' spa.prn written' write(44,445) 445 format(10x,' Gp/w ( f) Gp/w ( g )(velocity formalism)', 1/,'l(nm) ux uy uz') open(9,file='spg.prn') w=wmin-dw do 103 k=1,np w=w+dw do 104 i=1,3 do 104 j=1,3 al( i,j)=0.0d0 ali( i,j)=0.0d0 do 104 ie=1,n call setv(eau(ie),w,fr,fi,gamma,1) al( i,j) =al(i,j) +dv(ie,i)*r(ie,j)*fr 104 ali(i,j)=ali(i,j) +dv(ie,i)*r(ie,j)*fi write(9,69)w,-ali(1,1)-ali(2,2)-ali(3,3),al(1,1)+al(2,2)+al(3,3) write(44,441)w,(al(1,j),j=1,3),(ali(1,j),j=1,3) write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 103 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) close(9) write(6,*)' spg.prn written' write(44,447) 447 format(10x,' A ( f) A ( g) (velocity formalism)', 1/,'l(nm) ux uy uz') w=wmin-dw do 105 k=1,np w=w+dw do 105 l=1,3 do 106 i=1,3 do 106 j=1,3 al( i,j)=0.0d0 ali(i,j)=0.0d0 do 106 ie=1,n if(i.eq.j)then qua=(3.0d0*qv(ie,i,j)-qv(ie,1,1)-qv(ie,2,2)-qv(ie,3,3))/2.0d0 else qua=1.5d0*qv(ie,i,j) endif call setv(eau(ie),w,fr,fi,gamma,0) al(i,j) =al( i,j) +dv(ie,l)*qua*fr 106 ali(i,j)=ali(i,j) +dv(ie,l)*qua*fi if(l.eq.1)then write(44,441)w,(al(1,j),j=1,3),(ali(1,j),j=1,3) else write(44,442) (al(1,j),j=1,3),(ali(1,j),j=1,3) endif write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 105 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) close(44) write(6,*)' SOS.TTT written' c c recalculate dipole strengths c dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd if(lrds)then if(.not.lden)call report('lden needed for lrds') inquire(file='PX.MO.SCR.TXT',exist=lpx) if(.not.lpx)call report('*SCR.TXT files needed for lrds') allocate(bp(9*nmo2),bpd(9*nmo2)) call rwm(bp,nmo) if(lopen)call rwmd(bpd,nmo) write(6,*) write(6,*)'Transition moment check:' do 301 ie=1,n do 303 j=1,9 303 djg(j)=0.0d0 do 302 i=1,ni(ie) oa=aij(nmo2*(ie-1)+i) ob=bij(nmo2*(ie-1)+i) ciji=cij(nmo2*(ie-1)+i) do 302 j=1,9 302 djg(j)=djg(j)-ciji*bp(ob+nmo*(oa-1)+nmo2*(j-1)) do 312 i=1,nid(ie) oa=daij(nmo2*(ie-1)+i) ob=dbij(nmo2*(ie-1)+i) ciji=dij(nmo2*(ie-1)+i) do 312 j=1,9 312 djg(j)=djg(j)-ciji*bpd(ob+nmo*(oa-1)+nmo2*(j-1)) if(.not.lopen)then do 313 j=1,9 313 djg(j)=djg(j)*tq endif djg(4)=-djg(4)/eau(ie) djg(5)=-djg(5)/eau(ie) djg(6)=-djg(6)/eau(ie) write(6,3006)e(ie),(djg(j),j=1,9) 301 write(6,3006)e(ie),(dl(ie,j),j=1,3),(dv(ie,j),j=1,3) 3006 format(f6.1,9f8.4) if(lmcd)then c nuclear dipole moment: do 311 j=1,3 un(j)=0.0d0 do 311 l=1,nat 311 un(j)=un(j)+xau(j,l)*qz(l) c c electric dipole, grad and r x grad/2: do 310 j=1,9 u0(j)=0.0d0 u0d(j)=0.0d0 if(lopen)then do 327 i=1,na 327 u0(j)=u0(j)+bp(i+nmo*(i-1)+nmo2*(j-1)) do 329 i=1,nb 329 u0d(j)=u0d(j)+bpd(i+nmo*(i-1)+nmo2*(j-1)) else do 330 i=1,na 330 u0(j)=u0(j)+bp(i+nmo*(i-1)+nmo2*(j-1))*2.0d0 endif 310 continue write(6,3101)(un(j),j=1,3),(-u0(j)-u0d(j),j=1,3), 1 (-u0(j)-u0d(j)+un(j),j=1,3) 1 ,((-u0(j)-u0d(j)+un(j))*debye,j=1,3) 3101 format(/,' Dipole nuclear :',/,3f10.3,/, 1 ' electronic:',/,3f10.3,/, 1 ' total :',/,3f10.3,' debyes:',3f10.3,/) open(31,file='MCDV.TAB') open(32,file='MCDR.TAB') write(31,3007) 3007 format(' MCD spectrum in velocity formalism', 1 /,' n wavelength (nm) dipole strength (D^2)', 1 ' rotatory strength (cgs/10**-36)',/,80(1h-)) write(32,3108) 3108 format(' MCD spectrum in length formalism', 1 /,' n wavelength (nm) dipole strength (D^2)', 1 ' rotatory strength (cgs/10**-36)',/,80(1h-)) do 304 ie=1,n dsr=0.0d0 dsv=0.0d0 rsv1=0.0d0 rsv2=0.0d0 rsr1=0.0d0 rsr2=0.0d0 call getdipoles(ut,nm0,n,cij,nmo,nmo2,ni,bp,aij,bij,u0, 1dij,nid,bpd,daij,dbij,u0d,lopen) c djg=: " singlet (dsqrt(2)) do 321 j=1,9 djg(j)=0.0d0 do 305 i=1,ni(ie) oa=aij(nmo2*(ie-1)+i) ob=bij(nmo2*(ie-1)+i) ciji=cij(nmo2*(ie-1)+i) 305 djg(j)=djg(j)+ciji*bp(ob+nmo*(oa-1)+nmo2*(j-1)) do 322 i=1,nid(ie) oa=daij(nmo2*(ie-1)+i) ob=dbij(nmo2*(ie-1)+i) ciji=dij(nmo2*(ie-1)+i) 322 djg(j)=djg(j)+ciji*bpd(ob+nmo*(oa-1)+nmo2*(j-1)) 321 if(.not.lopen)djg(j)=djg(j)*tq do 320 ia=4,6 320 djg(ia)=-djg(ia)/eau(ie) do 307 iep=1,n if(iep.ne.ie)then ekg=eau(iep) ekj=eau(iep)-eau(ie) c dkg=: "" singlet (dsqrt(2)) do 331 j=1,9 dkg(j)=0.0d0 do 323 i=1,ni(iep) oa=aij(nmo2*(iep-1)+i) ob=bij(nmo2*(iep-1)+i) ciji=cij(nmo2*(iep-1)+i) 323 dkg(j)=dkg(j)+ciji*bp(ob+nmo*(oa-1)+nmo2*(j-1)) do 324 i=1,nid(iep) oa=daij(nmo2*(iep-1)+i) ob=dbij(nmo2*(iep-1)+i) ciji=dij(nmo2*(iep-1)+i) 324 dkg(j)=dkg(j)+ciji*bpd(ob+nmo*(oa-1)+nmo2*(j-1)) 331 if(.not.lopen)dkg(j)=dkg(j)*tq c djk=:"" do 309 j=1,9 309 djk(j)=0.0d0 do 308 ip=1,ni(iep) oc=aij(nmo2*(iep-1)+ip) od=bij(nmo2*(iep-1)+ip) cijip=cij(nmo2*(iep-1)+ip) do 308 i=1,ni(ie) oa=aij(nmo2*(ie-1)+i) ob=bij(nmo2*(ie-1)+i) ciji=cij(nmo2*(ie-1)+i)*cijip if(oa.ne.oc)then if(ob.eq.od)then if(.not.lopen)then do 3071 j=1,9 3071 djk(j)=djk(j)+ciji*bp(oa+nmo*(oc-1)+nmo2*(j-1))*tq else do 3074 j=1,9 3074 djk(j)=djk(j)+ciji*bp(oa+nmo*(oc-1)+nmo2*(j-1)) endif endif else if(ob.ne.od)then if(.not.lopen)then do 3072 j=1,9 3072 djk(j)=djk(j)+ciji*bp(ob+nmo*(od-1)+nmo2*(j-1))*tq else do 3075 j=1,9 3075 djk(j)=djk(j)+ciji*bp(ob+nmo*(od-1)+nmo2*(j-1)) endif else do 3073 j=1,3 3073 djk(j)=djk(j)+ciji*(u0(j)-bp(oa+nmo*(oa-1)+nmo2*(j-1)) 1 + bp(ob+nmo*(ob-1)+nmo2*(j-1))-un(j)) endif endif 308 continue do 325 ip=1,nid(iep) oc=daij(nmo2*(iep-1)+ip) od=dbij(nmo2*(iep-1)+ip) cijip=dij(nmo2*(iep-1)+ip) do 325 i=1,nid(ie) oa=daij(nmo2*(ie-1)+i) ob=dbij(nmo2*(ie-1)+i) ciji=dij(nmo2*(ie-1)+i)*cijip if(oa.ne.oc)then if(ob.eq.od)then do 3076 j=1,9 3076 djk(j)=djk(j)+ciji*bpd(oa+nmo*(oc-1)+nmo2*(j-1)) endif else if(ob.ne.od)then do 3077 j=1,9 3077 djk(j)=djk(j)+ciji*bpd(ob+nmo*(od-1)+nmo2*(j-1)) else do 3078 j=1,3 3078 djk(j)=djk(j)+ciji*(u0d(j)-bpd(oa+nmo*(oa-1)+nmo2*(j-1)) 1 +bpd(ob+nmo*(ob-1)+nmo2*(j-1))-un(j)) endif endif 325 continue c transform gradient to dipole: do 319 ia=1,3 djk(3+ia)=-djk(3+ia)/(-ekj) 319 dkg(3+ia)=-dkg(3+ia)/ekg do 318 ia=1,3 do 318 ib=1,3 do 318 ic=1,3 rsr1=rsr1+eps(ia,ib,ic)*dkg(6+ia)/ekg*djg(ib )*djk(ic ) rsv1=rsv1+eps(ia,ib,ic)*dkg(6+ia)/ekg*djg(ib+3)*djk(ic+3) rsr2=rsr2+eps(ia,ib,ic)*djk(6+ia)/ekj*djg(ib )*dkg(ic ) 318 rsv2=rsv2+eps(ia,ib,ic)*djk(6+ia)/ekj*djg(ib+3)*dkg(ic+3) endif 307 continue rsr3=0.0d0 rsv3=0.0d0 do 316 ia=1,3 do 316 ib=1,3 do 316 ic=1,3 rsv3=rsv3+eps(ia,ib,ic)*djg(6+ia)/eau(ie)*djg(ib+3) 1 *(ut(ic,ie)-u0(ic)-u0d(ic)) 316 rsr3=rsr3+eps(ia,ib,ic)*djg(6+ia)/eau(ie)*djg(ib ) 1 *(ut(ic,ie)-u0(ic)-u0d(ic)) dsr= djg(1)**2+djg(2)**2+djg(3)**2 dsv=(djg(4)**2+djg(5)**2+djg(6)**2)/eau(ie)**2 rsv=-(rsv1+rsv2+rsv3)/2.0d0 rsr=-(rsr1+rsr2+rsr3)/2.0d0 write(6,3003)ie,e(ie),dsv,rsv,rsv1,rsv2,rsv3,rsr,rsr1,rsr2,rsr3 3003 format(i3,10f9.2) write(31,3008)ie,e(ie),dsv,rsv 304 write(32,3008)ie,e(ie),dsr,rsr 3008 format(i4,f12.4,2f20.10,f8.1) write(31,3002) write(32,3002) close(31) close(32) write(6,3005) 3005 format(/,' State dipole sum cij**2:') write(6,3004)0,(un(ia)-u0(ia)-u0d(ia),ia=1,3),1.0d0 do 314 ie=1,n sum=0.0d0 do 317 i=1,ni(ie) 317 sum=sum+cij(nmo2*(ie-1)+i)**2 do 326 i=1,ni(ie) 326 sum=sum+dij(nmo2*(ie-1)+i)**2 314 write(6,3004)ie,(un(ia)-ut(ia,ie),ia=1,3),sum 3004 format(i6,4f10.3) endif deallocate(bp,bpd) endif c dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd c polarization density c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp if(lden)then if(lopen)call report('lden only for closed shells') if(wden.lt.0)then wmintem=wmin dwtem=dw nptem=np else wmintem=wden dwtem=0.0d0 nptem=1 endif wtem=wmintem-dwtem do 909 iw=1,nptem wtem=wtem+dwtem inquire(file='orbitals.cub',exist=lcub) if(.not.lcub)call report('orbitals.cub not present') open(8,file='orbitals.cub') open(9,file='ro_alpha.cub') read(8,2000)s80 write(9,2000)s80 read(8,*) write(9,901) 901 format(' SCF Total Density') read(8,*)nat,ax,ay,az nat=iabs(nat) write(9,2001)nat,ax,ay,az 2001 format(i5,3f12.6) read(8,*)grid(1),px write(9,2001)grid(1),px,0.0d0,0.0d0 read(8,*)grid(2),py,py write(9,2001)grid(2),0.0d0,py,0.0d0 read(8,*)grid(3),pz,pz,pz write(9,2001)grid(3),0.0d0,0.0d0,pz do 902 ia=1,nat read(8,2000)s80 read(s80(18:80),*)(x(i,ia),i=1,3) 902 write(9,2000)s80 read(8,905)j,(idum,i=1,j) 905 format(10i5) if(j.ne.nmo)call report('inconsistent number of orbitals') allocate(rox(grid(3)),orb(grid(3)*nmo),roy(grid(3)),roz(grid(3)), 1 ro(grid(3))) if(lw)then if(.not.lgeo.or.nat.lt.1)call report('Geometry unknown') do 2006 i=1,nat 2006 v(i)=0.0d0 endif db=(DD/bohr)**2 xp=ax-px do 906 ix=1,grid(1) xp=xp+px yp=ay-py do 906 iy=1,grid(2) yp=yp+py c do 904 i=1,grid(3) rox(i)=0.0d0 roy(i)=0.0d0 904 roz(i)=0.0d0 read(8,4000)((orb(nmo*(iz-1)+io),io=1,nmo),iz=1,grid(3)) c do 907 ie=1,n call setv(eau(ie),wtem,fr,fi,gamma,0) uux=dv(ie,1)*eau(ie)*fi uuy=dv(ie,2)*eau(ie)*fi uuz=dv(ie,3)*eau(ie)*fi do 907 j=1,ni(ie) ucx=cij(nmo**2*(ie-1)+j)*uux ucy=cij(nmo**2*(ie-1)+j)*uuy ucz=cij(nmo**2*(ie-1)+j)*uuz oa=aij(nmo**2*(ie-1)+j) ob=bij(nmo**2*(ie-1)+j) do 907 iz=1,grid(3) rox(iz)=rox(iz)+ucx*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) roy(iz)=roy(iz)+ucy*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) 907 roz(iz)=roz(iz)+ucz*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) c do 908 iz=1,grid(3) 908 ro(iz)=dsqrt(2.0d0*(rox(iz)**2+roy(iz)**2+roz(iz)**2)) if(lw)then do 2007 l=1,nat zp=az-pz do 2007 iz=1,grid(3) zp=zp+pz dist=((xp-x(1,l))**2+(yp-x(2,l))**2+(zp-x(3,l))**2)/db 2007 if(dist.lt.20.0d0)v(l)=v(l)+ro(iz)*exp(-dist) endif 906 write(9,4000)(ro(iz),iz=1,grid(3)) 4000 format(6E13.5) if(lw)then sum=0.0d0 do 2008 i=1,nat 2008 sum=sum+v(i) do 2009 i=1,nat 2009 v(i)=v(i)/sum write(6,6071)(100.0d0*v(i),i=1,nat) 6071 format(12f6.1) c recalculate to get rid of small contributions: allocate(ind(nat),jnd(nat)) call sort(v,ind,jnd,nat) call trimv(v,nat,nt) write(6,6071)(100.0d0*v(jnd(i)),i=1,nat) open(50,file='WEIGHTS.TXT') write(50,6071)(100.0d0*v(jnd(i)),i=1,nat) close(50) deallocate(ind,jnd) endif deallocate(rox,orb,roy,roz,ro) close(8) 909 close(9) write(6,*)' ro_alpha.cub written' endif c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp stop end subroutine trimv(v,nat,nt) implicit none real*8 v(*),t integer*4 nat,i,nt t=0.0d0 do 1 nt=1,nat t=t+v(nt) 1 if(t.ge.0.5d0)goto 2 2 if(nt.gt.nat)nt=nat do 4 i=1,nt 4 v(i)=v(i)/t do 3 i=nt+1,nat 3 v(i)=0.0d0 return end subroutine sort(v,ind,jnd,nat) implicit none real*8 v(*),t integer*4 ind(*),nat,i,j,it,jnd(*) do 1 i=1,nat jnd(i)=i 1 ind(i)=i do 2 i=1,nat-1 do 2 j=i+1,nat if(v(j).gt.v(i))then t=v(i) v(i)=v(j) v(j)=t it=ind(i) ind(i)=ind(j) ind(j)=it jnd(ind(j))=j jnd(ind(i))=i endif 2 continue return end c subroutine titles(ti,ext,out) character*(*) ti,ext character*80 out do 1 i=1,80 if (ti(i:i).eq.'.'.or.ti(i:i).eq.' ') then k=i-1 goto 2 endif 1 continue 2 out=ti(1:k)//ext return end subroutine setv(eau,w,fr,fi,g,ic) implicit none real*8 w,jm,fr,fi,g,t,wau,eau integer*4 ic c nm to au:= wau=(1.0d7/w)/219474.63d0 jm=((eau+wau)*(eau-wau))**2+(g*eau**2)**2 if(ic.eq.0)then t=eau else t=1.0d0 endif fr= (eau**2-wau**2) / jm * t fi= +g*eau**2 / jm * t fr=fr+fr fi=fi+fi return end subroutine getv(nm0,d,ds1,r0,angle,rs,r,i) implicit none integer*4 i,nm0 real*8 d(nm0,3),ds1,r0(*),angle,de,dm,ds,r(nm0,3),rs ds1=(d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3))*2.541765d0**2 rs=r0(i)*1.0d-4 c rs1=(d(i,1)*r(i,1)+d(i,2)*r(i,2)+d(i,3)*r(i,3)) c 1*2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 c c angle between electric and magnetic moment: de=dsqrt(d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3)) dm=dsqrt(r(i,1)*r(i,1)+r(i,2)*r(i,2)+r(i,3)*r(i,3)) ds= (d(i,1)*r(i,1)+d(i,2)*r(i,2)+d(i,3)*r(i,3)) angle=0.0d0 if(de*dm.gt.0.0d0)angle=acos(ds/de/dm)*180.0d0/4.0d0/atan(1.0d0) return end subroutine report(s) character*(*) s write(6,*)s stop 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 rwm(b,n) implicit none integer*4 i,j,k,n,n2 real*8 b(*) real*8 ,allocatable::p(:,:) character*2 ss(9) data ss/'PX','PY','PZ','VX','VY','VZ','LX','LY','LZ'/ allocate(p(n,n)) n2=n**2 do 1 i=1,9 open(38,file=ss(i)//'.MO.SCR.TXT') call rmtr38(p,n,n,1) close(38) do 1 j=1,n do 1 k=1,n 1 b(n2*(i-1)+n*(k-1)+j)=p(k,j) 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(9) data ss/'PBX','PBY','PBZ','VBX','VBY','VBZ','LBX','LBY','LBZ'/ allocate(p(n,n)) n2=n**2 do 1 i=1,9 open(38,file=ss(i)//'.MO.SCR.TXT') call rmtr38(p,n,n,1) close(38) do 1 j=1,n do 1 k=1,n 1 b(n2*(i-1)+n*(k-1)+j)=p(k,j) return end FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF(I.EQ.1)then IF(J.EQ.2.AND.K.EQ.3)then EPS= 1.0d0 else IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 endif else IF(I.EQ.2)then IF(J.EQ.3.AND.K.EQ.1)then EPS= 1.0d0 else IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 endif else IF(I.EQ.3)then IF(J.EQ.1.AND.K.EQ.2)then EPS= 1.0d0 else IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 endif endif endif endif RETURN END subroutine getdipoles(ut,nm0,n,cij,nmo,n2,ni,bp,aij,bij,u0, 1dij,nid,bpd,daij,dbij,u0d,lopen) c permanent dipoles of excited states implicit none integer*4 nm0,n,ie,j,ni(*),i,ip,oc,od,oa,ob,aij(*),bij(*),nmo, 1n2,nid(*),daij(*),dbij(*) real*8 ut(3,nm0),cij(*),bp(*),tq,cijip,ciji,u0(*),dij(*),bpd(*), 1u0d(*) logical lopen tq=dsqrt(2.0d0) do 2 ie=1,n do 2 j=1,3 c : ut(j,ie)=0.0d0 if(.not.lopen)then do 1 ip=1,ni(ie) oc=aij(n2*(ie-1)+ip) od=bij(n2*(ie-1)+ip) cijip=cij(n2*(ie-1)+ip) do 1 i=1,ni(ie) oa=aij(n2*(ie-1)+i) ob=bij(n2*(ie-1)+i) ciji=cij(n2*(ie-1)+i)*cijip if(oa.ne.oc)then if(ob.eq.od)ut(j,ie)=ut(j,ie)+ciji*bp(oc+nmo*(oa-1)+n2*(j-1))*tq else if(ob.ne.od)then ut(j,ie)=ut(j,ie)+ciji*bp(ob+nmo*(od-1)+n2*(j-1))*tq else ut(j,ie)=ut(j,ie)+ciji*(u0(j)-bp(oa+nmo*(oa-1)+n2*(j-1)) 1 + bp(ob+nmo*(ob-1)+n2*(j-1))) endif endif 1 continue else do 3 ip=1,ni(ie) oc=aij(n2*(ie-1)+ip) od=bij(n2*(ie-1)+ip) cijip=cij(n2*(ie-1)+ip) do 3 i=1,ni(ie) oa=aij(n2*(ie-1)+i) ob=bij(n2*(ie-1)+i) ciji=cij(n2*(ie-1)+i)*cijip if(oa.ne.oc)then if(ob.eq.od)ut(j,ie)=ut(j,ie)+ciji*bp(oc+nmo*(oa-1)+n2*(j-1)) else if(ob.ne.od)then ut(j,ie)=ut(j,ie)+ciji*bp(ob+nmo*(od-1)+n2*(j-1)) else ut(j,ie)=ut(j,ie)+ciji*(u0(j)-bp(oa+nmo*(oa-1)+n2*(j-1)) 1 + bp(ob+nmo*(ob-1)+n2*(j-1))) endif endif 3 continue do 4 ip=1,nid(ie) oc=daij(n2*(ie-1)+ip) od=dbij(n2*(ie-1)+ip) cijip=dij(n2*(ie-1)+ip) do 4 i=1,nid(ie) oa=daij(n2*(ie-1)+i) ob=dbij(n2*(ie-1)+i) ciji=dij(n2*(ie-1)+i)*cijip if(oa.ne.oc)then if(ob.eq.od)ut(j,ie)=ut(j,ie)+ciji*bpd(oc+nmo*(oa-1)+n2*(j-1)) else if(ob.ne.od)then ut(j,ie)=ut(j,ie)+ciji*bpd(ob+nmo*(od-1)+n2*(j-1)) else ut(j,ie)=ut(j,ie)+ciji*(u0d(j)-bpd(oa+nmo*(oa-1)+n2*(j-1)) 1 + bpd(ob+nmo*(ob-1)+n2*(j-1))) endif endif 4 continue endif 2 continue return end