program qvibr c simplified vibronic spectra implicit none integer*4 NQ,i,N,j,l,imax,m,ii,ia,mmax,iwr,ma,nat,NROOT, 1np,ix,np0,ip,vblock,nb parameter (mmax=20,nb=5) integer*4 sj(mmax+1),sit(mmax+3),ic(mmax+1) real*8 smax,CM,ww,kk,x,oo,on,xn,ai,tol1,is,aa,prod,ss, 1u0(3),m0(3),v0(3),q0(6),we,no,lmin,lmax,dl,u(3),r(3),FC,pp, 1rs,ds,t,fwcm,fwnm,ECM,eau,lw,oo1,tol2,Enm,tol3, 1GPT1(1,1),C1(1,1),D1(1),on1,qkf,jj,debye,clight, 1oqn1,oqn0,e00,dbig(nb),rbig(nb) c sij integer*4,allocatable:: nmax(:),na(:),ind(:),mo(:,:),mods(:,:), 1mors(:,:) real*8,allocatable:: A(:,:),B(:),C(:,:),D(:),E(:,:),JT(:,:),K(:), 1WP(:),W(:),sg(:,:),se(:,:),ov(:,:),dd(:),vv(:),rr(:),qq(:),ddi(:), 1vvi(:),rri(:),qqi(:),grad(:),s1(:),s2(:),oqn(:,:),G(:,:),GP(:,:), 1F(:,:),FI(:,:) logical lglg,lab,lq1 character*10 s10 write(6,6000) 6000 format(/, 1' Simplified vibronic spectra',/,/, 1' Input: VIBR.OPT options',/, 1' DUSCH.OUT normal mode transformations',/, 1' G.OUT Gaussian output with transition', 1'parameters',/, 1' Oputput: uvvibr.prn, ecdvibr.prn, the spectra',/) call readopt(lglg,lab,lq1,CM,tol1,tol2,lmin,lmax,np,fwcm, 1np0,e00,iwr,dl,vblock) sit=0 call nn(N,NQ) nat=N/3 allocate(JT(NQ,NQ),K(NQ),A(NQ,NQ),B(NQ),C(NQ,NQ),D(NQ),E(NQ,NQ), 1WP(NQ),W(NQ),sg(N,N),se(N,N),nmax(NQ),na(NQ),ov(NQ,mmax+1), 1dd(3*N),rr(3*N),vv(3*N),qq(6*N),ddi(3*N),rri(3*N),vvi(3*N), 1qqi(6*N),grad(N),ind(NQ),s1(np),s2(np),oqn(NQ,mmax+1), 1G(NQ,NQ),GP(NQ,NQ),F(NQ,NQ),FI(NQ,NQ),mo(NQ,mmax+1), 1mods(nb,NQ),mors(nb,NQ)) debye=2.541d0 clight=137.03599d0 mods=0 mors=0 dbig=0.0d0 rbig=0.0d0 oqn=0.0d0 ov=0.0d0 s1=0.0d0 s2=0.0d0 nmax=0 ov=0.0d0 is=0.0d0 c Q = J Q' + K call rdusch(NQ,JT,K,A,B,C,D,E,F,FI,W,WP,G,GP,'DUSCH.OUT') call rduschs(N,NQ,sg ,'DUSCH.OUT','Ground S-Matrix') call rduschs(N,NQ,se,'DUSCH.OUT','Excited S-Matri') c find which excited state is calculated call chkroot('G.OUT',NROOT) write(6,608)NROOT 608 format(' ROOT = ',i3) c read moments, energy: call r00('G.OUT',NROOT,u0,m0,v0,q0,we) c read moment derivatives: call rdd('G.OUT',nat,NROOT,u0,v0,m0,q0,dd,vv,rr,qq,grad) C transform into ground state normal modes call trafN(3,dd,ddi,nat,NQ,sg) call trafN(3,rr,rri,nat,NQ,sg) call trafN(6,qq,qqi,nat,NQ,sg) call trafN(3,vv,vvi,nat,NQ,sg) if(vblock.ne.0)write(6,607)vblock,W(vblock)*CM 607 format(i4,' blocked states, until',f10.1,' cm-1') c 111111111111111111111111111111111111111111111111111111111111 c one-dimensional integrals: do 1 i=vblock+1,NQ imax=1 smax=0.0d0 do 2 j=1,NQ if(dabs(JT(j,i)).gt.smax)then smax=dabs(JT(j,i)) imax=j endif 2 continue c sij=0.0d0 c do 3 l=1,N c sij=sij+sg(l,i)*se(l,j) c if(abs(sij).gt.smax)then c smax=abs(sij) c imax=j c endif c continue c do 2 j=1,NQ c sij=0.0d0 c do 3 l=1,N c sij=sij+sg(l,i)*se(l,j) c if(abs(sij).gt.smax)then c smax=abs(sij) c imax=j c endif c continue ww=W(i) pp=WP(imax) kk=K(i) x=ww*kk**2/2.0d0 oo=exp(-x/2.0d0) c <0|0'>: jj=abs(JT(imax,i)) c jj=1.0d0 t=jj*ww*jj+pp oo1=dsqrt(2.0d0*dsqrt(ww*pp)*jj/t)*exp(x*(jj*ww*jj/t-1.0d0)) C1(1,1)=2.0d0*pp/t-1.0d0 D1(1)=-2.0d0*dsqrt(pp)*ww*jj*kk/t GPT1(1,1)=GP(imax,imax) do 4 m=0,mmax do 42 ii=1,m 42 sj(ii)=1 c <_0|m>, xn = x^n/n!: xn=1.0d0 ai=0.0d0 do 41 ii=1,m ai=ai+1.0d0 41 xn=xn*x/ai no= dsqrt(xn)*oo if(mod(m,2).ne.0)then on=-no else on= no endif c <0|Q*k|*>=sqrt(h/(2*wk)) [ sqrt(vk) <0|k-1> + sqrt(vk+1) <0|k+1> ]: on1=FC(oo1,sj,m,np0,0,C1,D1,1) c <0|Q'|*> = <0|Q|*)-K <0|*>: oqn0=-0.5d0*kk*(1.0d0+dble(m)/x)*on oqn1=qkf(sj,m,1,oo1,sit,0,C1,D1,1,np0,GPT1) if(lab)then ov( i,m+1)=on1 oqn(i,m+1)=oqn1 else ov( i,m+1)=on oqn(i,m+1)=oqn0 endif 4 if(iwr.gt.1)write(6,601)m,on,on1,oqn0,oqn1 601 format(i4,6f12.6) c c select biggest,order by size ic=0 nmax(i)=0 4001 t=0.0d0 ip=0 do 43 ii=0,m if(dabs(ov(i,ii+1)).gt.dabs(t).and.ic(ii+1).eq.0)then t=ov(i,ii+1) ip=ii endif 43 continue if(ic(ip+1).eq.0.and.dabs(t).gt.tol1)then nmax(i)=nmax(i)+1 mo(i,nmax(i))=ip ic(ip+1)=1 goto 4001 endif ss=0.0d0 do 44 ii=1,nmax(i) on=ov(i,mo(i,ii)+1) if(iwr.gt.0)write(6,602)ii,mo(i,ii),on,oqn(i,mo(i,ii)+1) 602 format(2i4,2f12.6) 44 ss=ss+on*on if(nmax(i).eq.0)then do 441 ii=0,mmax 441 write(6,601)ii,ov(i,ii+1),oqn(i,ii+1) write(6,*)'mode ',i call report('no state left') endif is=is+log(dble(nmax(i))) ind(i)=imax if(iwr.gt.0)write(6,600)i,W(i)*CM,x,imax,smax,nmax(i),ss 600 format(i4,f11.2,e12.4,i4,e12.4,i3,f10.4) c renormalize: do 1 ii=1,nmax(i) 1 ov(i,mo(i,ii)+1)=ov(i,mo(i,ii)+1)/dsqrt(ss) c 111111111111111111111111111111111111111111111111111111111111 is=is/log(10.0d0) ma=int(is) aa=is-dble(ma) write(s10,771)ma 771 format(i10) do 8 ii=1,len(s10) 8 if(s10(ii:ii).ne.' ')goto 81 81 write(6,774)10.0d0**aa 774 format(/,f5.2,'x10^',$) do 711 i=ii,len(s10) 711 write(6,772)s10(i:i) 772 format(A1,$) write(6,7731)tol1 7731 format(' states, tol1:',e9.2,/) if(e00.gt.1.0d-40)then eau=e00 else eau=we endif c prod = <0|e> = <0|n1><0|n2> ... <0|nNQ> prod=1.0d0 do 801 i=vblock+1,NQ eau=eau+mo(i,1)*WP(ind(i)) 801 prod=prod*ov(i,mo(i,1)+1) write(6,603)prod,tol2 603 format(' Initial product:', e12.4,' tol2:',e9.2) tol3=tol2*prod j=0 ss=0.0d0 na=1 5 ia=vblock+1 if(NQ.eq.3)write(6,666)(mo(l,na(l)),l=1,NQ),prod 666 format(3i2,e12.4) j=j+1 ss=ss+prod*prod c <0|u|na> = u0 <0|na> + uq <0|q|na> c = -<0|m|e>: u= u0*prod r=-m0*prod if(lq1)then do 9 i=vblock+1,NQ t=prod*oqn(i,mo(i,na(i))+1)/ov(i,mo(i,na(i))+1) ip=3*(ind(i)-1) c du0/dQ <0|Qi|mi=na(i)+1> <0|na without i> do 9 ix=1,3 u(ix)=u(ix)+ddi(ix+ip)*t 9 r(ix)=r(ix)-rri(ix+ip)*t endif ds=u(1)*u(1)+u(2)*u(2)+u(3)*u(3) rs=u(1)*r(1)+u(2)*r(2)+u(3)*r(3) ECM=eau*CM Enm=1.0d7/ECM fwnm=fwcm*Enm/ECM call app(Enm,ds*Enm,rs*Enm,s1,s2,lmin,lmax,np,fwnm,lglg) call maxdd(nb,vblock,NQ,ds,rs,dbig,rbig,mods,mors,mo,na,mmax) 6 if(na(ia).lt.nmax(ia).and.abs(prod).gt.tol3)then na(ia)=na(ia)+1 prod=prod*ov(ia,mo(ia,na(ia))+1)/ov(ia,mo(ia,na(ia)-1)+1) eau=eau+dble(mo(ia,na(ia))-mo(ia,na(ia)-1))*WP(ind(ia)) if(abs(prod).lt.tol3)then prod=prod*ov(ia,mo(ia,1)+1)/ov(ia,mo(ia,na(ia))+1) eau=eau+dble(mo(ia,1)-mo(ia,na(ia)))*WP(ind(ia)) na(ia)=1 ia=ia+1 if(ia.le.NQ)goto 6 else goto 5 endif else prod=prod*ov(ia,mo(ia,1)+1)/ov(ia,mo(ia,na(ia))+1) eau=eau+dble(mo(ia,1)-mo(ia,na(ia)))*WP(ind(ia)) na(ia)=1 ia=ia+1 if(ia.le.NQ)goto 6 endif write(6,775)j,tol2,ss 775 format(i20,' states selected, tol2:',e12.4,/, 1' <0|n> = ',e12.4) s1=s1/ss*108.7d0*debye**2 s2=s2/ss*435.0d0*debye**2/clight open(45,file='uvvib.prn') open(46,file='ecdvib.prn') lw=lmin-dl do 10 i=1,np lw=lw+dl write(45,451)lw,s1(i) 10 write(46,451)lw,s2(i) 451 format(f12.2,e12.4) close(45) close(46) call bc('D',nb,vblock,NQ,dbig,mods,WP,CM,e00) call bc('R',nb,vblock,NQ,rbig,mors,WP,CM,e00) end subroutine rduschs(N,NQ,s,f,ff) IMPLICIT none character*(*) f,ff integer*4 NQ,N,ic,I,J real*8 s(N,N),SFAC character*80 s80 SFAC=0.0234280d0 open(20,file=f,status='old') ic=0 read(20,*)NQ 2 read(20,900,end=88,err=88)s80 900 format(a80) if(s80(2:16).eq.ff)then call rds(20,NQ,N,s,ic) DO 5 I=1,N DO 5 J=1,NQ 5 s(I,J)=s(I,J)*SFAC close(20) write(6,*)f write(6,*)'S read, ',NQ,' modes' return endif goto 2 88 close(20) write(6,*)f if(ic.eq.0)call report('S not found') end subroutine rds(io,NQ,N,M,ic) IMPLICIT none integer*4 NQ,N,N1,N3,io,LN,J,ic real*8 M(N,N) ic=ic+1 read(io,*) N1=1 1 N3=min(N1+4,NQ) read(io,*) DO 130 LN=1,N 130 READ(io,*)M(LN,N1),(M(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.NQ)GOTO 1 return end subroutine rdusch(NQ,JT,K,A,B,C,D,E,FF,FI,W,WP,G,GP,f) c Q = JQ' + K c JT = transposed J IMPLICIT none character*(*) f integer*4 NQ,ip(10),iw,i,IERR real*8 A(NQ,NQ),B(NQ),C(NQ,NQ),D(NQ),E(NQ,NQ),JT(NQ,NQ),K(NQ), 1WP(NQ),W(NQ),G(NQ,NQ),GP(NQ,NQ),FF(NQ,NQ),FI(NQ,NQ) real*8 ,allocatable::T(:,:),TP(:,:),J(:,:),EE(:,:) character*80 s80 allocate(T(NQ,NQ),TP(NQ,NQ),J(NQ,NQ),EE(NQ,NQ+NQ)) ip=0 iw=0 open(20,file=f,status='old') read(20,*) 2 read(20,900,end=88,err=88)s80 900 format(a80) if(s80(2:18).eq.'Duschinsky Matrix') 1 call rdm(20,NQ,JT,ip(1),s80) if(s80(2:13).eq.'Shift Vector') call rdv(20,NQ,K ,ip(2)) if(s80(2: 9).eq.'A Matrix' ) call rdm(20,NQ,A ,ip(3),s80) if(s80(2: 9).eq.'B Vector' ) call rdv(20,NQ,B ,ip(4)) if(s80(2: 9).eq.'C Matrix' ) call rdm(20,NQ,C ,ip(5),s80) if(s80(2: 9).eq.'D Vector' ) call rdv(20,NQ,D ,ip(6)) if(s80(2: 9).eq.'E Matrix' ) call rdm(20,NQ,E ,ip(7),s80) if(s80(2: 9).eq.'Ground w' ) call rdw(20,NQ,W ,iw) if(s80(2:10).eq.'Excited w') call rdw(20,NQ,WP,iw) goto 2 88 close(20) if(ip(1).eq.0)call report('Duschinsky Matrix not found') if(ip(2).eq.0)call report('Shift Vector not found') if(ip(3).eq.0)call report('A Matrix not found') if(ip(4).eq.0)call report('B Vector not found') if(ip(5).eq.0)call report('C Matrix not found') if(ip(6).eq.0)call report('D Vector not found') if(ip(7).eq.0)call report('E Matrix not found') if(iw.lt.2)call report('Frequencies not found') G=0.0d0 GP=0.0d0 do 1 i=1,NQ G(i,i)=W(i) 1 GP(i,i)=WP(i) call mt(J,JT,NQ) call mm(T,G,J,NQ) call mm(TP,JT,T,NQ) call ms(FF,TP,GP,NQ) call INV(FF,FI,NQ,EE,IERR) if(IERR.ne.0)call report('Inversion error') write(6,*)'DUSCH.OUT read' return end subroutine rdv(io,N,V,ic) IMPLICIT none integer*4 N,io,LN,ic real*8 V(N) ic=ic+1 read(io,*) read(io,*) DO 130 LN=1,N 130 READ(io,*)V(LN),V(LN) return end subroutine rdm(io,NQ,M,ic,s) IMPLICIT none integer*4 NQ,N1,N3,io,LN,J,ic real*8 M(NQ,NQ) character*(*) s write(6,*)s(1:len(s)/2) ic=ic+1 read(io,*) N1=1 1 N3=min(N1+4,NQ) read(io,*) DO 130 LN=1,NQ 130 READ(io,*)M(LN,N1),(M(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.NQ)GOTO 1 return end subroutine rdw(io,N,e,iw) IMPLICIT none integer*4 i,io,iw,N real*8 e(N),CM read(io,*) read(io,100)(e(i),i=1,N) 100 format(6f11.2) CM=219474.63d0 do 1 i=1,N 1 e(i)=e(i)/CM iw=iw+1 return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine nn(N,NQ) implicit none integer*4 NQ,i,N character*80 s80 N=0 open(9,file='DUSCH.OUT',status='old') read(9,*)NQ 1 read(9,80,end=99,err=99)s80 80 format(a80) if(s80(2:16).eq.'Ground S-Matrix')then read(9,*) read(9,*) 2 read(9,7,end=99,err=99)i 7 format(i7) if(i.eq.0)goto 99 N=N+1 goto 2 endif goto 1 99 close(9) write(6,*)NQ,' modes,',N/3,' atoms' return end subroutine trafN(N,dd,ddi,nat,NQ,s) implicit none integer*4 N,nat,NQ,ix,iq,ia,xa,iqi,ii real*8 dd(*),ddi(*),s(3*nat,3*nat) do 1 ix=1,N do 1 iq=1,NQ iqi=ix+N*(iq-1) ddi(iqi)=0.0d0 do 2 ia=1,nat do 2 xa=1,3 ii=xa+3*(ia-1) c s-opposite order of normal modes: 2 ddi(iqi)=ddi(iqi)+dd(ix+N*(ii-1))*s(ii,iq) 1 continue return end subroutine app(ECM,r,a,reram,reroa,wrin,wrax,npx,fwhh,lglg) implicit none real*8 ECM,wrin,wrax,reram(*),reroa(*),fwhh,w,dw,t,sf, 1r,a integer*4 npx,i logical lglg c Make the scale similar to usual new6/tabprnf way c dw=(wrax-wrin)/(npx-1) w=wrin-dw do 1 i=1,npx w=w+dw t=sf(fwhh,lglg,w,ECM) reram(i)=reram(i)+t*r 1 reroa(i)=reroa(i)+t*a return end function sf(d,g,x,x0) implicit none logical g real*8 pi,spi,sf,d,x,x0,dd pi =3.14159265358979d0 spi=1.77245385090552d0 dd=((x-x0)/d)**2 if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.1.0d10)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif return end subroutine chkroot(f,NROOT) implicit none integer*4 NROOT,found,nd,i character*(*) f character*80 s nd=0 found=0 open(9,file=f,status='old') 1 read(9,900,end=88,err=88)s 900 format(a80) if(s(2:15).eq.'Excited State '.and.s(16:16).ne.'s')then do 20 i=1,79 20 if(s(i:i) .eq.':') read(s(15:i-1),*)nd endif if(s(2:29).eq.'This state for optimization ')found=nd goto 1 88 close(9) write(6,600)found,nd 600 format(i6,' of',i6,' excited states in ',$) do 21 i=1,len(f) 21 if(f(i:i).ne.' ')write(6,601)f(i:i) 601 format(a1,$) write(6,*) NROOT=found return end subroutine r00(f,NROOT,u0,m0,v0,q0,eau) c read transition dipoles from excited state freq calc implicit none integer*4 ln,i,NROOT,ic,ix,nd,i1,i2,i3,i4,i5 real*8 u0(*),m0(3),v0(3),q0(6),ev,eau character*(*) f character*80 s80 real*8,allocatable::v16(:,:) c write(6,*) write(6,*)' R00: reading '//f write(6,*) open(9,file=f,status='old') ln=0 i1=0 i2=0 i3=0 i4=0 i5=0 c where it starts: 1 read(9,900,end=88,err=88)s80 900 format(a80) ln=ln+1 if(s80(1:2).eq.' #')then ic=0 do 3 i=1,len(s80)-4 if(s80(i:i+4).eq.' freq'.or.s80(i:i+4).eq.' freq')ic=ic+1 3 if(s80(i:i+2).eq.' td' .or.s80(i:i+2).eq.' TD' )ic=ic+1 if(ic.eq.2)then write(6,*)'TD starts at line ',ln goto 4 endif endif goto 1 88 close(9) call report('TD freq not found') 4 read(9,900,end=78,err=78)s80 if(s80(2:65).eq.'Ground to excited state transition electric dipol 1e moments (Au):'.and.i1.eq.0)then write(6,900)s80 do 5 i=1,NROOT 5 read(9,*) read(9,*)u0(1),(u0(ix),ix=1,3) i1=i1+1 endif if(s80(2:65).eq.'Ground to excited state transition velocity dipol 1e moments (Au):'.and.i2.eq.0)then write(6,900)s80 do 11 i=1,NROOT 11 read(9,*) read(9,*)v0(1),(v0(ix),ix=1,3) i2=i2+1 endif if(s80(2:65).eq.'Ground to excited state transition magnetic dipol 1e moments (Au):'.and.i3.eq.0)then write(6,900)s80 do 14 i=1,NROOT 14 read(9,*) read(9,*)m0(1),(m0(ix),ix=1,3) i3=i3+1 endif if(s80(2:69).eq.'Ground to excited state transition velocity quadr 1upole moments (Au):'.and.i4.eq.0)then write(6,900)s80 do 17 i=1,NROOT 17 read(9,*) read(9,*)q0(1),q0(1),q0(4),q0(6),q0(2),q0(3),q0(5) i4=i4+1 endif if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s'. 1and.i5.eq.0)then do 20 i=1,79 if(s80(i:i) .eq.':') read(s80(15:i-1),*)nd 20 if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev if(nd.eq.NROOT)then write(6,900)s80 eau=ev/27.211384205943d0 i5=i5+1 endif endif c analytical frequency output: c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA if(s80(2:31).eq.'Electronic transition elements')then c read NROOT moments write(6,900)s80 allocate(v16(16,NROOT)) call ru(9,v16,16,NROOT) do 21 ix=1,3 u0(ix )=v16(ix+ 1,NROOT) v0(ix )=v16(ix+ 4,NROOT) 21 m0(ix )=v16(ix+ 7,NROOT) q0(1)=v16(11,NROOT) q0(2)=v16(14,NROOT) q0(3)=v16(15,NROOT) q0(4)=v16(12,NROOT) q0(5)=v16(16,NROOT) q0(6)=v16(13,NROOT) c our order 1 2 3 4 5 6 c our order xx xy xz yy yz zz c 1 2 3 4 5 6 c 11 12 13 14 15 16 c gaussian # xx yy zz xy xz yz deallocate(v16) write(6,*)'Analytical transitions read' endif c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA goto 4 78 close(9) if(nd.eq.0)call report('Energy not found') v0 = -v0 / eau q0 = -q0 / eau m0 = m0 / 2.0d0 c looks like that we have this situation: c Gaussian gives <0|r|e>, , <0|-r x grad|e> c and c To transform it to what we want we can use c hbar^2 c <0|r|e> = ------- <0|grad|e> c m Ee0 c and c hbar^2 c <0|rr|e> = ------- <0|r grad + grad r|e> c m E0 c Eeo=Ee-E0 c so that (in atomic units hbar=m=1) c <0|r|e>(in velocity) = - / Ee0 c Im <0|m|e>= Im <0|-r x grad/2|e> = <0|-r x grad|e> / 2 c <0|rr|e> = - / Ee0 c return end subroutine rdd(f,nat,NROOT,u0,v0,m0,q0,d,v,a,q,grad) c read transition dipole derivatives from excited state freq calc implicit none integer*4 nat,ia,xa,i,NROOT,ii,ix,nd,ntr,j,ir real*8 d(9*nat),u0(*),m0(3),a(9*nat), 1v(9*nat),v0(3),q(18*nat),q0(6),ev,eau,grad(3*nat) character*(*) f character*80 s80 real*8,allocatable::v16(:,:) c write(6,*) write(6,*)' Rdd: reading '//f write(6,*) c total number of roots (exc. states) calculated: ntr=0 open(9,file=f,status='old') 1 read(9,900,end=88,err=88)s80 900 format(a80) if(s80(38:59).eq.'Forces (Hartrees/Bohr)')then read(9,*) read(9,*) do 2 i=1,nat 2 read(9,*)(grad(ix+3*(i-1)),ix=1,2),(grad(ix+3*(i-1)),ix=1,3) endif if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then do 20 i=1,79 if(s80(i:i) .eq.':') read(s80(15:i-1),*)nd 20 if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev if(nd.eq.NROOT)then eau=ev/27.211384205943d0 endif endif if(s80(2:31).eq.'Electronic transition elements')then c read NROOT moments allocate(v16(16,NROOT)) call ru(9,v16,16,NROOT) do 21 ix=1,3 u0(ix)=v16(ix+ 1,NROOT) v0(ix)=v16(ix+ 4,NROOT) 21 m0(ix)=v16(ix+ 7,NROOT) q0(1 )=v16(11,NROOT) q0(2 )=v16(14,NROOT) q0(3 )=v16(15,NROOT) q0(4 )=v16(12,NROOT) q0(5 )=v16(16,NROOT) q0(6 )=v16(13,NROOT) deallocate(v16) write(6,*)'Analytical transitions read' endif if(s80(2:34).eq.'Electronic Transition Derivatives')then allocate(v16(16,3 + 3*nat)) call ru(9,v16,16,3 + 3*nat) do 22 ia=1,nat do 22 xa=1,3 ii=xa+3*(ia-1) do 221 ix=1,3 d(ix +3*(ii-1))=v16(ix+ 1,3+ii) v(ix +3*(ii-1))=v16(ix+ 4,3+ii) 221 a(ix +3*(ii-1))=v16(ix+ 7,3+ii) q(1+6*(ii-1))=v16(11,3+ii) q(2+6*(ii-1))=v16(14,3+ii) q(3+6*(ii-1))=v16(15,3+ii) q(4+6*(ii-1))=v16(12,3+ii) q(5+6*(ii-1))=v16(16,3+ii) 22 q(6+6*(ii-1))=v16(13,3+ii) deallocate(v16) endif c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA c checkpoint style c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC if(s80(1:9).eq.'ETran sym')read(s80(50:61),*)ntr if(s80(1:18).eq.'ETran state values')then if(ntr.eq.0)call report('nstates not found') allocate(v16(1,16*ntr+48+16*3*nat)) read(9,*)v16 ii=0 do 237 ir=1,ntr ii=ii+1 do 231 ix=1,3 ii=ii+1 231 if(ir.eq.NROOT)u0(ix)=v16(1,ii) do 232 ix=1,3 ii=ii+1 232 if(ir.eq.NROOT)v0(ix)=v16(1,ii) do 233 ix=1,3 ii=ii+1 233 if(ir.eq.NROOT)m0(ix)=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(1 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(4 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(6 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(2 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(3 )=v16(1,ii) ii=ii+1 237 if(ir.eq.NROOT)q0(5 )=v16(1,ii) ii=ii+48 do 238 j=1,3*nat ii=ii+1 do 234 ix=1,3 ii=ii+1 234 d(ix+3*(j-1))=v16(1,ii) do 235 ix=1,3 ii=ii+1 235 v(ix+3*(j-1))=v16(1,ii) do 236 ix=1,3 ii=ii+1 236 a(ix+3*(j-1))=v16(1,ii) ii=ii+1 q(1+6*(j-1))=v16(1,ii) ii=ii+1 q(4+6*(j-1))=v16(1,ii) ii=ii+1 q(6+6*(j-1))=v16(1,ii) ii=ii+1 q(2+6*(j-1))=v16(1,ii) ii=ii+1 q(3+6*(j-1))=v16(1,ii) ii=ii+1 238 q(5+6*(j-1))=v16(1,ii) deallocate(v16) write(6,*)'Analytical transition derivatives read from chk' endif c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC goto 1 88 close(9) if(nd.eq.0)call report('Energy not found') v = -v/eau q = -q/eau a = a/2.0d0 v0 = -v0/eau q0 = -q0/eau m0 = m0/2.0d0 return end subroutine ru(io,A,n,m) c matrix n x m implicit none REAL*8 A(n,m) INTEGER*4 io,n,m,N1,N3,LN,J N1=1 1 N3=MIN(N1+4,m) read(io,*) DO 130 LN=1,N 130 read(io,*)A(LN,N1),(A(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.m)GOTO 1 return end subroutine ms(A,B,C,NQ) c matrix sum A = B + C implicit none integer*4 NQ,i,j real*8 A(NQ,NQ),B(NQ,NQ),C(NQ,NQ) do 1 i=1,NQ do 1 j=1,NQ 1 A(i,j)=B(i,j)+C(i,j) return end subroutine mt(AT,A,NQ) c matrix transposition AT = A^T implicit none integer*4 NQ,i,j real*8 A(NQ,NQ),AT(NQ,NQ) do 1 i=1,NQ do 1 j=1,NQ 1 AT(i,j)=A(j,i) return end subroutine INV(a,ai,n,e,IERR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(n,n),a(n,n),e(n,2*n) C TOL=1.0d-10 10000 IERR=0 DO 1 ii=1,n DO 1 jj=1,n e(ii,jj)=a(ii,jj) 1 e(ii,jj+n)=0.0D0 do 13 ii=1,n 13 e(ii,ii+n)=1.0D0 c DO 2 ii=1,n-1 if (ABS(e(ii,ii)).LE.TOL) then DO 3 io=ii+1,n 3 if (ABS(e(io,ii)).GT.TOL) goto 11 IERR=1 write(6,*)ii write(6,*)'tol = ',tol tol=tol*0.50d0 if(tol.gt.1.0d-20)goto 10000 RETURN c 11 CONTINUE DO 4 kk=1,2*n w=e(ii,kk) e(ii,kk)=e(io,kk) 4 e(io,kk)=w ENDIF eii=e(ii,ii) DO 5 jj=ii+1,n e1=e(jj,ii)/eii DO 6 kk=ii+1, 2*n 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e1 5 e(jj,ii)=0.0D0 2 CONTINUE c DO 7 i2=n,2,-1 eii=e(i2,i2) DO 7 j2=i2-1,1,-1 e1=e(j2,i2)/eii DO 9 kk=1, n 9 e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1 7 e(j2,i2)=0.0d0 c DO 10 ii=1,n ei=1.0d0/e(ii,ii) DO 12 jj=1,n 12 ai(ii,jj)=e(ii,jj+n)*ei 10 CONTINUE c RETURN END subroutine mm(A,B,C,NQ) c matrix multiplication A = B x C implicit none integer*4 NQ,i,j,k real*8 A(NQ,NQ),B(NQ,NQ),C(NQ,NQ) do 1 i=1,NQ do 1 j=1,NQ A(i,j)=0.0d0 do 1 k=1,NQ 1 A(i,j)=A(i,j)+B(i,k)*C(k,j) return end c ========================= function FC(fac,si,Nexc,np0,iwr,C,D,NQ) c FC = Franc Condon factor for <0|si> c fac = <0|0*> c np0 ... dimension of working buffer implicit none integer*4 si(*),Nexc,np,iex,ii,nu,jj,ip,NQ,ic,jold, 1nuj,jc,kk,iwr,np0,iq,jq real*8 FC,fac,D(*),C(NQ,NQ),pini real*8,allocatable::p(:) integer*4,allocatable::NN(:),ie(:,:),it(:) np=1 allocate(p(np0),NN(np0),ie(np0,Nexc),it(np0)) p(np)=1.0d0 NN(np)=Nexc if(iwr.gt.1)write(6,604)(si(iex),iex=1,NExc) 604 format(/,'<0|e>, e = ',10i3,' requested:') do 102 iex=1,NN(np) 102 ie(np,iex)=si(iex) c expand reccurent formula into a sum and shrink back to <0|0> 777 do 101 ii=1,np c term ii - reduce excitation one by one: pini=p(ii) if(iwr.gt.1)then write(6,6071)ii,p(ii) 6071 format(' term ',i2,':',f10.6) write(6,6072)(ie(ii,iex),iex=1,NN(ii)) 6072 format(20i3) endif do 101 iex=1,NN(ii) iq=ie(ii,iex) if(iq.gt.0)then c <0| si_nu> reduce to <0| i-1>, <0|i-2>, <0|i-1,j-1>(j<>i) nu=0 do 1032 jj=1,NN(ii) 1032 if(ie(ii,jj).eq.iq)nu=nu+1 c write term <0|v'-1_iex> into string it: ip=0 ic=0 do 1031 jj=1,NN(ii) if(iq.eq.ie(ii,jj).and.ic.lt.1)then ic=ic+1 else ip=ip+1 it(ip)=ie(ii,jj) endif 1031 continue c call digest(np,ie,np0,Nexc,p,NN,it,NN(ii)-1, 1 pini*D(iq)/dsqrt(dble(2*nu))) if(iwr.gt.1)write(6,606)(it(jj),jj=1,NN(ii)-1) 606 format('digest ',10i3) if(nu.gt.1)then c write term <0|v'-2_iex> into string it: ic=0 ip=0 do 1033 jj=1,NN(ii) if(ie(ii,jj).eq.iq.and.ic.lt.2)then ic=ic+1 else ip=ip+1 it(ip)=ie(ii,jj) endif 1033 continue call digest(np,ie,np0,Nexc,p,NN,it,NN(ii)-2, 1 pini*dsqrt(dble(nu-1)/dble(nu))*C(iq,iq)) if(iwr.gt.1)write(6,606)(it(jj),jj=1,NN(ii)-2) endif jold=0 do 106 jj=1,NN(ii) jq=ie(ii,jj) if(jq.ne.iq.and.jq.ne.jold)then nuj=0 do 1034 kk=1,NN(ii) 1034 if(jq.eq.ie(ii,kk))nuj=nuj+1 c write term <0|v'-1_iex-1_j, j<>iex> into it: ic=0 jc=0 ip=0 do 1035 kk=1,NN(ii) if(ie(ii,kk).eq.iq.and.ic.lt.1)then ic=ic+1 else if(ie(ii,kk).eq.jq.and.jc.lt.1)then jc=jc+1 else ip=ip+1 it(ip)=ie(ii,kk) endif endif 1035 continue call digest(np,ie,np0,Nexc,p,NN,it,NN(ii)-2, 1 pini*dsqrt(dble(nuj)/dble(nu))*C(iq,jq)) if(iwr.gt.1)write(6,606)(it(kk),kk=1,NN(ii)-2) endif 106 jold=jq if(iwr.gt.1)then write(6,*)'np now ',np do 108 jj=1,np 108 write(6,605)jj,NN(jj),p(jj),(ie(jj,kk),kk=1,NN(jj)) endif c eliminate the old term and start over do 105 jj=ii,np-1 p(jj)=p(jj+1) NN(jj)=NN(jj+1) do 105 kk=1,NN(jj) 105 ie(jj,kk)=ie(jj+1,kk) np=np-1 if(iwr.gt.1)then write(6,*)'np after elimination ',np do 107 jj=1,np 107 write(6,605)jj,NN(jj),p(jj),(ie(jj,kk),kk=1,NN(jj)) 605 format(i3,i2,' ',g9.3,10i3) endif goto 777 endif 101 continue if(np.ne.1)call report('np <> 1') FC=p(1)*fac return end function qkf(sj,Nj,kk,fac,sit,iwr,C,D,N,np0,GP) c <0|Qk|*>=sqrt(h/(2wk))(sqrt(vk)<0|*-1k>+sqrt(vk+1)<0|*+1k>) implicit none integer*4 sj(*),sit(*),Nj,nuk,kk,iwr,ic,ip,ik,N,np0 real*8 qkf,pm,FC,fac,C(N,N),D(*),pp,GP(N,N) nuk=0 do 16 ic=1,Nj 16 if(sj(ic).eq.kk)nuk=nuk+1 if(nuk.gt.0)then ip=0 ik=0 do 17 ic=1,Nj if(sj(ic).eq.kk.and.ik.eq.0)then ik=ik+1 else ip=ip+1 sit(ip)=sj(ic) endif 17 continue pm=dsqrt(dble(nuk))*FC(fac,sit,Nj-1,np0,iwr,C,D,N) else pm=0.0d0 endif do 18 ic=1,Nj 18 sit(ic)=sj(ic) sit(Nj+1)=kk pp=dsqrt(dble(nuk+1))*FC(fac,sit,Nj+1,np0,iwr,C,D,N) qkf=(pm+pp)/dsqrt(2.0d0*GP(kk,kk)) return end subroutine digest(np,ie,np0,LEXCL,p,NN,it,iexc,T) implicit none integer*4 je,np,np0,LEXCL,NN(*),it(*),ie(np0,LEXCL),jj,iexc, 1jeje real*8 T,p(*) c does it exist already within 1 ... np?: je=jeje(np,iexc,it,NN,ie,np0,LEXCL) if(je.ne.0)then c the term already exist in the expansion as je^th - just add it p(je)=p(je)+T else c term not found - add as new term np=np+1 if(np.gt.np0)then write(6,*)np,np0 call report('digest - too many terms') endif p(np)=T NN(np)=iexc do 1 jj=1,iexc 1 ie(np,jj)=it(jj) endif return end function jeje(np,nx,it,NN,ee,np0,LEXCL) implicit none integer*4 np,nx,it(*),NN(*),np0,LEXCL,ee(np0,LEXCL),je,jeje,jj,ie je=0 do 104 jj=1,np if(nx.ne.NN(jj))goto 104 do 1041 ie=1,NN(jj) 1041 if(it(ie).ne.ee(jj,ie))goto 104 je=jj goto 1042 104 continue 1042 jeje=je return end subroutine readopt(lglg,lab,lq1,CM,tol1,tol2,lmin,lmax,np,fwcm, 1np0,e00,iwr,dl,vblock) implicit none integer*4 np,np0,iwr,vblock logical lglg,lab,lex,lq1 real*8 CM,tol1,lmin,lmax,fwcm,e00,dl,e00cm,tol2 character*4 key,k2 CM=219474.63d0 c use Gaussian bands (and not Lorantzian): lglg=.true. c overlap tolerance tol1=1.0d-4 c product tolerance tol2=1.0d-4 c minimal wavelength / nm: lmin=50.0d0 c maximal wavelength / nm: lmax=4000.0d0 c include dipole derivatives (Herzberg-Teller term): lq1=.true. c number of points np=3951 c peak width in cm-1: fwcm=160.0d0 c buffer for FC integral calculation: np0=10000 c 0 - 0' energy, in atomic units: e00cm=16744.549d0 c extended output: iwr=1 c use 1D Duschinsky transformation lab=.false. c number of blocked modes: vblock=0 inquire(file='QVIBR.OPT',exist=lex) if(lex)then write(6,*)'QVIBR.OPT found' open(9,file='QVIBR.OPT') 1 read(9,80,end=99,err=99)key 80 format(A4) if(key(1:3).eq.'E00')read(9,*)e00cm if(key(1:3).eq.'IWR')read(9,*)iwr if(key(1:3).eq.'LAB')read(9,*)lab if(key(1:3).eq.'LQ1')read(9,*)lq1 if(key(1:2).eq.'NP')read(9,*)np if(key.eq.'TOL1')read(9,*)tol1 if(key.eq.'VBLO')read(9,*)vblock if(key.eq.'TOL2')read(9,*)tol2 if(key.eq.'LGLG')read(9,*)lglg if(key.eq.'LMIN')read(9,*)lmin if(key.eq.'LMAX')read(9,*)lmax if(key.eq.'FWCM')read(9,*)fwcm backspace 9 read(9,80)k2 if(k2.eq.key)call report(k2//' unknown option') goto 1 99 close(9) endif e00=e00cm/CM dl=(lmax-lmin)/dble(np-1) return end subroutine maxdd(nb,b,NQ,d,r,dbig,rbig,mods,mors,mo,na,mmax) implicit none integer*4 nb,b,NQ,mods(nb,NQ),mors(nb,NQ),i,j,mmax, 1na(NQ),mo(NQ,mmax+1) real*8 d,r,dbig(nb),rbig(nb) do 1 i=1,nb if(d.gt.dbig(i))then do 2 j=nb,i+1,-1 2 dbig(j)=dbig(j-1) dbig(i)=d do 3 j=b+1,NQ 3 mods(i,j)=mo(j,na(j)) goto 11 endif 1 continue 11 continue do 4 i=1,nb if(dabs(r).gt.dabs(rbig(i)))then do 5 j=nb,i+1,-1 5 rbig(j)=rbig(j-1) rbig(i)=r do 6 j=b+1,NQ 6 mors(i,j)=mo(j,na(j)) return endif 4 continue return end subroutine bc(s,nb,vblock,NQ,dbig,mods,WP,CM,e00) implicit none integer*4 nb,NQ,mods(nb,NQ),vblock,j,i,ii,jj,ic real*8 dbig(nb),WP(NQ),e,CM,e00 character*1 s character*10 s10 write(6,600)s 600 format(' Biggest contributions ',a1,':') do 1 i=1,nb e=e00 do 21 j=vblock+1,NQ 21 e=e+WP(j)*dble(mods(i,j)) write(6,605)i,1.0d7/(e*CM),s,dbig(i) 605 format(i3,':',F12.2,' nm, ',A1,':',e12.4,' |0> -> |',$) ic=0 do 2 j=vblock+1,NQ if(mods(i,j).gt.0)then ic=ic+1 write(s10,100)j 100 format(i10) do 3 ii=1,len(s10) 3 if(s10(ii:ii).ne.' ')goto 4 4 do 5 jj=ii,len(s10) 5 write(6,602)s10(jj:jj) 602 format(A1,$) write(6,603) 603 format('^',$) write(s10,100)mods(i,j) do 6 ii=1,len(s10) 6 if(s10(ii:ii).ne.' ')goto 7 7 do 8 jj=ii,len(s10) 8 write(6,602)s10(jj:jj) endif 2 continue if(ic.eq.0) write(6,604) 604 format('0',$) 1 write(6,606) 606 format('>') return end