program vibpol implicit none integer*4 nm,nat,i,j,ntr,np real*8 wmin,wmax,g,w1min,w1max,debye,pi,CM, 1gammaau real*8 ,allocatable::u(:,:),m(:,:),uv(:,:),qv(:,:,:),w(:),eau(:) debye=2.541765d0 CM=219474.0d0 pi=4.0d0*datan(1.0d0) write(6,600) 600 format(' vibrational polarizabilities',/,/, 1 ' Input: F.INP',/, 1 ' VIBPOL.OPT',/, 1 ' Output: SOS.TTT') call readopt(w1min,w1max,wmin,wmax,np,g) open(9,file='F.INP') read(9,*)nm,nat,nat c transition dipoles allocate(u(nm,3),m(nm,3),w(nm),uv(nm,3),qv(nm,3,3),eau(nm)) call rtd(nm,nat,w,u,m,w1min,w1max) qv=0.0d0 c computation of polarizabilities from transition dipoles: do 1 i=1,nm eau(i)=w(i)/CM do 1 j=1,3 m(i,j)=m(i,j)/debye 1 u(i,j)=u(i,j)/debye uv=u gammaau=g/CM ntr=nm call cpolar('SOS.TTT','spav.prn','spal.prn','spgv.prn','spgl.prn', 1wmax,wmin,np,debye,pi,eau,gammaau,u,uv,m, 2qv,ntr,nm) end subroutine rtd(nm,nat,w,u,m,w1min,w1max) implicit none integer*4 nm,nat,i,j real*8 u(nm,3),m(nm,3),w1min,w1max,w(nm) do 1 i=1,nat+1+nm*nat+1 1 read(9,*) read(9,*)(w(i),i=1,nm) do 2 i=1,nm read(9,*)(u(i,j),j=1,3),(m(i,j),j=1,3) if(w(i).lt.w1min.or.w(i).gt.w1max)then do 3 j=1,3 u(i,j)=0.0d0 3 m(i,j)=0.0d0 endif 2 continue close(9) return end subroutine setv(w0,w,fr,fi,G) implicit none real*8 w,jm,fr,fi,g,w0,a,ff ff= (w0-w)*(w0+w) jm= ff**2+(G*w0)**2 c 2(w^2-w0^2) /[(w^2-w0^2)^2 + G^2 w0^2] c 2G /[(w^2-w0^2)^2 + G^2 w0^2] fr=2.0d0*ff / jm fi=2.0d0* G / jm return end subroutine readopt(w1min,w1max,wmin,wmax,np,g) implicit none integer*4 np real*8 wmin,wmax,w1min,w1max,g character*4 key c number of points/frequencies: np=300 c output frequency interval: wmin =1400.0d0 wmax =1700.0d0 c input frequency interval: w1min=1400.0d0 w1max=1700.0d0 c gamma/bandwidth: g=20.0d0 open(9,file='VIBPOL.OPT') 1 read(9,90,end=99,err=99)key 90 format(a4) if(key.eq.'WMIN')read(9,*)wmin if(key.eq.'1MIN')read(9,*)w1min if(key.eq.'WMAX')read(9,*)wmax if(key.eq.'1MAX')read(9,*)w1max if(key.eq.'NPOI')read(9,*)np if(key.eq.'GAMM')read(9,*)g goto 1 99 close(9) return end c ============================================================ subroutine cpolar(st,sav,sal,sgv,sgl,wmax,wmin,np,debye,pi,eau, 1gammaau,dl,dv,r,qv,ntr,n) implicit none integer*4 np,i,j,k,l,ie,ntr,n real*8 wmax,wmin,dw,cab,debye,pi,ccd,w,wau,al(3,3),ali(3,3), 1alit(3,3),allv(3,3),alilv(3,3),eau(*),gammaau,t,clight, 1dl(ntr,3),dv(ntr,3),r(ntr,3),qua,qv(ntr,3,3),fr,fi,CM character*(*) st,sav,sal,sgv,sgl CM=219474.0d0 write(6,6009) 6009 format(/, 1 ' Cpolar',/, 1 ' ^^^^^^',/) dw=(wmax-wmin)/dble(np-1) c constants making tensor traces epsilon c (factor of two missing for cab - do not know why): cab=debye**2 *108.7d0/pi c ccd=2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36/pi clight=137.036d0 ccd=debye**2*4.0d0*108.7d0/pi/clight open(44,file=st) c velocity-based trace: open(9,file=sav) c dipole-based trace: open(91,file=sal) c d and r spectra: open(81,file='DS.PRN') open(82,file='RS.PRN') write(44,446)np,wmin,wmax 446 format(i6,2f12.2) write(44,443) 443 format(' Dynamic polarizabilities') c alpha start ----------------------------------------------------- write(44,444) 444 format(10x,' Alpha ( f ) Alpha ( g) length-length ', 1'and length-velocity : ', 1/,'l(nm) Re ullx Re ully Re ullz Im ullx Im ully', 1' Im ullz Re ulvx Re ulvy Re ulvz Im ulvx Im ulvy', 1' Im ulvz') w=wmin-dw do 102 k=1,np w=w+dw wau=w/CM do 101 i=1,3 do 101 j=1,3 al(i,j)=0.0d0 ali(i,j)=0.0d0 alit(i,j)=0.0d0 allv(i,j)=0.0d0 alilv(i,j)=0.0d0 do 101 ie=1,n call setv(eau(ie),wau,fr,fi,gammaau) c a_ab = 2 sum wjn (1/(wjn**2-w**2)) Re ua ub c length-length: t=dl(ie,i)*dl(ie,j)*eau(ie) al(i,j) =al(i,j) +t*fr ali(i,j)=ali(i,j)+t*fi c length-velocity: t=dl(ie,i)*dv(ie,j)*eau(ie) allv(i,j) =allv(i,j) +t*fr alilv(i,j)=alilv(i,j)+t*fi c velocity-velocity (for spectra calculation only): 101 alit(i,j)=alit(i,j)+dv(ie,i)*dv(ie,j)*fi*eau(ie) write(9,69)w, 1wau**2*cab*(alit(1,1)+alit(2,2)+alit(3,3)), 1 al( 1,1)+al( 2,2)+al( 3,3) write(91,69)w, 1wau**2*cab*(ali(1,1)+ali(2,2)+ali(3,3)), 1 al(1,1)+al( 2,2)+al( 3,3) 69 format(f7.1,2g15.4) write(81,69)w,wau**2*cab*(ali(1,1)+ali(2,2)+ali(3,3)) write(44,441)w,(al( 1,j),j=1,3),(ali( 1,j),j=1,3), 1 (allv(1,j),j=1,3),(alilv(1,j),j=1,3) write(44,442) (al( 2,j),j=1,3),(ali( 2,j),j=1,3), 1 (allv(2,j),j=1,3),(alilv(2,j),j=1,3) 102 write(44,442) (al( 3,j),j=1,3),(ali( 3,j),j=1,3), 1 (allv(3,j),j=1,3),(alilv(3,j),j=1,3) 441 format(f8.2,12f12.4) 442 format(8x ,12f12.4) close(9) write(6,*)sav//' (velocity) written' close(91) write(6,*)sal//' (length ) written' c alpha end ----------------------------------------------------- c Gp start ----------------------------------------------------- write(44,445) 445 format(10x,' Gp/w ( f) Gp/w ( g )', 1/,'l(nm) ux uy uz') open(9 ,file=sgv) open(91,file=sgl) w=wmin-dw do 103 k=1,np w=w+dw wau=w/CM do 104 i=1,3 do 104 j=1,3 al( i,j)=0.0d0 ali( i,j)=0.0d0 do 104 ie=1,n call setv(eau(ie),wau,fr,fi,gammaau) c length: t=dl(ie,i)*r(ie,j) al( i,j) =al(i,j) +t*fr ali(i,j)=ali(i,j) +t*fi c velocity: 104 alit(i,j)=alit(i,j) +dv(ie,i)*r(ie,j)*fi write(91,69)w,( -ali(1,1)-ali( 2,2)-ali( 3,3))*wau**3*ccd, 1 al(1,1)+al( 2,2)+al( 3,3) write(82,69)w,( ali(1,1)+ali( 2,2)+ali( 3,3))*wau**3*ccd*clight write(9 ,69)w,(-alit(1,1)-alit(2,2)-alit(3,3))*wau**3*ccd, 1 al( 1,1)+al( 2,2)+al( 3,3) write(44,441)w,(al(1,j),j=1,3),(ali(1,j),j=1,3) write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 103 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) close(9) close(91) close(81) close(82) write(6,*)sgv//' written' write(6,*)sgl//' written' c Gp end ----------------------------------------------------- c A start ----------------------------------------------------- write(44,447) 447 format(10x,' A ( f) A ( g) ', 1/,'l(nm) ux uy uz') w=wmin-dw do 105 k=1,np w=w+dw wau=w/CM do 105 l=1,3 do 106 i=1,3 do 106 j=1,3 al( i,j)=0.0d0 ali(i,j)=0.0d0 do 106 ie=1,n if(i.eq.j)then qua=(3.0d0*qv(ie,i,j)-qv(ie,1,1)-qv(ie,2,2)-qv(ie,3,3))/2.0d0 else qua=1.5d0*qv(ie,i,j) endif call setv(eau(ie),wau,fr,fi,gammaau) al(i,j) =al( i,j) +dl(ie,l)*qua*fr*eau(ie) 106 ali(i,j)=ali(i,j) +dl(ie,l)*qua*fi*eau(ie) if(l.eq.1)then write(44,441)w,(al(1,j),j=1,3),(ali(1,j),j=1,3) else write(44,442) (al(1,j),j=1,3),(ali(1,j),j=1,3) endif write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 105 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) close(44) c A end ----------------------------------------------------- write(6,*)st//' written' return end c ============================================================