program polderdip implicit none integer*4 i,iu,n,nat,n3,iud,np,iw,n2,nst 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),GPv(3,3), 1nmstart,nmend,dnm,evcm,nm real*8,allocatable::etd(:,:),ald(:,:,:),gpd(:,:,:),ad(:,:,:,:), 1qd(:,:,:),x(:),alpha2(:),b2(:),aG(:),bA(:),bG(:),ge(:), 1Aw(:,:,:,:),GPw(:,:,:),alphaw(:,:,:),alphavw(:,:,:), 1aldw(:,:,:,:),gpdw(:,:,:,:),adw(:,:,:,:,:),GPvw(:,:,:), 1aldvw(:,:,:,:),aldv(:,:,:), 1als(:,:,:),gs(:,:,:),as(:,:,:,:) integer*4,allocatable::lst(:) character*80 filename logical lint,lstat,lnm,lsep,ldog,usedd call readopt(width,excnm,evmin,evmax,np,lint,lstat,lnm, 1lsep,ldog,CM,evcm,w,g,dev,nmstart,nmend,dnm,nat,usedd) 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),GPvw(np,3,3),etd(16,n3+3), 1qd(n3,3,3),x(n3),ge(n3), 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,GPv,A, 1b2,alpha2,aG,bA,bG,alphaw,GPw,Aw,alphavw,GPvw,ald,aldv,gpd,ad, 1aldw,aldvw,gpdw,adw) call dnst(nst) allocate(lst(nst)) open(9,file='TD.LST') do 1 i=1,nst iu=0 iud=0 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,usedd) write(6,605)filename(1:50),i,n 605 format(1x,a50,2i3) c redefine quadrupoles: call qqm(q,qm) call qqmd(etd,qd,n3) c polarizabilities: iu=iu+1 f1=2.0d0*fi(wi,w,g) f2=-2.0d0*wi*fi(wi,w,g) c add this state contribution to polarizabilities, for w call add(f1,f2,d,v,m,qm,alpha,alphav,GP,GPv,A) 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,GPvw,Aw) endif c derivatives: iud=iud+1 c add this state contribution to polarizability derivatives call addd(d,v,m,qm,n3,ald,aldv,gpd,ad,etd,qd,lstat,wi,w,g,ge) 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) if(lsep) call drop(n3,i,wi,w,g,d,m,etd,qd,qm,ge) 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 wral(alphav,'POLV.TTT',w) call wra(a,w) call wrg(GP,'GP.TTT',w) call wrg(GPv,'GPv.TTT',w) if(ldog)call roadog(n3,x,ald,aldv,gpd,ad) if(lint.and.ldog)call roadogw(np,n3,x,aldw,aldvw,gpdw,adw) 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 call WRITETTT('FILE.TTT',N3,ald,aldv,gpd,ad,w) 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 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')read(s80(15:18),*)actual if(s80(2:28).eq.'This state for optimization')then r=actual goto 88 endif goto 1 88 close(8) if(r.eq.0)call report('Cannot determine excited state') rstate=r return end subroutine gd(io,filename,wi,d,r,v,m,q,n3,etd,x,ge,lst,usedd) 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(*) real*8 wi,d(3),ev,v(3),m(3),q(6),etd(16,n3+3),x(*),bohr,ge(*),gn logical usedd character*14 s8(8) data s8/ 'geometry ', 'dipoles ', 'v-dipoles ', 1 'm-dipoles ', 'quadrupoles ', 'Excited state ', 1 'derivatives ', 'gradient '/ 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 c determine the differentiated state in file io: r=rstate(filename) lst(io)=r do 101 I=1,io-1 101 if(lst(I).eq.r)call report('state defined multiple times') 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:'.OR. 2 s80(20:40).EQ.'Standard orientation:'.OR. 2 s80(26:46).EQ.'Standard 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 read the last dipole moment. etc: if(s80(2:65).eq.'Ground to excited state transition electric dipol 1e moments (Au):'.and.r.gt.0)then do 5 i=1,r 5 read(8,*) read(8,*)d(1),(d(ix),ix=1,3) ie(2)=ie(2)+1 endif if(s80(2:65).eq.'Ground to excited state transition velocity dipol 1e moments (Au):'.and.r.gt.0)then do 11 i=1,r 11 read(8,*) read(8,*)v(1),(v(ix),ix=1,3) ie(3)=ie(3)+1 endif if(s80(2:65).eq.'Ground to excited state transition magnetic dipol 1e moments (Au):'.and.r.gt.0)then do 14 i=1,r 14 read(8,*) read(8,*)m(1),(m(ix),ix=1,3) ie(4)=ie(4)+1 endif if(s80(2:69).eq.'Ground to excited state transition velocity quadr 1upole moments (Au):'.and.r.gt.0)then do 17 i=1,r 17 read(8,*) read(8,*)q(1),(q(ix),ix=1,6) c XX YY ZZ XY XZ YZ ie(5)=ie(5)+1 endif 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 if(s80(2:34).eq.'Electronic Transition Derivatives'. 1 and.ie(7).eq.0)then write(6,*)s80(2:34) 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 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 82 i=1,8 82 write(6,607)s8(i),ie(i) 607 format(a14,i4,' x') do 81 i=1,8 if(ie(i).eq.0)then if(i.ne.7.or.usedd)call report(' Incomplete output') endif 81 continue do 4 ix=1,3 m(ix)=m(ix)/2.0d0 4 v(ix)=-v(ix)/wi do 6 ix=1,6 6 q(ix)=-q(ix)/wi if(usedd)then do 8 i=1,n3 do 7 ix=1,3 etd(4 +ix,3+i)=-etd(4 +ix,3+i)/wi 7 etd(7 +ix,3+i)= etd(7 +ix,3+i)/2.0d0 do 8 ix=1,6 8 etd(10+ix,3+i)=-etd(10+ix,3+i)/wi else etd=0.0d0 if(ie(7).eq.0)then write(6,*)'Dipole derivatives not used' else write(6,*)'Dipole derivatives set to zero' endif endif 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,GD,ad,fr) 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) character*(*) fn 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)', 1 ' gradient alternative:') 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(ii,I,J),J=1,3),(GD(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(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(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 add(f1,f2,d,v,m,qm,alpha,alphav,GP,GPv,A) implicit none integer*4 i,j real*8 f1,f2,d(3),v(3),m(3),qm(3,3),GP(3,3),alpha(3,3),A(3,3,3), 1GPv(3,3),alphav(3,3) 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)+f2*d(i)*v(j) GP( i,j)=GP( i,j)+f1*d(i)*m(j) GPv( i,j)=GPv( i,j)+f1*v(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,GPvw,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),GPvw(np,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) GPvw( iw,ix,jx) =GPvw(iw,ix,jx)+f1*v(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) 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 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)+f3*ux*uy*gd aldv(ii,ix,jx)=aldv(ii,ix,jx)+f2*(uxd*vy+ux*vyd)+f3*ux*vy*gd ald( n3+ii,ix,jx)=ald (n3+ii,ix,jx)+g2*(uxd*uy+ux*uyd)+g3*ux*uy*gd aldv(n3+ii,ix,jx)=aldv(n3+ii,ix,jx)+g2*(uxd*vy+ux*vyd)+g3*ux*vy*gd c gpd: last index - magnetic gpd( ii,ix,jx)=gpd( ii,ix,jx)+f1*(uxd*my+ux*myd)+f4*ux*my*gd gpd(n3+ii,ix,jx)=gpd(n3+ii,ix,jx)+g1*(uxd*my+ux*myd)+g4*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))+f3*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))+g3*ux*qm(jx,kx)*gd 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 subroutine addd(d,v,m,qm,n3,ald,aldv,gpd,ad,etd,qd,lstat,wi,w,g, 1ge) 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, 1g1,g2,g3,g4 logical lstat call setf(1,f1,f2,f3,f4,g1,g2,g3,g4,wi,w,g) call addx(f1,f2,f3,f4,g1,g2,g3,g4, 1d,v,m,qm,n3,ald,aldv,gpd,ad,etd,qd,ge) 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) 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) implicit none integer*4 np,n3,iw,ic logical lnm 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) 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) 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 call axx(d,v,m,qm,np,n3,ald,alv,gpd,ad,etd,qd,g, 1evmin,dev,wi,1,lnm,ge,dnm) if(lstat)call axx(d,v,m,qm,np,n3,ald,alv,gpd,ad,etd,qd,g, 1evmin,dev,wi,2,lnm,ge,dnm) return end subroutine readopt(width,excnm,evmin,evmax,np,lint,lstat, 1lnm,lsep,ldog,CM,evcm,w,g,dev,nmstart,nmend,dnm,nat,usedd) implicit none real*8 width,excnm,evmin,evmax,CM,evcm,w,g,dev,nmstart,nmend, 1dnm integer*4 np,nat logical lint,lex,lstat,lnm,lsep,ldog,usedd character*80 s80,filename CM=219474.0d0 evcm=8065.73d0 c electronic band width in cm-1 width=700.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. 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.'WIDTH')read(9,*)width if(s80(1:5).eq.'EXCNM')read(9,*)excnm 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.'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 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 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) 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 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,g,a,w) 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 open(9,file=f) write(9,900)' Polarizability',w 900 format(a15,' w = ',F12.6,/, 1' 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 roadog(n3,x,ald,aldv,G,Q) implicit none integer*4 n3,a,b,c,L,IE,i 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(2,3,3),R(3) DO 1 L=0,n3/3-1 DO 1 IE=1,3 i=3*L+IE do 2 a=1,3 R(a)=X(a+3*L) do 2 b=1,3 d(1,a,b)=aldv( i,a,b)-ald( i,a,b) 2 d(2,a,b)=aldv(n3+i,a,b)-ald(n3+i,a,b) DO 1 a=1,3 G(i,a,1)=G(i,a,1)+0.5d0*(R(2)*d(1,a,3)-R(3)*d(1,a,2)) G(i,a,2)=G(i,a,2)+0.5d0*(R(3)*d(1,a,1)-R(1)*d(1,a,3)) G(i,a,3)=G(i,a,3)+0.5d0*(R(1)*d(1,a,2)-R(2)*d(1,a,1)) G(n3+i,a,1)=G(n3+i,a,1)+0.5d0*(R(2)*d(2,a,3)-R(3)*d(2,a,2)) G(n3+i,a,2)=G(n3+i,a,2)+0.5d0*(R(3)*d(2,a,1)-R(1)*d(2,a,3)) G(n3+i,a,3)=G(n3+i,a,3)+0.5d0*(R(1)*d(2,a,2)-R(2)*d(2,a,1)) DO 1 b=1,3 Q(i,a,b,b)=Q(i,a,b,b)+R(1)*d(1,a,1)+R(2)*d(1,a,2)+R(3)*d(1,a,3) Q(n3+i,a,b,b)= 1Q(n3+i,a,b,b) +R(1)*d(2,a,1)+R(2)*d(2,a,2)+R(3)*d(2,a,3) DO 1 c=1,3 Q( i,a,b,c)=Q( i,a,b,c)-1.5d0*(R(b)*d(1,a,c)+R(c)*d(1,a,b)) 1 Q(n3+i,a,b,c)=Q(n3+i,a,b,c)-1.5d0*(R(b)*d(2,a,c)+R(c)*d(2,a,b)) return end subroutine roadogw(np,n3,x,aldw,aldvw,Gw,Qw) 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) 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,G,Q) 1 call tofrom(j,1,np,n3,aldw,aldvw,Gw,Qw,ald,aldv,G,Q) return end subroutine wra(a,w) implicit none real*8 a(3,3,3),w integer*4 ix,jx,kx open(9,file='A.TTT') write(9,900)' A, dipole-quad',w 900 format(a15,' w = ',F12.6,/, 1' 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 open(9,file=s) write(9,901)' G-tensor ',w 901 format(a15,' w = ',F12.6,/, 1' 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 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,GPv,A, 1b2,alpha2,aG,bA,bG,alphaw,GPw,Aw,alphavw,GPvw,ald,aldv,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), 1GPv(3,3),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),GPvw(np,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 GPv=0.0d0 A=0.0d0 alphaw=0.0d0 alphavw=0.0d0 GPw=0.0d0 GPvw=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 gpd=0.0d0 ad=0.0d0 aldw=0.0d0 aldvw=0.0d0 gpdw=0.0d0 adw=0.0d0 return end subroutine dnst(nst) integer*4 nst character*80 filename nst=0 open(9,file='TD.LST') 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 TD.LST') if(filename.eq.' ')call report('nic') return end