program polderdip implicit none integer*4 i,iu,n,nat,n3,iud,np,iw,n2,nst,ideg,MQ,nse real*8 alpha(3,3),wi,w,CM,d(3),g,A(3,3,3),GP(3,3), 1v(3),m(3),q(6),fi,qm(3,3),width,excnm,f1,f2,evmin, 1evmax,ev,dev,wau,alphav(3,3),cum(6),radapt, 1nmstart,nmend,dnm,evcm,nm,wdeg real*8,allocatable::etd(:,:),ald(:,:,:),gpd(:,:,:),ad(:,:,:,:), 1qd(:,:,:),x(:),alpha2(:),b2(:),aG(:),bA(:),bG(:),ge(:),gr(:), 1Aw(:,:,:,:),GPw(:,:,:),alphaw(:,:,:),alphavw(:,:,:), 1aldw(:,:,:,:),gpdw(:,:,:,:),adw(:,:,:,:,:), 1aldvw(:,:,:,:),aldv(:,:,:),aldj(:,:,:),S(:,:),E(:), 1als(:,:,:),gs(:,:,:),as(:,:,:,:),we(:),esc(:),dsc(:),rsc(:), 1aldQ(:,:,:,:),aldvQ(:,:,:,:),aldjQ(:,:,:,:), 1gpQ(:,:,:,:),adQ(:,:,:,:,:) integer*4,allocatable::lst(:) character*80 filename logical lint,lstat,lnm,lsep,ldog,usedd,usegr,lgr,lddeg,ladapt, 1escale,gswitch,lsum character*80 fgr,fgg call readopt(width,excnm,evmin,evmax,np,lint,lstat,lnm, 1lsep,ldog,CM,evcm,w,g,dev,nmstart,nmend,dnm,nat,usedd, 1usegr,fgr,lgr,lddeg,wdeg,fgg,ladapt,radapt,escale,gswitch, 1lsum) n3=3*nat n2=2*n3 allocate(b2(np),alpha2(np),aG(np),bA(np),bG(np), 1alphaw(np,3,3),GPw(np,3,3),Aw(np,3,3,3), 1alphavw(np,3,3),etd(16,n3+3),S(n3,n3),E(n3), 1qd(n3,3,3),x(n3),ge(n3),gr(n3),aldj(n2,3,3), 1ald(n2,3,3),aldv(n2,3,3),aldw(np,n2,3,3),aldvw(np,n2,3,3), 1gpd(n2,3,3),gpdw(np,n2,3,3),ad(n2,3,3,3),adw(np,n2,3,3,3)) c Polarizability derivatives,1..n3 real,n3+1..2n3 imaginary: c ald .. alpha c aldv .. velocity form c aldw,aldvw .. in frequency interval c gpd,dpdw .. G tensor c ad,adw .. A tensor call initt(np,n3,alpha,alphav,GP,A, 1b2,alpha2,aG,bA,bG,alphaw,GPw,Aw,alphavw,ald,aldv,aldj,gpd,ad, 1aldw,aldvw,gpdw,adw) c read ground-state gradient or set to zero: call rgsg(n3,gr,lgr,fgr) if(ladapt)then c read S-matrix call READS(N3,S,E,MQ) allocate(aldQ(MQ,n2,3,3),aldvQ(MQ,n2,3,3),aldjQ(MQ,n2,3,3), 1 gpQ(MQ,n2,3,3),adQ(MQ,n2,3,3,3)) aldQ=0.0d0 aldvQ=0.0d0 aldjQ=0.0d0 gpQ=0.0d0 adQ=0.0d0 endif c determine number of states: call dnst(nst,'TD.LST') allocate(lst(nst),we(nst)) if(escale)then c determine number of scale states: call dnst(nse,'ESC.LST') if(nse.ne.nst)call report('nse<>nst') allocate(esc(nst),dsc(nst),rsc(nst)) c read new wavelengths dipole and rotational strengths open(8,file='ESC.LST') do 2 i=1,nst 2 read(8,*)esc(i),dsc(i),rsc(i) close(8) endif c if lddeg read energies already here: call gen(nst,we,lddeg) iu=0 iud=0 open(9,file='TD.LST') do 1 i=1,nst read(9,90)filename 90 format(a80) c read derivatives of transition dipoles and gradient: call gd(i,filename,wi,d,n,v,m,q,n3,etd,x,ge,lst,lsum) c scale electronic variables: if(escale)call escala(nst,n3,i,wi,d,v,m,q,etd,esc,dsc,rsc) c redefine quadrupoles: call qqm(q,qm) call qqmd(etd,qd,n3) c polarizabilities: iu=iu+1 c add this state contribution to polarizabilities, for w,real part call addr(d,v,m,qm,alpha,alphav,GP,A,wi,w,g) c the same for frequency interval: if(lint)then ev=evmin-dev nm=nmstart+dnm do 101 iw=1,np call evset(lnm,ev,nm,dnm,dev) wau=ev*evcm/CM f1=2.0d0*fi(wi,wau,g) f2=-2.0d0*wi*fi(wi,wau,g) 101 call addw(iw,f1,f2,d,v,m,qm,np,alphaw,alphavw,GPw,Aw) endif if(lddeg)then ideg=0 do 12 iw=1,nst 12 if(dabs(we(iw)-wi).lt.wdeg/219474.5d0)ideg=ideg+1 if(ideg.gt.1)then write(6,612)ideg 612 format(i6,' x degeneracy detected, derivatives skipped') goto 1 endif endif c derivatives: iud=iud+1 c add this state contribution to polarizability derivatives call addd(d,v,m,qm,n3,ald,aldv,aldj,gpd,ad,etd,qd,lstat,wi,w, 1g,gr,ge,usegr,usedd) c record false traces of the tensors and derivatives: call wcum(cum,alpha,GP,A,ald,gpd,ad,n2,i) c the same for frequency interval: if(lint)call adddw(d,v,m,qm,np,n3,aldw,aldvw,gpdw, 1adw,etd,qd,lstat,wi,g,evmin,dev,lnm,ge,dnm,usegr) if(lsep) call drop(n3,i,wi,w,g,d,m,etd,qd,qm,ge,gswitch) if(np.ne.1)then write(6,*)np stop endif c derivatives for each mode if(ladapt)call addya(MQ,wi,w,g,E,qd,d,v,m,qm,n3,aldQ,aldvQ, 1aldjQ,gpQ,adQ,etd,radapt) 1 continue close(9) write(6,606)nst,iu,iud,width,excnm 606 format(i4,' states,',i4,' used,',i4,' used for derivatives',/, 1' width (cm-1): ',f12.2,' excnm (nm): ',f12.2,/) call wral(alpha,'POL.TTT',w) call wralv(alphav,'POLV.TTT',w) call wra(A,w) c G/w limit: GP=GP/w call wrg(GP,'GP.TTT',w) if(ldog)then write(6,*)'Distributed orihin gauge' call roadog(n3,x,ald,aldv,aldj,gpd,ad,w) if(ladapt)call dogq(n3,MQ,x,aldQ,aldvQ,aldjQ,gpQ,adQ,w) if(lint)call roadogw(np,n3,x,aldw,aldvw,gpdw,adw,w) endif if(lstat)then c add static limit of the polarizabilities allocate(als(N3,3,3),gs(N3,3,3),as(N3,3,3,3)) call READTTT(N3,als,as,gs) call adds1(N3,als,gs,as,ald,aldv,gpd,ad) if(lint)call addsw(N3,np,als,gs,as,aldw,aldvw,gpdw,adw) endif c G/w limit: gpd=gpd/w if(ladapt)then gpQ=gpq/w CALL TRTQ(S,MQ,N3,aldQ,aldvQ,gpQ,adQ,E) endif call WRITETTT('FILE.TTT',N3,ald,aldv,aldj,gpd,ad,w,gswitch) if(lint)then call mk(np,alpha2,aG,b2,bG,bA,alphaw,GPw,Aw) call wrw(alpha2,np,'alpha2.prn',evmin,dev,nmstart,dnm,lnm) call wrw(aG,np,'aG.prn',evmin,dev,nmstart,dnm,lnm) call wrw(b2,np,'b2.prn',evmin,dev,nmstart,dnm,lnm) call wrw(bG,np,'bG.prn',evmin,dev,nmstart,dnm,lnm) call wrw(bA,np,'bA.prn',evmin,dev,nmstart,dnm,lnm) call WRITETTTW(np,N3,aldw,aldvw,gpdw,adw,evmin,dev,nmstart, 1 dnm,lnm) endif call wrscale(np,evmin,dev,lnm,dnm,nmstart) end SUBROUTINE addsw(N3,np,als,gs,as,aldw,aldvw,gpdw,adw) IMPLICIT none INTEGER*4 n3,np,iw,ii,J,I,K real*8 als(N3,3,3),Gs(N3,3,3),as(N3,3,3,3), 1aldw(np,2*n3,3,3),gpdw(np,2*n3,3,3),adw(np,2*n3,3,3,3), 1aldvw(np,2*n3,3,3) do 1 iw=1,np do 1 ii=1,N3 do 1 i=1,3 do 1 j=1,3 aldw (iw,ii,i,j )=aldw (iw,ii,i,j )+als(ii,i,j ) aldvw(iw,ii,i,j )=aldvw(iw,ii,i,j )+als(ii,i,j ) gpdw (iw,ii,i,j )=gpdw (iw,ii,i,j )+gs (ii,i,j ) do 1 k=1,3 1 adw (iw,ii,i,j,k)=adw (iw,ii,i,j,k)+as (ii,i,j,k) return end SUBROUTINE adds1(N3,als,gs,as,ald,alv,GD,ad) IMPLICIT none INTEGER*4 n3,ii,J,I,K real*8 ald(2*N3,3,3),GD(2*N3,3,3),ad(2*N3,3,3,3) real*8 als(N3,3,3),Gs(N3,3,3),as(N3,3,3,3) real*8 alv(2*N3,3,3) do 1 ii=1,N3 do 1 i=1,3 do 1 j=1,3 ald(ii,i,j )=ald(ii,i,j )+als(ii,i,j ) alv(ii,i,j )=alv(ii,i,j )+als(ii,i,j ) gd (ii,i,j )=gd (ii,i,j )+gs (ii,i,j ) do 1 k=1,3 1 ad (ii,i,j,k)=ad (ii,i,j,k)+as (ii,i,j,k) return end SUBROUTINE READTTT(N,ALPHA,A,G) IMPLICIT none integer*4 N,I,NAT,L,IX,J,K,IIND real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3) NAT=N/3 open(15,file='STATIC.TTT',status='old') READ(15,*) READ(15,*) READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3) CLOSE(15) RETURN END C subroutine wrw(v,n,s,xi,dx,nmstart,dnm,lnm) implicit none real*8 v(*),x,xi,dx,nm,nmstart,dnm integer*4 n,i character*(*) s logical lnm open(12,file=s) x=xi-dx nm=nmstart+dnm do 1 i=1,n call evset(lnm,x,nm,dnm,dx) if(lnm)then write(12,120)nm,v(i) 120 format(f12.2,' ',g12.4) else write(12,121)x,v(i) 121 format(f12.6,' ',g12.4) endif 1 continue close(12) write(6,*)s//' written' return end subroutine wrscale(np,xi,dx,lnm,dnm,nmstart) implicit none real*8 xi,dx,ev,nm,evcm,dnm,nmstart integer*4 np,i logical lnm evcm=8065.73d0 open(12,file='xev.prn') open(13,file='xnm.prn') ev=xi-dx nm=nmstart+dnm do 1 i=1,np call evset(lnm,ev,nm,dnm,dx) nm=1.0d7/(ev*evcm) write(12,120)ev 1 write(13,121)nm 120 format(f12.4) 121 format(f12.2) close(12) close(13) write(6,*)'scales written' return end function rstate(filename) implicit none character*(*) filename integer*4 rstate,r,actual,ic character*(80) s80 r=0 actual=0 open(8,file=filename) 1 read(8,80,end=88,err=88)s80 80 format(a80) if(s80(2:14).eq.'Excited State')then do 2 ic=1,len(s80) 2 if(s80(ic:ic).eq.':')goto 3 3 read(s80(15:ic-1),*)actual endif if(s80(2:28).eq.'This state for optimization')then r=actual goto 88 endif goto 1 88 close(8) if(r.eq.0) 1call report(filename//': cannot determine excited state') rstate=r return end subroutine gd(io,filename,wi,d,r,v,m,q,n3,etd,x,ge,lst,lsum) c d transition dipole moment c v transition dipole moment, velocity c m transition magnetic dipole moment c q transition quadrupole moment, probably mixed c edt derivatives c ge gradient implicit none character*(*) filename character*(80) s80 integer*4 ie(8),i,ix,nd,r,n3,io,rstate,lst(*),j,k,ip real*8 wi,d(3),ev,v(3),m(3),q(6),etd(16,n3+3),x(n3),bohr,ge(*),gn, 1t(3),s(3) character*14 s8(8) data s8/ 'geometry ', 'dipoles ', 'v-dipoles ', 1 'm-dipoles ', 'quadrupoles ', 'Excited state ', 1 'derivatives ', 'gradient '/ integer ii(5) data ii/1,4,7,10,13/ character*10 nn(5) data nn/'r-dipole ','v-dipole ','m-dipole ','quadrupole', 1'quadrupole'/ character*1 xyz(3) data xyz/'X','Y','Z'/ logical lsum ie=0 c ie: 1:geometry c 2:dipoles c 3:v-dipoles c 4:m-dipoles c 5:quadrupoles c 6:Excited state c 7:derivatives c 8:gradient bohr=0.529177d0 v=0.0d0 m=0.0d0 q=0.0d0 etd=0.0d0 wi=0.0d0 c c determine the differentiated state in file io: r=rstate(filename) do 107 ix=len(filename),1,-1 107 if(filename(ix:ix).ne.' ')goto 1022 1022 write(6,*) do 102 j=1,ix 102 write(6,6052)filename(j:j) 6052 format(a1,$) write(6,6051)io,r 6051 format(': number',i6,', state',i6,' found') lst(io)=r do 101 I=1,io-1 if(lst(I).eq.r)then write(6,666)r,I 666 format(' State ',i6,' defined already in file #',i6) call report('state defined multiple times') endif 101 continue c need Input orientation here! open(8,file=filename) 1 read(8,80,end=88,err=88)s80 80 format(a80) c IF(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:')THEN DO 2004 I=1,4 2004 READ(8,*) DO 2005 I=1,n3/3 2005 READ(8,*)(x(ix+3*(I-1)),ix=1,3),(x(ix+3*(I-1)),ix=1,3) DO 2006 I=1,n3 2006 x(i)=x(i)/bohr ie(1)=ie(1)+1 ENDIF c Input orientation: if(s80(2:31).eq.'Electronic transition elements'.and.r.gt.0) 1call gete(r,8,d,v,m,q,ie) if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s'. 1 and.r.gt.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.r)then wi=ev/27.211384205943d0 ie(6)=ie(6)+1 endif endif c Input orientation: if(s80(2:34).eq.'Electronic Transition Derivatives'. 1 and.ie(7).eq.0)then call rmtrio(8,etd,16,n3+3) c probably: c lines: 1 energy c 2 - 4 el dipole c 5 - 7 el dipole gradient c 8 -10 mag dipole c 11-16 quadrupole c columns: 1-3 normal modes c 4-3nat+3 Cartesians ie(7)=ie(7)+1 endif c Input orientation: if(s80(38:59).eq.'Forces (Hartrees/Bohr)')then read(8,*) read(8,*) do 105 i=1,n3/3 105 read(8,*)(ge(ix+3*(i-1)),ix=1,2),(ge(ix+3*(i-1)),ix=1,3) gn=0.0d0 do 106 i=1,n3 gn=gn+ge(i)**2 106 ge(i)=-ge(i) gn=dsqrt(gn/dble(n3)) write(6,6001)ge(1),ge(n3),gn 6001 format(' G1 GN |G|:',3G12.4) ie(8)=ie(8)+1 endif goto 1 88 close(8) do 83 i=1,8 write(6,607)s8(i),ie(i) 83 if(mod(i,2).eq.0)write(6,*) 607 format(a14,i4,' x ',$) do 81 i=1,8 if(ie(i).eq.0)then if(i.ne.7)call report(' Incomplete output') endif 81 continue v=-v/wi q=-q/wi m= m/2.0d0 do 4 ix=1,3 do 4 i=1,n3 etd(4 +ix,3+i)=-etd(4 +ix,3+i)/wi etd(10+ix,3+i)=-etd(10+ix,3+i)/wi etd(13+ix,3+i)=-etd(13+ix,3+i)/wi 4 etd(7 +ix,3+i)= etd(7 +ix,3+i)/2.0d0 if(lsum)then do 31 ip=1,2 do 3 k =1,5 do 3 j =1,3 write(6,608)nn(k),xyz(j) 608 format(1X,A8,': sum rule for shift ',A1,':') s=0.0d0 do 2 ix=1,3 t(ix)= etd(ii(k)+ix,3+j+3*(1-1)) do 2 i=2,n3/3 2 s(ix)=s(ix)-etd(ii(k)+ix,3+j+3*(i-1)) 3 write(6,609)t,s 609 format(' First atom:',3f12.6,/, 1 ' Rest sum :',3f12.6) write(6,*) 31 if(ip.eq.1)call project(n3,x,etd) endif return end subroutine project(n3,x,etd) implicit none integer*4 n3,ia,ix,n6,i,j,k real*8 x(n3),etd(16,n3+3) real*8,allocatable::a(:,:),t(:),r(:,:) integer ii(5) data ii/1,4,7,10,13/ allocate(a(n3,n3),t(n3),r(3,n3/3)) n6=6 do 1 ia=1,n3/3 do 1 ix=1,3 1 r(ix,ia)=x(ix+3*(ia-1)) call DOMA(a,r,n3/3,n6) do 2 k=1,5 do 2 ix=1,3 t=0.0d0 do 3 i=1,n3 do 3 j=1,n3 3 t(i)=t(i)+a(i,j)*etd(ii(k)+ix,3+j) do 2 i=1,n3 2 etd(ii(k)+ix,3+i)=t(i) write(6,*)' Translation, rotation projected.' return end subroutine rmtrio(io,A,n,m) implicit none integer*4 n,N1,N3,io,LN,J,M real*8 A(n,m) 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 function fi(wi,w,g) c 2 2 c w - wi c ------------------ c 2 2 2 2 2 c (wi -w ) + wi g implicit none real*8 wi,w,g,fi,t t=w**2-wi**2 fi=t/(t**2+(wi*g)**2) return end function gi(wi,w,g) c c - wi g c ------------------ c 2 2 2 2 2 c (wi -w ) + wi g implicit none real*8 wi,w,g,gi,t t=wi**2-w**2 gi=-wi*g/(t**2+(wi*g)**2) return end function fg(wi,w,g) c 2 2 2 2 2 c wi g - (wi - w ) c ------------------------------ c 2 2 2 2 2 2 c [(wi -w ) + wi g ] implicit none real*8 wi,w,g,fg,di di=wi**2-w**2 fg=((wi*g)**2-di**2) / (di**2+(wi*g)**2)**2 return end function gg(wi,w,g) c 2 2 c (wi - w ) wi g c ------------------------------ c 2 2 2 2 2 2 c [(wi -w ) + wi g ] implicit none real*8 wi,w,g,gg,di di=wi**2-w**2 gg=di*g*wi / (di**2+(wi*g)**2)**2 return end subroutine getnat(fo,nat) implicit none character*(*) fo integer*4 nat,I,l,q character*80 s80 nat=0 open(2,file=fo) 1 read(2,2000,end=99,err=99)s80 2000 format(a80) IF(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:'.OR. 2 s80(20:40).EQ.'Standard orientation:'.OR. 2 s80(26:46).EQ.'Standard orientation:')THEN DO 2004 I=1,4 2004 READ(2,*) l=0 2005 READ(2,2000)s80 IF(s80(2:4).NE.'---')THEN l=l+1 BACKSPACE 2 READ(2,*)q,q IF(q.EQ.-1)l=l-1 GOTO 2005 ENDIF nat=l ENDIF goto 1 99 close(2) return end subroutine qqm(q,qm) implicit none real*8 q(6),qm(3,3) qm(1,1)=(3.0d0*q(1)-q(1)-q(2)-q(3))/2.0d0 qm(2,2)=(3.0d0*q(2)-q(1)-q(2)-q(3))/2.0d0 qm(3,3)=(3.0d0*q(3)-q(1)-q(2)-q(3))/2.0d0 qm(1,2)=q(4)*1.5d0 qm(2,1)=q(4)*1.5d0 qm(1,3)=q(5)*1.5d0 qm(3,1)=q(5)*1.5d0 qm(2,3)=q(6)*1.5d0 qm(3,2)=q(6)*1.5d0 return end subroutine qqmd(etd,qd,n3) implicit none real*8 qd(n3,3,3),etd(16,n3+3),q(6),qm(3,3) integer*4 j,ii,n3,ix,jx do 1 ii=1,n3 do 2 j=1,6 2 q(j)=etd(10+j,3+ii) call qqm(q,qm) do 1 ix=1,3 do 1 jx=1,3 1 qd(ii,ix,jx)=qm(ix,jx) return end SUBROUTINE WRITETTT(fn,N3,ald,alv,aldj,GD,ad,fr,gswitch) IMPLICIT none INTEGER*4 n3,L,IX,ii,J,I,K,NAT real*8 ald(2*N3,3,3),GD(2*N3,3,3),ad(2*N3,3,3,3),fr real*8 alv(2*N3,3,3),aldj(2*N3,3,3) character*(*) fn logical gswitch NAT=N3/3 OPEN(2,FILE=fn) WRITE(2,2000)NAT,fr 2000 FORMAT(' ROA tensors, cartesian derivatives',/, 1I4,' atoms, freq. ',f11.6/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 ii=3*(L-1)+IX 1 WRITE(2 ,2001)L,IX,(ald(ii,I,J),J=1,3),(ald(ii+n3,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,6g15.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') if(gswitch)then DO 6 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 6 L=1,NAT DO 6 IX=1,3 ii=3*(L-1)+IX 6 WRITE(2 ,2001)L,IX,(GD(n3+ii,I,J),J=1,3),(GD(ii,I,J),J=1,3) else DO 2 I=1,3 WRITE(2,2003)I DO 2 L=1,NAT DO 2 IX=1,3 ii=3*(L-1)+IX 2 WRITE(2 ,2001)L,IX,(GD(ii,I,J),J=1,3),(GD(n3+ii,I,J),J=1,3) endif WRITE(2,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 ii=3*(L-1)+IX 3 WRITE(2,2007)L,IX,(ad(ii,I,J,K),K=1,3), 1(ad(n3+ii,I,J,K),K=1,3),L,IX,I,J 2007 FORMAT(I5,1H ,I1,6g15.7,' ',4i3) write(2,*) write(2,*)'Alpha v:' DO 4 I=1,3 WRITE(2,2002)I DO 4 L=1,NAT DO 4 IX=1,3 ii=3*(L-1)+IX 4 WRITE(2 ,2001)L,IX,(alv(ii,I,J),J=1,3),(alv(ii+n3,I,J),J=1,3) write(2,*) write(2,*)'Alpha jv:' DO 5 I=1,3 WRITE(2,2002)I DO 5 L=1,NAT DO 5 IX=1,3 ii=3*(L-1)+IX 5 WRITE(2 ,2001)L,IX,(aldj(ii,I,J),J=1,3),(aldj(ii+n3,I,J),J=1,3) write(6,*)fn//' written' CLOSE(2) RETURN END SUBROUTINE WRITETTTW(np,N3,ald,aldv,gd,ad,evmin,dev,nmstart, 1dnm,lnm) IMPLICIT none INTEGER*4 np,n3,L,IX,ii,J,I,K,NAT,iw real*8 ald(np,2*N3,3,3),gd(np,2*N3,3,3),ad(np,2*N3,3,3,3), 1fr,evmin, 1dev,ev,aldv(np,2*N3,3,3),evcm,nmstart,dnm,nm logical lnm evcm=8065.73d0 NAT=N3/3 OPEN(2,FILE='FILEW.TTT') ev=evmin-dev nm=nmstart+dnm do 4 iw=1,np call evset(lnm,ev,nm,dnm,dev) fr=ev*evcm WRITE(2,2000)NAT,fr,ev,1.0d7/fr 2000 FORMAT(' ROA tensors, cartesian derivatives',/, 1I4,' atoms, freq. ',f11.2,' cm-1',f11.5,' ev,',f10.2,' nm'/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 ii=3*(L-1)+IX 1 WRITE(2 ,2001)L,IX,(ald(iw,ii,I,J),J=1,3), 1(ald(iw,ii+n3,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,6g15.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz), v-form:') DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 ii=3*(L-1)+IX 2 WRITE(2 ,2001)L,IX,(GD(iw,ii,I,J),J=1,3),(GD(iw,n3+ii,I,J),J=1,3) WRITE(2,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 ii=3*(L-1)+IX 3 WRITE(2,2007)L,IX,(ad(iw,ii,I,J,K),K=1,3), 1(ad(iw,ii+n3,I,J,K),K=1,3),L,IX,I,J 2007 FORMAT(I5,1H ,I1,6g15.7,' ',4i3) write(2,*) write(2,*)' Alpha v:' DO 4 I=1,3 WRITE(2,2002)I DO 4 L=1,NAT DO 4 IX=1,3 ii=3*(L-1)+IX 4 WRITE(2 ,2001)L,IX,(aldv(iw,ii,I,J),J=1,3), 1(aldv(iw,ii+n3,I,J),J=1,3) CLOSE(2) RETURN END subroutine addr(d,v,m,qm,alpha,alphav,GP,A,wi,w,g) implicit none integer*4 i,j real*8 d(3),v(3),m(3),qm(3,3),GP(3,3),alpha(3,3),A(3,3,3), 1alphav(3,3),wi,w,g,ff,gg,f2,f1 ff=(wi-w)/((wi-w)**2+g**2) gg=(wi+w)/((wi+w)**2+g**2) f2=ff+gg f1=ff-gg do 2 i=1,3 do 2 j=1,3 alpha( i,j)=alpha( i,j)+f2*d(i)*d(j) alphav( i,j)=alphav( i,j)+f1*d(i)*v(j)*wi/w GP( i,j)=GP( i,j)+f1*d(i)*m(j) A(i,j,1) =A(i,j,1)+f2*d(i)*qm(j,1) A(i,j,2) =A(i,j,2)+f2*d(i)*qm(j,2) 2 A(i,j,3) =A(i,j,3)+f2*d(i)*qm(j,3) return end subroutine addw(iw,f1,f2,d,v,m,qm,np,alphaw,alphavw,GPw,Aw) implicit none integer*4 np,ix,jx,kx,iw real*8 f1,f2,d(*),m(*),qm(3,3),alphaw(np,3,3),alphavw(np,3,3), 1GPw(np,3,3),Aw(np,3,3,3),v(*) do 101 ix=1,3 do 101 jx=1,3 GPw( iw,ix,jx) =GPw( iw,ix,jx)+f1*d(ix)*m(jx) alphaw( iw,ix,jx)=alphaw( iw,ix,jx)+f2*d(ix)*d(jx) alphavw(iw,ix,jx)=alphavw(iw,ix,jx)+f2*d(ix)*v(jx) do 101 kx=1,3 101 Aw(iw,ix,jx,kx) =Aw(iw,ix,jx,kx )+f2*d(ix)*qm(jx,kx) return end subroutine addx(f1,f2,f3,f4,g1,g2,g3,g4, 1d,v,m,qm,n3,ald,aldv,gpd,ad,etd,qd,ge,usegr) implicit none integer*4 ii,n3,ix,jx,kx real*8 f1,f2,f3,g1,g2,g3,g4, 1f4,d(*),v(*),m(*),qm(3,3),ux,uxd,etd(16,n3+3), 1ald(2*n3,3,3),aldv(2*n3,3,3), 1gpd(2*n3,3,3),ad(2*n3,3,3,3),qd(n3,3,3),uy,uyd,gd, 1ge(*),vy,vyd,my,myd,y3,x3,x4,y4 logical usegr if(usegr)then y3=f3 x3=g3 y4=f4 x4=g4 else y3=0.0d0 x3=0.0d0 y4=0.0d0 x4=0.0d0 endif do 3 ii=1,n3 gd=ge(ii) do 3 ix=1,3 ux =d (ix) uxd=etd(ix+1,ii+3) do 3 jx=1,3 uy =d (jx) vy =v (jx) my =m (jx) uyd=etd(jx+1,ii+3) vyd=etd(jx+4,ii+3) myd=etd(jx+7,ii+3) ald( ii,ix,jx)=ald( ii,ix,jx)+f2*(uxd*uy+ux*uyd)+y3*ux*uy*gd aldv(ii,ix,jx)=aldv(ii,ix,jx)+f2*(uxd*vy+ux*vyd)+y3*ux*vy*gd ald( n3+ii,ix,jx)=ald (n3+ii,ix,jx)+g2*(uxd*uy+ux*uyd)+x3*ux*uy*gd aldv(n3+ii,ix,jx)=aldv(n3+ii,ix,jx)+g2*(uxd*vy+ux*vyd)+x3*ux*vy*gd c gpd: last index - magnetic gpd( ii,ix,jx)=gpd( ii,ix,jx)+f1*(uxd*my+ux*myd)+y4*ux*my*gd gpd(n3+ii,ix,jx)=gpd(n3+ii,ix,jx)+g1*(uxd*my+ux*myd)+x4*ux*my*gd do 3 kx=1,3 c ad - first index - dipole ad( ii,ix,jx,kx)=ad( ii,ix,jx,kx)+f2* 1(uxd*qm(jx,kx)+ux*qd(ii,jx,kx))+y3*ux*qm(jx,kx)*gd 3 ad(n3+ii,ix,jx,kx)=ad(n3+ii,ix,jx,kx)+g2* 1(uxd*qm(jx,kx)+ux*qd(ii,jx,kx))+x3*ux*qm(jx,kx)*gd return end subroutine addy(wi,w,g, 1d,v,m,qm,n3,ald,aldv,aldj,gp,ad,etd,qd,gr,ge,usegr,usedd) implicit none integer*4 ii,n3,ix,jx,X real*8 wi,w,g,d(*),v(*),m(*),qm(3,3),ux,uxd,etd(16,n3+3), 1ald(2*n3,3,3),aldv(2*n3,3,3),fv,fw,aldj(2*n3,3,3), 1gp(2*n3,3,3),ad(2*n3,3,3,3),qd(n3,3,3),uy,uyd,gd,fsi,ff,fdi,d2i, 1ge(*),vy,vyd,my,myd,gg,dg,fs2,d2,fs,fd,u,p,fi,fd2,fm,dgi,gr(*) logical usegr,usedd u=(wi-w)**2+g**2 p=(wi+w)**2+g**2 ff=(wi-w)/u gg=(wi+w)/p fi=-g/u-g/p fm=-g/u+g/p fs=ff+gg fs2=ff**2+gg**2 fd2=ff**2-gg**2 fd=ff-gg fsi=2.0d0*g*((wi-w)/u**2+(wi+w)/p**2) fdi=2.0d0*g*((wi-w)/u**2-(wi+w)/p**2) fv=fd*wi/w fw=fm*wi/w do 1 ii=1,n3 c transition dipole derivative part: if(usedd)then do 2 ix=1,3 ux =d (ix) uxd=etd(ix+1,ii+3) do 2 jx=1,3 uy =d (jx) vy =v (jx) my =m (jx) uyd=etd(jx+1,ii+3) vyd=etd(jx+4,ii+3) myd=etd(jx+7,ii+3) c real parts of frequency functions: ald( ii,ix,jx)=ald( ii,ix,jx)+(uxd*uy+ux*uyd)*fs c aldv(ii,ix,jx)=aldv(ii,ix,jx)+(uxd*vy+ux*vyd)*fs aldv(ii,ix,jx)=aldv(ii,ix,jx)+(uxd*vy+ux*vyd)*fv aldj(ii,ix,jx)=aldj(ii,ix,jx)+(uxd*vy+ux*vyd)*fs c aldj(ii,ix,jx)=aldj(ii,ix,jx)+(uxd*vy+ux*vyd)*fv c gp: last index - magnetic gp(ii,ix,jx)=gp(ii,ix,jx)+(uxd*my+ux*myd)*fd c ad - first index - dipole ad(ii,ix,jx,1)=ad(ii,ix,jx,1)+(uxd*qm(jx,1)+ux*qd(ii,jx,1))*fs ad(ii,ix,jx,2)=ad(ii,ix,jx,2)+(uxd*qm(jx,2)+ux*qd(ii,jx,2))*fs ad(ii,ix,jx,3)=ad(ii,ix,jx,3)+(uxd*qm(jx,3)+ux*qd(ii,jx,3))*fs c complex: X=ii+n3 ald( X ,ix,jx)=ald( X ,ix,jx)+(uxd*uy+ux*uyd)*fi c aldv(X ,ix,jx)=aldv(X ,ix,jx)+(uxd*vy+ux*vyd)*fi aldv(X ,ix,jx)=aldv(X ,ix,jx)+(uxd*vy+ux*vyd)*fw c aldj(X ,ix,jx)=aldj(X ,ix,jx)+(uxd*vy+ux*vyd)*fw aldj(X ,ix,jx)=aldj(X ,ix,jx)+(uxd*vy+ux*vyd)*fi gp(X,ix,jx)=gp(X,ix,jx)+(uxd*my+ux*myd)*fm ad(X,ix,jx,1)=ad(X,ix,jx,1)+(uxd*qm(jx,1)+ux*qd(ii,jx,1))*fi ad(X,ix,jx,2)=ad(X,ix,jx,2)+(uxd*qm(jx,2)+ux*qd(ii,jx,2))*fi 2 ad(X,ix,jx,3)=ad(X,ix,jx,3)+(uxd*qm(jx,3)+ux*qd(ii,jx,3))*fi endif if(usegr)then c ground state gradient: gg=gr(ii) c excited state gradient: gd=ge(ii) c difference of gradients(1/(wi-w)^2+1/(wi+w)^`2): dg =(gd-gg)*fs2 d2 =(gd-gg)*fd2 dgi=(gd-gg)*fsi d2i=(gd-gg)*fdi c d3r=(ff**2+gg**2+(ff+gg)/wi)*(gd-gg) c d3i=-(2.0d0*(ff/u+gg/p)+(1.0d0/u**2+1.0d0/p**2)/wi)*g*(gd-gg) do 3 ix=1,3 ux =d (ix) do 3 jx=1,3 uy =d (jx) vy =v (jx) my =m (jx) ald( ii,ix,jx)=ald( ii,ix,jx)-ux*uy*dg c aldv(ii,ix,jx)=aldv(ii,ix,jx)-ux*vy*dg aldv(ii,ix,jx)=aldv(ii,ix,jx)-ux*vy*d2*wi/w c aldj(ii,ix,jx)=aldj(ii,ix,jx)-ux*vy*d3r aldj(ii,ix,jx)=aldj(ii,ix,jx)-ux*vy*dg gp(ii,ix,jx)=gp(ii,ix,jx)-ux*my*d2 ad(ii,ix,jx,1)=ad(ii,ix,jx,1)-ux*qm(jx,1)*dg ad(ii,ix,jx,2)=ad(ii,ix,jx,2)-ux*qm(jx,2)*dg ad(ii,ix,jx,3)=ad(ii,ix,jx,3)-ux*qm(jx,3)*dg X=ii+n3 ald( X,ix,jx)=ald( X,ix,jx)-ux*uy*dgi c aldv(X,ix,jx)=aldv(X,ix,jx)-ux*vy*dgi aldv(X,ix,jx)=aldv(X,ix,jx)-ux*vy*d2i*wi/w c aldj(X,ix,jx)=aldj(X,ix,jx)-ux*vy*d2i*d3i aldj(X,ix,jx)=aldj(X,ix,jx)-ux*vy*dgi gp(X,ix,jx)=gp(X,ix,jx)-ux*my*d2i ad(X,ix,jx,1)=ad(X,ix,jx,1)-ux*qm(jx,1)*dgi ad(X,ix,jx,2)=ad(X,ix,jx,2)-ux*qm(jx,2)*dgi 3 ad(X,ix,jx,3)=ad(X,ix,jx,3)-ux*qm(jx,3)*dgi endif 1 continue return end subroutine setf(ic,f1,f2,f3,f4,g1,g2,g3,g4,wi,w,g) implicit none integer*4 ic real*8 f1,f2,f3,f4,g1,g2,g3,g4,wi,w,g,fi,fg,gi,gg if(ic.eq.1)then f1= 2.0d0*fi(wi,w,g) f2=-2.0d0*wi*fi(wi,w,g) f3= 2.0d0*(wi**2+w**2)*fg(wi,w,g) f4=-4.0d0*fg(wi,w,g) g1=-2.0d0*gi(wi,w,g) g2=2.0d0*wi*gi(wi,w,g) g3=4.0d0*(wi**2+w**2)*gg(wi,w,g) g4=-8.0d0*wi*gg(wi,w,g) else f1=-2.0d0*fi(wi,0.0d0,0.0d0) f2= 2.0d0*wi*fi(wi,0.0d0,0.0d0) f3=-2.0d0*(wi**2+w**2)*fg(wi,0.0d0,0.0d0) f4= 4.0d0*fg(wi,0.0d0,0.0d0) g1= 2.0d0*gi(wi,0.0d0,0.0d0) g2=-2.0d0*wi*gi(wi,0.0d0,0.0d0) g3=-4.0d0*(wi**2+w**2)*gg(wi,0.0d0,0.0d0) g4= 8.0d0*wi*gg(wi,0.0d0,0.0d0) endif return end c ============================================================ subroutine addd(d,v,m,qm,n3,ald,aldv,aldj,gpd,ad,etd,qd,lstat, 1wi,w,g,gr,ge,usegr,usedd) implicit none integer*4 n3 real*8 f1,f2,d(3),v(3),m(3),qm(3,3),etd(16,n3+3), 1ald(2*n3,3,3),gpd(2*n3,3,3),ad(2*n3,3,3,3),qd(n3,3,3), 1aldv(2*n3,3,3),wi,w,g,ge(n3),f3,f4,gr(n3),aldj(2*n3,3,3), 1g1,g2,g3,g4 logical lstat,usegr,usedd call setf(1,f1,f2,f3,f4,g1,g2,g3,g4,wi,w,g) call addy(wi,w,g, 1d,v,m,qm,n3,ald,aldv,aldj,gpd,ad,etd,qd,gr,ge,usegr,usedd) if(lstat)then c delete the contribution of this transition to static part call setf(0,f1,f2,f3,f4,g1,g2,g3,g4,wi,w,g) call addx(f1,f2,f3,f4,g1,g2,g3,g4, 1 d,v,m,qm,n3,ald,aldv,gpd,ad,etd,qd,ge,usegr) endif return end subroutine axx(d,v,m,qm,np,n3,ald,alv,gpd, 1ad,etd,qd,g,evmin,dev,wi,ic,lnm,ge,dnm,usegr) implicit none integer*4 np,n3,iw,ic logical lnm,usegr real*8 f1,f2,d(*),v(*),m(*),qm(3,3),etd(16,n3+3), 1ald(np,2*n3,3,3),gpd(np,2*n3,3,3),ad(np,2*n3,3,3,3), 1qd(n3,3,3),wi,alv(np,2*n3,3,3),ev,wau,CM,evmin,dev,g,evcm, 1nm,nmstart,dnm,ge(*),f3,f4,g1,g2,g3,g4 real*8,allocatable::al1(:,:,:),gg1(:,:,:),a1(:,:,:,:),av1(:,:,:) allocate(al1(2*n3,3,3),gg1(2*n3,3,3),a1(2*n3,3,3,3),av1(2*n3,3,3)) CM=219474.0d0 evcm=8065.73d0 ev=evmin-dev nmstart=1.0d7/(evcm*evmin) nm=nmstart+dnm do 1 iw=1,np call evset(lnm,ev,nm,dnm,dev) wau=ev*evcm/CM call setf(ic,f1,f2,f3,f4,g1,g2,g3,g4,wi,wau,g) call tofrom(iw,0,np,n3,ald,alv,gpd,ad,al1,av1,gg1,a1) call addx(f1,f2,f3,f4,g1,g2,g3,g4, 1d,v,m,qm,n3,al1,av1,gg1,a1,etd,qd,ge,usegr) 1 call tofrom(iw,1,np,n3,ald,alv,gpd,ad,al1,av1,gg1,a1) return end subroutine adddw(d,v,m,qm,np,n3,ald,alv,gpd, 1ad,etd,qd,lstat,wi,g,evmin,dev,lnm,ge,dnm,usegr) implicit none integer*4 np,n3 real*8 d(*),v(*),m(*),qm(3,3),etd(16,n3+3),ge(*), 1ald(np,2*n3,3,3),gpd(np,2*n3,3,3),ad(np,2*n3,3,3,3),qd(n3,3,3), 1alv(np,2*n3,3,3),wi,g,evmin,dev,dnm logical lstat,lnm,usegr call axx(d,v,m,qm,np,n3,ald,alv,gpd,ad,etd,qd,g, 1evmin,dev,wi,1,lnm,ge,dnm,usegr) if(lstat)call axx(d,v,m,qm,np,n3,ald,alv,gpd,ad,etd,qd,g, 1evmin,dev,wi,2,lnm,ge,dnm,usegr) return end subroutine readopt(width,excnm,evmin,evmax,np,lint,lstat, 1lnm,lsep,ldog,CM,evcm,w,g,dev,nmstart,nmend,dnm,nat,usedd, 1usegr,fgr,lgr,lddeg,wdeg,fgg,ladapt,radapt,escale,gswitch, 1lsum) implicit none integer*4 np,nat real*8 width,excnm,evmin,evmax,CM,evcm,w,g,dev,nmstart,nmend, 1dnm,wdeg,radapt logical lint,lex,lstat,lnm,lsep,ldog,usedd,usegr,lgr,lddeg,ladapt, 1escale,gswitch,lsum character*80 s80,filename,fgr,fgg,st CM=219474.0d0 evcm=8065.73d0 c scale electronic enegries, dipole and rotational strengths escale=.false. c adapt for every fundamental transition: ladapt=.false. c adaptation multiplier: radapt=1.0d0 c CP FILE.TTT: fgg='CP.TTT' c electronic band width in cm-1 width=0.0d0 c excitation frequency in nm excnm=532.0d0 c excitation freq limit for scan, in eV evmin=4.0d0 c excitation freq limit for scan, in eV evmax=9.0d0 c number of exc frequencies on the scan np=500 c scan an interval of exc. frequencies: lint=.false. c distributed origin gauge for G: ldog=.true. c read static derivatives: lstat=.false. c scan linear in nm: lnm=.false. c separate ttt files for each transition lsep=.false. c use transition derivatives usedd=.true. c use gradient usegr=.true. c file with ground state gradient fgr='GROUND.OUT' c use ground state gradient: lgr=.false. c ignore degenerate transitions: lddeg=.false. c the limit (cm-1) for degenerate transitions: wdeg=100.0d0 c put G' in FILE.TTT second gswitch=.false. c derivative sum rules check lsum=.false. inquire(file='POLDERDIP.OPT',exist=lex) if(lex)then open(9,file='POLDERDIP.OPT') 1111 read(9,90,end=991,err=991)s80 90 format(a80) if(s80(1:5).eq.'ADAPT')read(9,*)ladapt if(s80(1:5).eq.'RADAP')read(9,*)radapt if(s80(1:5).eq.'ESCAL')read(9,*)escale if(s80(1:5).eq.'WIDTH')read(9,*)width if(s80(1:5).eq.'EXCNM')read(9,*)excnm if(s80(1:5).eq.'GSWIT')read(9,*)gswitch if(s80(1:5).eq.'INTER')read(9,*)lint if(s80(1:5).eq.'EVMIN')read(9,*)evmin if(s80(1:5).eq.'EVMAX')read(9,*)evmax if(s80(1:5).eq.'NPOIN')read(9,*)np if(s80(1:4).eq.'LDOG' )read(9,*)ldog if(s80(1:5).eq.'FILEG')read(9,*)fgg if(s80(1:5).eq.'STATI')read(9,*)lstat if(s80(1:5).eq.'NMLIN')read(9,*)lnm if(s80(1:5).eq.'SEPAR')read(9,*)lsep if(s80(1:5).eq.'USEDD')read(9,*)usedd if(s80(1:5).eq.'SUMRU')read(9,*)lsum if(s80(1:5).eq.'USEGR')read(9,*)usegr if(s80(1:5).eq.'GROUN')read(9,*)fgr if(s80(1:5).eq.'LGROU')read(9,*)lgr if(s80(1:5).eq.'LDDEG')read(9,*)lddeg if(s80(1:3).eq.'END')goto 991 backspace 9 read(9,90)st if(st.eq.s80)then write(6,*)st call report('Unknown option') endif goto 1111 991 close(9) endif w=1.0d7/excnm/CM g=width/CM nmstart=1.0d7/(evcm*evmin) nmend =1.0d7/(evcm*evmax) if(np.gt.1)then dev=(evmax -evmin)/dble(np-1) dnm=(nmstart-nmend)/dble(np-1) else dev=0.0d0 dnm=0.0d0 endif if(.not.lint)np=1 nat=0 open(9,file='TD.LST') read(9,90,end=99,err=99)filename call getnat(filename,nat) 99 close(9) if(nat.eq.0)then call report('no atoms') else write(6,*)nat,' atoms' endif write(6,600)excnm 600 format(' wexc:',f10.2,' nm') return end subroutine mk(np,alpha2,aG,b2,bG,bA,alw,GPw,Aw) implicit none integer*4 np,iw,i,j,ig,id real*8 alpha2(np),aG(*),b2(np),bG(np),bA(np),spal,spG, 1alw(np,3,3),GPw(np,3,3),Aw(np,3,3,3) b2=0.0d0 bA=0.0d0 bG=0.0d0 do 1 iw=1,np spal=(alw(iw,1,1)+alw(iw,2,2)+alw(iw,3,3))/3.0d0 spG =(GPw(iw,1,1)+GPw(iw,2,2)+GPw(iw,3,3))/3.0d0 do 103 i=1,3 do 103 j=1,3 b2(iw)=b2(iw)+1.5d0*alw(iw,i,j)*alw(iw,i,j) bG(iw)=bG(iw)+1.5d0*alw(iw,i,j)*GPw(iw,i,j) ig=i+1 if(ig.gt.3)ig=1 id=ig+1 if(id.gt.3)id=1 103 bA(iw)=bA(iw)+0.5d0*alw(iw,i,j)*(Aw(iw,ig,id,j)-Aw(iw,id,ig,j)) alpha2(iw)=spal*spal aG( iw)=spal*spG b2(iw)=b2(iw)-4.5d0*alpha2(iw) 1 bG(iw)=bG(iw)-4.5d0*aG(iw) return end subroutine evset(lnm,ev,nm,dnm,dev) implicit none logical lnm real*8 ev,nm,dnm,dev,evcm evcm=8065.73d0 if(lnm)then nm=nm-dnm ev=(1.0d7/nm)/evcm else ev=ev+dev endif return end subroutine drop(n3,io,wi,w,gg,d,m,etd,qd,qm,ge,gswitch) implicit none integer*4 i,j,k,io,is,n3,ii real*8 wi,w,gg,f1,f2,ux,uxd,d(*),etd(16,n3+3), 1m(*),qd(n3,3,3),qm(3,3),ge(*),f3,f4,gd,g1,g2,g3,g4 real*8,allocatable::al(:,:,:),g(:,:,:),a(:,:,:,:) character*10 s10 logical gswitch write(s10,200)io 200 format(i10) do 1 is=1,len(s10) 1 if(s10(is:is).ne.' ')goto 2 2 allocate(al(2*n3,3,3),g(2*n3,3,3),a(2*n3,3,3,3)) al=0.0d0 g=0.0d0 a=0.0d0 call setf(1,f1,f2,f3,f4,g1,g2,g3,g4,wi,w,gg) do 3 ii=1,n3 gd=ge(ii) do 3 i=1,3 ux =d( i) uxd=etd(i+1,ii+3) do 3 j=1,3 al(ii ,i,j) =al(ii ,i,j)+f2*(uxd*d(j)+ux*etd(j+1,ii+3)) 1-f3*ux*d(j)*gd al(ii+n3,i,j) =al(ii+n3,i,j)+g2*(uxd*d(j)+ux*etd(j+1,ii+3)) 1-g3*ux*d(j)*gd g( ii ,i,j) =g( ii ,i,j)+f1*(uxd*m(j)+ux*etd(j+7,ii+3)) 1-f4*ux*m(j)*gd g( ii+n3,i,j) =g( ii+n3,i,j)+g1*(uxd*m(j)+ux*etd(j+7,ii+3)) 1-g4*ux*m(j)*gd do 3 k=1,3 a( ii ,i,j,k)=a(ii ,i,j,k)+f2*(uxd*qm(j,k)+ux*qd(ii,j,k)) 1-f3*ux*qm(j,k)*gd 3 a( ii+n3,i,j,k)=a(ii+n3,i,j,k)+g2*(uxd*qm(j,k)+ux*qd(ii,j,k)) 1-g3*ux*qm(j,k)*gd call WRITETTT(s10(is:len(s10))//'.ttt',n3,al,al,al,g,a,w,gswitch) return end subroutine report(r) character*(*) r write(6,*)r stop end subroutine wral(a,f,w) implicit none real*8 a(3,3),w character*(*) f integer*4 ix,jx character*1 x(3) data x/'X','Y','Z'/ open(9,file=f) write(9,900)' Polarizability',w 900 format(a15,' w = ',F12.6) do 901 ix=1,3 do 901 jx=1,ix 901 write(9,902)x(ix),x(jx) 902 format(8x,2a1,$) write(9,*) c XX XY YY XZ YZ ZZ write(9,600)((a(ix,jx),jx=1,ix),ix=1,3) 600 format(6f10.4) close(9) return end subroutine wralv(a,f,w) implicit none real*8 a(3,3),w character*(*) f integer*4 ix,jx character*1 x(3) data x/'X','Y','Z'/ open(9,file=f) write(9,901)' Polarizability',w 901 format(a15,' w = ',F12.6,' length-velocity, last index v') do 900 ix=1,3 do 900 jx=1,3 900 write(9,902)x(ix),x(jx) 902 format(8x,2a1,$) write(9,*) c XX XY XZ YX YY YZ ZX ZY ZZ write(9,600)((a(ix,jx),jx=1,3),ix=1,3) 600 format(9f10.4) close(9) return end subroutine wra(a,w) implicit none real*8 a(3,3,3),w integer*4 ix,jx,kx character*1 x(3) data x/'X','Y','Z'/ open(9,file='A.TTT') write(9,900)' A, dipole-quad',w 900 format(a15,' w = ',F12.6) do 901 ix=1,3 do 901 jx=1,ix 901 write(9,902)x(ix),x(jx) 902 format(8x,2a1,$) write(9,*) c XX XY YY XZ YZ ZZ write(9,600)(((a(kx,ix,jx),jx=1,ix),ix=1,3),kx=1,3) 600 format(6f10.4) close(9) return end subroutine wrg(g,s,w) implicit none real*8 g(3,3),w integer*4 ix,jx character*(*) s character*1 x(3) data x/'X','Y','Z'/ open(9,file=s) write(9,901)' G-tensor ',w 901 format(a15,' w = ',F12.6) do 900 ix=1,3 do 900 jx=1,3 900 write(9,902)x(ix),x(jx) 902 format(8x,2a1,$) write(9,*) c XX XY XZ YX YY YZ ZX ZY ZZ write(9,601)((g(ix,jx),jx=1,3),ix=1,3) 601 format(9f10.4) close(9) return end subroutine dogq(n3,M,x,aldQ,alvQ,aljQ,GQ,QQ,w) implicit none integer*4 n3,a,b,c,i,M,iq real*8 x(*),aldQ(M,2*n3,3,3),alvQ(M,2*n3,3,3),GQ(M,2*n3,3,3), 1QQ(M,2*n3,3,3,3),w,aljQ(M,2*n3,3,3), 1ald(2*n3,3,3),alv(2*n3,3,3),G(2*n3,3,3), 1Q(2*n3,3,3,3),alj(2*n3,3,3) do 1 iq=1,M do 2 i=1,n3+n3 do 2 a=1,3 do 2 b=1,3 ald(i,a,b)=aldQ(iq,i,a,b) alv(i,a,b)=alvQ(iq,i,a,b) alj(i,a,b)=aljQ(iq,i,a,b) G (i,a,b)=GQ (iq,i,a,b) do 2 c=1,3 2 Q (i,a,b,c)=QQ (iq,i,a,b,c) call roadog(n3,x,ald,alv,alj,G,Q,w) do 1 i=1,n3+n3 do 1 a=1,3 do 1 b=1,3 aldQ(iq,i,a,b)=ald(i,a,b) alvQ(iq,i,a,b)=alv(i,a,b) aljQ(iq,i,a,b)=alj(i,a,b) GQ (iq,i,a,b)=G (i,a,b) do 1 c=1,3 1 QQ (iq,i,a,b,c)=Q (i,a,b,c) return end subroutine roadog(n3,x,ald,aldv,aldj,G,Q,w) implicit none integer*4 n3,a,b,c,L,IE,i,ir real*8 x(*),ald(2*n3,3,3),aldv(2*n3,3,3),G(2*n3,3,3), 1Q(2*n3,3,3,3),d(3),R(3),w,s,aldj(2*n3,3,3) DO 1 L=0,n3/3-1 R(1)=X(1+3*L) R(2)=X(2+3*L) R(3)=X(3+3*L) DO 1 IE=1,3 c ir: real and imaginary parts: do 1 ir=0,1 i=3*L+IE+n3*ir DO 1 a=1,3 d(1)=aldv(i,a,1)-ald(i,a,1) d(2)=aldv(i,a,2)-ald(i,a,2) d(3)=aldv(i,a,3)-ald(i,a,3) c G,Q: first index electric dipole G(i,a,1)=G(i,a,1)+0.5d0*w*(R(2)*d(3)-R(3)*d(2)) G(i,a,2)=G(i,a,2)+0.5d0*w*(R(3)*d(1)-R(1)*d(3)) G(i,a,3)=G(i,a,3)+0.5d0*w*(R(1)*d(2)-R(2)*d(1)) d(1)=aldj(i,a,1)-ald(i,a,1) d(2)=aldj(i,a,2)-ald(i,a,2) d(3)=aldj(i,a,3)-ald(i,a,3) s=R(1)*d(1)+R(2)*d(2)+R(3)*d(3) DO 1 b=1,3 Q(i,a,b,b)=Q(i,a,b,b)+s DO 1 c=1,3 1 Q(i,a,b,c)=Q(i,a,b,c)-1.5d0*(R(b)*d(c)+R(c)*d(b)) return end subroutine roadogw(np,n3,x,aldw,aldvw,Gw,Qw,w) implicit none integer*4 np,n3,j real*8 x(*),aldw(np,2*n3,3,3),aldvw(np,2*n3,3,3),Gw(np,2*n3,3,3), 1Qw(np,2*n3,3,3,3),ald(2*n3,3,3),aldv(2*n3,3,3),G(2*n3,3,3), 1Q(2*n3,3,3,3),w real*8,allocatable::aldj(:,:,:) allocate(aldj(2*n3,3,3)) aldj=0.0d0 DO 1 j=1,np call tofrom(j,0,np,n3,aldw,aldvw,Gw,Qw,ald,aldv,G,Q) call roadog(n3,x,ald,aldv,aldj,G,Q,w) 1 call tofrom(j,1,np,n3,aldw,aldvw,Gw,Qw,ald,aldv,G,Q) return end subroutine tofrom(j,ic,np,n3,aldw,aldvw,Gw,Qw,ald,aldv,G,Q) implicit none integer*4 ic,np,n3,j,l,a,b,c real*8 aldw(np,2*n3,3,3),aldvw(np,2*n3,3,3),Gw(np,2*n3,3,3), 1Qw(np,2*n3,3,3,3),ald(2*n3,3,3),aldv(2*n3,3,3),G(2*n3,3,3), 1Q(2*n3,3,3,3) if(ic.eq.0)then do 1 l=1,2*n3 do 1 a=1,3 do 1 b=1,3 ald( l,a,b)=aldw( j,l,a,b) aldv(l,a,b)=aldvw(j,l,a,b) G( l,a,b)=Gw( j,l,a,b) do 1 c=1,3 1 Q( l,a,b,c)=Qw( j,l,a,b,c) else do 2 l=1,2*n3 do 2 a=1,3 do 2 b=1,3 aldw( j,l,a,b)=ald( l,a,b) aldvw(j,l,a,b)=aldv(l,a,b) Gw( j,l,a,b)=G( l,a,b) do 2 c=1,3 2 Qw( j,l,a,b,c)=Q( l,a,b,c) endif return end subroutine initt(np,n3,alpha,alphav,GP,A, 1b2,alpha2,aG,bA,bG,alphaw,GPw,Aw,alphavw,ald,aldv,aldj,gpd,ad, 1aldw,aldvw,gpdw,adw) implicit none integer*4 np,n3 real*8 alpha(3,3),alphav(3,3),A(3,3,3),GP(3,3), 1b2(np),alpha2(np),aG(np),bA(np),bG(np), 1alphaw(np,3,3),GPw(np,3,3),Aw(np,3,3,3), 1alphavw(np,3,3),aldj(2*n3,3,3), 1ald(2*n3,3,3),aldv(2*n3,3,3),gpd(2*n3,3,3),ad(2*n3,3,3,3), 1aldw(np,2*n3,3,3),gpdw(np,2*n3,3,3),adw(np,2*n3,3,3,3), 1aldvw(np,2*n3,3,3) alpha=0.0d0 alphav=0.0d0 GP=0.0d0 A=0.0d0 alphaw=0.0d0 alphavw=0.0d0 GPw=0.0d0 Aw=0.0d0 b2=0.0d0 alpha2=0.0d0 AG=0.0d0 bA=0.0d0 bG=0.0d0 ald=0.0d0 aldv=0.0d0 aldj=0.0d0 gpd=0.0d0 ad=0.0d0 aldw=0.0d0 aldvw=0.0d0 gpdw=0.0d0 adw=0.0d0 return end subroutine dnst(nst,s) integer*4 nst,i character*80 filename character*(*) s nst=0 open(9,file=s) 1 read(9,90,end=99,err=99)filename 90 format(a80) nst=nst+1 goto 1 99 close(9) write(6,601)nst 601 format(i6,' lines in ',$) do 2 i=1,len(s) 2 write(6,602)s(i:i) 602 format(a1,$) write(6,*) c to avoid compiler messages: if(filename.eq.' ')call report('nic') return end subroutine rgsg(n3,gr,lgr,fgr) implicit none integer*4 n3,i,ix real*8 gr(n3),gn character*(*) fgr character*80 s logical lgr if(lgr)then open(8,file=fgr,status='old') 1 read(8,90,end=99,err=99)s 90 format(a80) if(s(38:59).eq.'Forces (Hartrees/Bohr)')then read(8,*) read(8,*) do 2 i=1,n3/3 2 read(8,*)(gr(ix+3*(i-1)),ix=1,2),(gr(ix+3*(i-1)),ix=1,3) gr=-gr gn=0.0d0 do 3 i=1,n3 3 gn=gn+gr(i)**2 gn=dsqrt(gn/dble(n3)) write(6,6001)gr(1),gr(n3),gn 6001 format(' Ground gradient G1 GN |G|:',3G12.4) close(8) goto 88 endif goto 1 99 close(8) call report('ground state gradient not found') else gr=0.0d0 endif 88 return end subroutine wcum(cum,alpha,GP,AA,ald,gpd,ad,n2,ie) implicit none integer*4 n2,a,b,c,i,ie real*8 cum(6),alpha(3,3),GP(3,3),AA(3,3,3),ald(n2,3,3), 1gpd(n2,3,3),ad(n2,3,3,3) cum=0.0d0 cum(1)=(alpha(1,1)+alpha(2,2)+alpha(3,3))/3.0d0 cum(2)=( GP(1,1)+ GP(2,2)+ GP(3,3))/3.0d0 do 1 a=1,3 do 1 b=1,3 do 1 c=1,3 1 cum(3)=cum(3)+2.0d0*AA(a,a,b)*(AA(b,c,c)+2.0d0*AA(c,b,c)) 1 +AA(a,b,b)*(AA(a,c,c)+2.0d0*AA(c,a,c)) 1 +2.0d0*AA(a,b,c)*(AA(a,b,c)+AA(b,a,c)+AA(c,a,b)) cum(3)=cum(3)/105.0d0 do 2 i=1,n2/2 cum(4)=cum(4)+ald(i,1,1 )**2+ald(i,2,2 )**2+ald(i,3,3 )**2 cum(5)=cum(5)+gpd(i,1,1 )**2+gpd(i,2,2 )**2+gpd(i,3,3 )**2 2 cum(6)=cum(6)+ad(i,1,1,1)**2+ad(i,2,2,2)**2+ad(i,3,3,3)**2 write(6,600)ie,cum,ald(1,1,1),gpd(1,1,1),ad(1,1,1,1) 600 format(' ##',i6,9e12.4) return end subroutine gen(nst,we,lddeg) implicit none integer*4 nst,i,j,r,rstate,nd real*8 we(nst),ev logical lddeg character*80 filename,s80 if(lddeg)then open(9,file='TD.LST') do 1 i=1,nst read(9,90)filename 90 format(a80) c determine the differentiated state: r=rstate(filename) open(8,file=filename) 11 read(8,80,end=88,err=88)s80 80 format(a80) if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s'. 1 and.r.gt.0)then do 20 j=1,79 if(s80(j:j) .eq.':') read(s80(15:j-1),*)nd 20 if(s80(j:j+1).eq.'eV')read(s80(j-10:j-2),*)ev if(nd.eq.r)we(i)=ev/27.211384205943d0 endif goto 11 88 close(8) 1 continue close(9) else we=0.0d0 endif return end subroutine gete(r,io,d,v,m,q,ie) implicit none integer*4 r,io,i1,i2,i,l,ie(8) real*8 d(3),v(3),m(3),q(6),x i1=-4 11 i1=i1+5 i2=min(i1+4,r) read(io,*) do 1 l=1,16 read(io,*)(x,i=i1-1,i2) if(i2.eq.r)then if(l.ge. 2.and.l.le. 4)d(l- 1)=x if(l.ge. 5.and.l.le. 7)v(l- 4)=x if(l.ge. 8.and.l.le.10)m(l- 7)=x if(l.ge.11 )q(l-10)=x endif 1 continue if(i2.lt.r)goto 11 do 2 i=2,5 2 ie(i)=ie(i)+1 return end SUBROUTINE READS(N3,S,E,M) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION S(N3,N3),E(N3) CM=219470.0d0 SFAC=0.0234280d0 S=0.0d0 E=0.0d0 OPEN(4,FILE='F.INP',STATUS='OLD') read(4,*)M,NAT,NAT do 1 i=1,NAT+1 1 read(4,*) DO 2 I=1,NAT DO 2 J=1,M 2 read(4,*)(s(3*(i-1)+ix,j),ix=1,2),(s(3*(i-1)+ix,j),ix=1,3) read(4,*) READ(4,4000)(E(I),I=1,M) 4000 FORMAT(6F11.6) close(4) c convert to atomic units: DO 5 J=1,M E(J)=E(J)/CM DO 5 I=1,N3 5 S(I,J)=S(I,J)*SFAC RETURN end subroutine addya(MQ,wi,w,g,E,qd,d,v,m,qm,n3,aldQ,aldvQ, 1aldjQ,gpQ,adQ,etd,radapt) implicit none integer*4 ii,n3,ix,jx,X,MQ,iq real*8 wi,w,g,d(3),v(3),m(3),qm(3,3),ux,uxd,etd(16,n3+3),E(n3), 1aldQ(MQ,2*n3,3,3),aldvQ(MQ,2*n3,3,3),aldjQ(MQ,2*n3,3,3), 1gpQ(MQ,2*n3,3,3),adQ(MQ,2*n3,3,3,3),uy,uyd,wo,vy,vyd,my,myd, 1fsq1,fsq2,fsq,fdq,wp,fvq,fiq1,fiq2,fmq,fiq,fwq,qd(n3,3,3), 1radapt do 1 iq=1,MQ c scattered frequency: wp=w-E(iq) c effective electronic frequency for this mode: wo=wi+radapt*E(iq) fsq1=(wo-w )/((wo-w )**2+g**2) fsq2=(wo+wp)/((wo+wp)**2+g**2) fsq=fsq1+fsq2 fdq=fsq1-fsq2 fvq=fdq*wo/w fiq1=-g/((wo-w )**2+g**2) fiq2=-g/((wo+wp)**2+g**2) fiq=fiq1+fiq2 fmq=fiq1-fiq2 fwq=fmq*wo/w do 1 ii=1,n3 do 1 ix=1,3 ux =d (ix) uxd=etd(ix+1,ii+3) do 1 jx=1,3 uy =d (jx) vy =v (jx) my =m (jx) uyd=etd(jx+1,ii+3) vyd=etd(jx+4,ii+3) myd=etd(jx+7,ii+3) c real parts of frequency functions: aldQ( iq,ii,ix,jx)=aldQ( iq,ii,ix,jx)+(uxd*uy+ux*uyd)*fsq aldvQ(iq,ii,ix,jx)=aldvQ(iq,ii,ix,jx)+(uxd*vy+ux*vyd)*fvq aldjQ(iq,ii,ix,jx)=aldjQ(iq,ii,ix,jx)+(uxd*vy+ux*vyd)*fsq c aldj(ii,ix,jx)=aldj(ii,ix,jx)+(uxd*vy+ux*vyd)*fv c gp: last index - magnetic gpQ(iq,ii,ix,jx)=gpQ(iq,ii,ix,jx)+(uxd*my+ux*myd)*fdq c ad - first index - dipole adQ(iq,ii,ix,jx,1)=adQ(iq,ii,ix,jx,1) 1+(uxd*qm(jx,1)+ux*qd(ii,jx,1))*fsq adQ(iq,ii,ix,jx,2)=adQ(iq,ii,ix,jx,2) 1+(uxd*qm(jx,2)+ux*qd(ii,jx,2))*fsq adQ(iq,ii,ix,jx,3)=adQ(iq,ii,ix,jx,3) 1+(uxd*qm(jx,3)+ux*qd(ii,jx,3))*fsq c complex: X=ii+n3 aldQ( iq,X,ix,jx)=aldQ( iq,X,ix,jx)+(uxd*uy+ux*uyd)*fiq aldvQ(iq,X,ix,jx)=aldvQ(iq,X,ix,jx)+(uxd*vy+ux*vyd)*fwq aldjQ(iq,X,ix,jx)=aldjQ(iq,X,ix,jx)+(uxd*vy+ux*vyd)*fiq gpQ( iq,X,ix,jx)=gpQ( iq,X,ix,jx)+(uxd*my+ux*myd)*fmq adq( iq,X,ix,jx,1)=adQ(iq,X,ix,jx,1) 1+(uxd*qm(jx,1)+ux*qd(ii,jx,1))*fiq adQ( iq,X,ix,jx,2)=adQ(iq,X,ix,jx,2) 1+(uxd*qm(jx,2)+ux*qd(ii,jx,2))*fiq 1 adQ( iq,X,ix,jx,3)=adQ(iq,X,ix,jx,3) 1+(uxd*qm(jx,3)+ux*qd(ii,jx,3))*fiq return end SUBROUTINE TRTQ(S,M,N,alQ,avQ,gdQ,aQ,E) IMPLICIT none integer*4 M,N,iq,ii,ix,iy,iz,i real*8 S(N,N),alQ(M,2*N,3,3),avQ(M,2*N,3,3),gdQ(M,2*N,3,3), 1aQ(M,2*N,3,3,3),E(N),si real*8,allocatable::anm(:,:,:,:),alnm(:,:,:),gnm(:,:,:), 1avnm(:,:,:),anmi(:,:,:,:),alnmi(:,:,:),gnmi(:,:,:), 1avnmi(:,:,:) allocate(anm( N,3,3,3),alnm( N,3,3),gnm( N,3,3),avnm( N,3,3)) allocate(anmi(N,3,3,3),alnmi(N,3,3),gnmi(N,3,3),avnmi(N,3,3)) anm=0.0d0 alnm=0.0d0 avnm=0.0d0 gnm=0.0d0 anmi=0.0d0 alnmi=0.0d0 avnmi=0.0d0 gnmi=0.0d0 c c real: do 1 iq=1,M do 1 ii=1,N si=S(ii,iq) do 1 ix=1,3 do 1 iy=1,3 alnm(iq,ix,iy)=alnm(iq,ix,iy)+si*alQ(iq,ii,ix,iy) avnm(iq,ix,iy)=avnm(iq,ix,iy)+si*avQ(iq,ii,ix,iy) gnm( iq,ix,iy)=gnm( iq,ix,iy)+si*gdQ(iq,ii,ix,iy) do 1 iz=1,3 1 Anm(iq,ix,iy,iz)=Anm(iq,ix,iy,iz)+si*aQ(iq,ii,ix,iy,iz) c c imaginary: do 2 iq=1,M do 2 ii=1,N i=ii+N si=S(ii,iq) do 2 ix=1,3 do 2 iy=1,3 alnmi(iq,ix,iy)=alnmi(iq,ix,iy)+si*alQ(iq,i,ix,iy) avnmi(iq,ix,iy)=avnmi(iq,ix,iy)+si*avQ(iq,i,ix,iy) gnmi( iq,ix,iy)=gnmi( iq,ix,iy)+si*gdQ(iq,i,ix,iy) do 2 iz=1,3 2 Anmi(iq,ix,iy,iz)=Anmi(iq,ix,iy,iz)+si*aQ(iq,i,ix,iy,iz) call WTQi(N,M,alnm,anm,gnm,avnm,alnmi,anmi,gnmi,avnmi,E) return end C ========================================== SUBROUTINE WTQi(N,M,ALPHA,A,G,av,ALPHAi,Ai,Gi,avi,E) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),av(N,3,3), 1ALPHAi(N,3,3),Gi(N,3,3),Ai(N,3,3,3),avi(N,3,3),E(N) c CM=219470.0d0 open(22,file='FILE.Q.TTT') WRITE(22,2003)M 2003 FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2 ' mode e(cm-1) jx jy jz') DO 1 II=1,M write(22,221)II,E(II)*CM 221 format(i5,f12.2) DO 1 I=1,3 1 WRITE(22,2001)I,(ALPHA(II,I,J),J=1,3),(ALPHAi(II,I,J),J=1,3) 2001 FORMAT(I3,6e14.6) WRITE(22,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' mode e(cm-1) jx(Bx) jy(By) jz(Bz)') DO 2 II=1,M write(22,221)II,E(II)*CM DO 2 I=1,3 2 WRITE(22,2001)I,(G(II,I,J),J=1,3),(Gi(II,I,J),J=1,3) c last inner index magnetic WRITE(22,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' mode e(cm-1) kx ky kz') DO 3 II=1,M write(22,221)II,E(II)*CM DO 3 I=1,3 DO 3 J=1,3 3 WRITE(22,2002)I,J,(A(II,I,J,K),K=1,3),(Ai(II,I,J,K),K=1,3) 2002 FORMAT(2I3,6e14.6) WRITE(22,2009) 2009 FORMAT('The velocity form of alpha',/, 1' mode e(cm-1) vx vy vz') DO 5 II=1,M write(22,221)II,E(II)*CM DO 5 I=1,3 5 WRITE(22,2001)I,(av(II,I,J),J=1,3),(avi(II,I,J),J=1,3) close(22) RETURN END subroutine escala(nst,n3,i,wi,d,v,m,q,etd,esc,dsc,rsc) implicit none integer*4 i,a,ix,nst,n3 real*8 d(3),m(3),q(6),etd(16,n3+3),fd,fr,esc(nst),dsc(nst), 1rsc(nst),lold,dold,rold,debye,CM,wi,v(3),clight debye=2.541580253D0 clight=137.03599d0 CM=219474.0d0 lold=1.0d7/(wi*CM) wi=(1.0d7/esc(i))/CM dold=(d(1)*d(1)+d(2)*d(2)+d(3)*d(3))*debye**2 rold=(d(1)*m(1)+d(2)*m(2)+d(3)*m(3))*debye**2/clight write(6,600)i,lold,esc(i),dold,dsc(i),-rold,rsc(i) 600 format(' Transition',i4,': e:',f12.1,' nm ->',f12.1,' nm',/, 1 ' ',4x,' d:',e12.3,' D ->',e12.3,' D ',/, 1 ' ',4x,' r:',e12.3,' D ->',e12.3,' D') if(dabs(dold).gt.1.0d-9)then fd=dsqrt(dsc(i)/dold) else call report('Zero dipole strength, cannot scale') endif if(dabs(rold).gt.1.0d-9)then if(dabs(dsc(i)).gt.1.0d-9)then fr=dsqrt(dold/dsc(i))*rsc(i)/rold else fr=0.0d0 endif else call report('Zero rotational strength, cannot scale') endif d=d*fd v=v*fd m=m*fr q=q*fd do 4 ix=1,3 do 4 a=1,n3 etd(1 +ix,3+a)=etd(1 +ix,3+a)*fd etd(4 +ix,3+a)=etd(4 +ix,3+a)*fd etd(7 +ix,3+a)=etd(7 +ix,3+a)*fr etd(10+ix,3+a)=etd(10+ix,3+a)*fd 4 etd(13+ix,3+a)=etd(13+ix,3+a)*fd return end c ==================================== SUBROUTINE DOMA(A,C,NAT,N6) c N6 - number of coordinates to be projected c 1..3 translations, 4..6 rotations IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION C(3,NAT),A(3*NAT,3*NAT) real*8,allocatable::TEM(:,:) allocate(TEM(3*NAT,3*NAT)) N=NAT*3 AMACH=1.0d-12 DO 12 I=1,N DO 12 J=1,N 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 A(3*I-2,4)=-C(2,I) A(3*I-1,4)= C(1,I) A(3*I-1,5)=-C(3,I) A(3*I ,5)= C(2,I) A(3*I-2,6)= C(3,I) 1 A(3*I ,6)=-C(1,I) CALL ORT(A,N,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1000.0d0*AMACH)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE DO 2 I=1,N DO 2 J=1,N TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END c =============================== SUBROUTINE ORT(A,NDIM,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(NDIM,NDIM) AMACH=1.0d-13 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/SQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP RETURN END