program grroa c c read Gaussian RROA output c implicit none integer*4 nq,nf,npx,ns0,i,nw,nwr,iargc,nws,nwe,nm,nn parameter (ns0=18) real*8 alr1(3,3),ali1(3,3),gr1(3,3),gi1(3,3), 1ar1(3,3,3),ai1(3,3,3),gcr1(3,3),gci1(3,3), 1acr1(3,3,3),aci1(3,3,3),avr1(3,3),avi1(3,3), 1THRC,wnm,wcm,cm,wau,wmin,wmax,fwhh,kelvin,ff logical lg,luseg,lusea,lstatic,ldo(ns0) real*8,allocatable::e(:),alr(:,:,:),ali(:,:,:),gr(:,:,:), 1gi(:,:,:),ar(:,:,:,:),ai(:,:,:,:),gcr(:,:,:),bf(:), 1gci(:,:,:),acr(:,:,:,:),aci(:,:,:,:),sr(:,:,:), 1ls(:,:,:),gs(:,:,:),as(:,:,:,:),ws(:) character*80 s npx=3951 THRC=1.0d-12 wnm=532.0d0 wcm=1.0d7/wnm cm=219470.0d0 wau=wcm/cm write(6,6001) 6001 format(' Reads Gaussian RROA output') call readopt('GRROA.OPT',luseg,lusea,lstatic,ff,ldo,fwhh, 1wmin,wmax,lg,kelvin) nm=1 allocate(ls(1,3,3),gs(1,3,3),as(1,3,3,3),ws(1)) if(iargc().ne.1)call report('Usage: grroa ') if(lstatic)then call rst(1,nm,ls,gs,as,ws,0) deallocate(ls,gs,as,ws) allocate(ls(nm,3,3),gs(nm,3,3),as(nm,3,3,3),ws(nm)) call rst(nm,nm,ls,gs,as,ws,1) endif call getarg(1,s) nq=nf(s) c maximal dimension in case some modes not assigned: nn=nq+nm allocate(e(nn),alr(nn,3,3),ali(nn,3,3),gr(nn,3,3),gi(nn,3,3), 1ar(nn,3,3,3),ai(nn,3,3,3),gcr(nn,3,3),gci(nn,3,3), 1acr(nn,3,3,3),aci(nn,3,3,3),bf(nn)) alr=0.0d0 ali=0.0d0 gr=0.0d0 gi=0.0d0 gcr=0.0d0 gci=0.0d0 ar=0.0d0 ai=0.0d0 acr=0.0d0 aci=0.0d0 bf=1.0d0 e=0.0d0 call gw(s,nw) if(nw.gt.1)then write(6,6002)nw 6002 format(' Which (1 -',i2,') frequency to read (negative=all)? ',$) read(5,*)nwr else nwr=1 endif if(nwr.lt.0)then nws=1 nwe=nw else nws=nwr nwe=nwr endif allocate(sr(npx,2,ns0)) do 3 nwr=nws,nwe call ge(nn,s,e,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci,wau,nwr,bf) if(lstatic)call adds(nn,nq,e,alr,gi,ar,gci,acr,nm,ls,gs,as,ws,ff, 1wau) wcm=wau*cm wnm=1.0d7/wcm call initab(wnm,ldo) sr=0.0d0 open(78,file='INV.TXT') do 1 i=1,nq call bz2(nn,i,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci, 1alr1,ali1,gr1,gi1,ar1,ai1,gcr1,gci1,acr1,aci1,avr1,avi1) 1 call wrram3(wau,0.0d0,e(i),THRC,i,wnm,wmin,wmax,npx,lg,fwhh, 1alr1,ali1,gr1,gi1,gcr1,gci1,ar1,ai1,acr1,aci1,avr1,avi1,sr, 1lusea,luseg,bf(i),ldo) close(78) call closetab(ldo) 3 call c33(wmin,wmax,npx,sr,kelvin,wnm,ldo) end c ======================================== subroutine report(s) character*(*) s write(6,*)s stop end c ======================================== subroutine readopt(s,luseg,lusea,lstatic,ff,ldo,fwhh, 1wmin,wmax,lg,kelvin) implicit none character*(*) s character*11 key integer*4 ns0,ie,i parameter (ns0=18) logical lex,lusea,luseg,lstatic,ldo(ns0),lg real*8 ff,fwhh,wmin,wmax,kelvin character*11 st(ns0) data st/'ICP_0 ','ICPx_90 ','ICPz_90 ','ICPs_90 ', 1 'ICPu_90 ','ICP_180 ','SCP_0 ','SCPx_90 ', 1 'SCPz_90 ','SCPs_90 ','SCPu_90 ','SCP_180 ', 1 'DCPI_0 ','DCPI_90 ','DCPI_180 ','DCPII_0 ', 1 'DCPII_90 ','DCPII_180 '/ wmin=50.0d0 wmax=4000.0d0 fwhh=10.0d0 lg=.false. kelvin=300.0d0 ldo=.true. inquire(file=s,exist=lex) luseg=.true. lusea=.true. lstatic=.false. ff=1.0d0 if(lex)then write(6,*)s open(9,file=s) 1 read(9,90,end=2,err=2)key write(6,90)key 90 format(a11) if(key(1:4).eq.'WMIN')read(9,*)wmin if(key(1:4).eq.'WMAX')read(9,*)wmax if(key(1:4).eq.'TEMP')read(9,*)kelvin if(key(1:4).eq.'GAUS')read(9,*)lg if(key(1:4).eq.'FWHH')read(9,*)fwhh if(key(1:4).eq.'USEA')read(9,*)lusea if(key(1:4).eq.'USEG')read(9,*)luseg if(key(1:4).eq.'STAT')read(9,*)lstatic if(key(1:4).eq.'FACT')read(9,*)ff do 4 ie=1,len(key) 4 if(key(ie:ie).eq.' ')goto 3 3 ie=ie-1 do 5 i=1,ns0 5 if(key(1:ie).eq.st(i)(1:ie))read(9,*)ldo(i) goto 1 2 close(9) endif return end subroutine ap3(is,ns0,ECM,a,s,wrin,wrax,npx,fwhh,lglg) implicit none integer*4 ns0,npx,i,is real*8 ECM,wrin,wrax,fwhh,w,dw,sf,a(2),s(npx,2,ns0),t logical lglg dw=(wrax-wrin)/(npx-1) w=wrin-dw do 1 i=1,npx w=w+dw t=sf(fwhh,lglg,w,ECM) s(i,1,is)=s(i,1,is)+t*a(1) 1 s(i,2,is)=s(i,2,is)+t*a(2) return end c ======================================== subroutine c33(wrin,wrax,npx,s,T,wnm,ldo) implicit none integer*4 npx,i,ns0,ii,ie,is parameter (ns0=18) real*8 wrin,wrax,s(npx,2,ns0),w,dw,T,wnm character*11 st(ns0),s11 logical ldo(ns0) data st/'ICP_0 ','ICPx_90 ','ICPz_90 ','ICPs_90 ', 1 'ICPu_90 ','ICP_180 ','SCP_0 ','SCPx_90 ', 1 'SCPz_90 ','SCPs_90 ','SCPu_90 ','SCP_180 ', 1 'DCPI_0 ','DCPI_90 ','DCPI_180 ','DCPII_0 ', 1 'DCPII_90 ','DCPII_180 '/ write(s11,111)wnm 111 format(f11.0) do 3 is=1,len(s11) 3 if(s11(is:is).ne.' ')goto 4 c4 tempcm=0.6950d0*T 4 do 100 ii=1,ns0 if(ldo(ii))then do 101 ie=len(st(ii)),1,-1 101 if(st(ii)(ie:ie).ne.' ')goto 102 102 open(45,file=s11(is:11)//st(ii)(1:ie)//'.ram.prn') dw=(wrax-wrin)/(npx-1) w=wrin-dw do 1 i=1,npx w=w+dw 1 write(45,4545)w,s(i,1,ii) c write(45,4545)w,s(i,1,ii)/(1.0d0-exp(-w/tempcm)) 4545 format(f12.3,' ',g12.4) close(45) open(45,file=s11(is:11)//st(ii)(1:ie)//'.roa.prn') dw=(wrax-wrin)/(npx-1) w=wrin-dw do 2 i=1,npx w=w+dw 2 write(45,4545)w,s(i,2,ii) c write(45,4545)w,s(i,2,ii)/(1.0d0-exp(-w/tempcm)) endif 100 continue return end c ============================================================== subroutine wrram3(w,ei,ef,THRC,nt,EXCNM,wrin,wrax,npx,lglg,fwhh, 1alr,ali,gtr,gti,gtcr,gtci,atr,ati,atcr,atci,avr,avi,sr, 1lusea,luseg,bf,ldo) c EXCNM ... excitation frequency in nm implicit none integer*4 nt,npx,ns0,ni0, 1ia,b,c,id,e,ii,is parameter (ns0=18,ni0=13) c ns0 : number of experimental setups c ni0 : number of invariants real*8 ef,ei,fwhh,THRC,AMU,BOHR,ECM,EXCNM,wrax,wrin,w,clight,CM, 1YDY,YDX,sr(npx,2,ns0),gpisvejc,EXCA,roa1,ram1,tr,ti,a(2),eps, 1inv(ni0),co(ni0+2,ns0),aaar,aaai,gtaar,gtaai,gtcaar,gtcaai, 1alr(3,3),ali(3,3),gtr(3,3),gti(3,3),gtcr(3,3),gtci(3,3), 1atr(3,3,3),ati(3,3,3),atcr(3,3,3),atci(3,3,3), 1avr(3,3),avi(3,3), 1alsr(3,3),alsi(3,3),gtsr(3,3),gtsi(3,3),gtcsr(3,3),gtcsi(3,3), 1avsr(3,3),avsi(3,3), 1alar(3,3),alai(3,3),gtar(3,3),gtai(3,3),gtcar(3,3),gtcai(3,3), 1avar(3,3),avai(3,3), 1ear(3,3),eapr(3,3),easr(3,3),eapsr(3,3),eaar(3,3),eapar(3,3), 1eai(3,3),eapi(3,3),easi(3,3),eapsi(3,3),eaai(3,3),eapai(3,3), 1eacr(3,3),eapcr(3,3),eacsr(3,3),eapcsr(3,3),eacar(3,3), 1eapcar(3,3),eaci(3,3),eapci(3,3),eacsi(3,3),eapcsi(3,3), 1eacai(3,3),eapcai(3,3),bf character*9 si(ni0) logical ldebug,lglg,lusea,luseg,ldo(ns0) data si/'a2 ','bs(a)2 ','ba(a)2 ','aG ', 1 'bs(G)2 ','ba(G)2 ','bs(A)2 ','ba(A)2 ', 1 'aGc ','bs(Gc)2 ','ba(Gc)2 ','bs(Ac)2 ', 1 'ba(Ac)2 '/ c prefix RAM, prefix ROA c a2 bs(a)2 ba(a)2 aG bs(G)2 ba(G)2 bs(A)2 c ba(A)2 aG bs(G)2 ba(G)2 bs(A)2 ba(A)2 c script: ^ ^ ^ ^ ^ c c co: Raman prefactor, ROA prefactor,3 Raman invariants, c 10 ROA invariants, all for ns0=18 spectral kinds c 1 1 2 ICP(0o): Nafie data co/ 4.0d0, 8.0d0, 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, 7.0d0, 5.0d0, 1.0d0, 1 -1.0d0,-45.0d0, 5.0d0, -5.0d0, -3.0d0, 1.0d0, c 2 3 4 ICPx(90o): Nafie 1 4.0d0, 2.0d0, 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, 7.0d0, 5.0d0, 1.0d0, 1 -1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, c 3 5 6 ICPz(90o): Nafie 1 8.0d0, 4.0d0, 1 0.0d0, 3.0d0, 5.0d0, 0.0d0, 3.0d0, 5.0d0, -1.0d0, 1 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, c 4 7 8 ICP*(90o): (magic) Nafie 1 13.3333d0,6.666d0, 1 9.0d0, 2.0d0, 2.0d0, 9.0d0, 2.0d0, 2.0d0, 0.0d0, 1 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, c 5 9 10 ICPu(90o): Nafie 1 4.0d0, 4.0d0, 1 45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0, 1 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, c 6 11 12 ICP(180o): Nafie 1 4.0d0, 8.0d0, 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, 7.0d0, 5.0d0, 1.0d0, 1 -1.0d0, 45.0d0, -5.0d0, 5.0d0, 3.0d0, -1.0d0, c 7 13 14 SCP(0o): nafie 1 4.0d0, 8.0d0, 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, -5.0d0, 5.0d0, -3.0d0, 1 -1.0d0,-45.0d0, -7.0d0, -5.0d0, 1.0d0, 1.0d0, c 8 15 16 SCPx(90o): nafie 1 2.0d0, 4.0d0, 1 45.0d0, 7.0d0, 5.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0,-45.0d0, -7.0d0, -5.0d0, 1.0d0, 1.0d0, c 9 17 18 SCPz(90o): Nafie 1 4.0d0, 8.0d0, 1 0.0d0, 3.0d0, 5.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0, 0.0d0, -3.0d0, -5.0d0, -1.0d0, -1.0d0, c 10 19 20 SCP*(90o): Nafie 1 13.3333d0,6.666d0, 1 9.0d0, 2.0d0, 2.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0, -9.0d0, -2.0d0, -2.0d0, 0.0d0, 0.0d0, c 11 21 22 SCPu(90o): Nafie 1 4.0d0, 4.0d0, 1 45.0d0, 13.0d0, 15.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0,-45.0d0,-13.0d0,-15.0d0, -1.0d0, -1.0d0, c 12 23 24 SCP(180o):Nafie 1 4.0d0, 8.0d0, 1 45.0d0, 7.0d0, 5.0d0,-45.0d0, 5.0d0, -5.0d0, 3.0d0, 1 1.0d0,-45.0d0, -7.0d0, -5.0d0, 1.0d0, 1.0d0, c 13 25 26 DCPI(0o): Nafie 1 4.0d0, 8.0d0, 1 45.0d0, 1.0d0, 5.0d0, 45.0d0, 1.0d0, 5.0d0, -1.0d0, 1 -1.0d0,-45.0d0, -1.0d0, -5.0d0, -1.0d0, +1.0d0, c 14 27 28 DCPI(90o): Nafie 1 2.0d0, 2.0d0, 1 45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0, 1 +1.0d0,-45.0d0,-13.0d0,-15.0d0, -1.0d0, -1.0d0, c 15 29 30 DCPI(180o): Nafie 1 24.0d0, 16.0d0, 1 0.0d0, 1.0d0, 0.0d0, 0.0d0, 3.0d0, 0.0d0, 1.0d0, 1 0.0d0, 0.0d0, -3.0d0, 0.0d0, 1.0d0, 0.0d0, c 16 31 32 DCPII(0o): Nafie 1 24.0d0, 16.0d0, 1 0.0d0, 1.0d0, 0.0d0, 0.0d0, 3.0d0, 0.0d0, 1.0d0, 1 0.0d0, 0.0d0, 3.0d0, 0.0d0, -1.0d0, 0.0d0, c 17 33 34 DCPII(90o): 1 2.0d0, 2.0d0, 1 45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0, 1 1.0d0, 45.0d0, 13.0d0, 15.0d0, 1.0d0, 1.0d0, c 18 35 36 DCPII(180o): 1 4.0d0, 8.0d0, 1 45.0d0, 1.0d0, 5.0d0, 45.0d0, 1.0d0, 5.0d0, -1.0d0, 1 -1.0d0, 45.0d0, 1.0d0, 5.0d0, 1.0d0, -1.0d0/ c prefix a2 bs(a)2 ba(a)2 aG bs(G)2 ba(G)2 c bs(A)2 ba(A)2 aG bs(G)2 ba(G)2 bs(A)2 ba(A)2 c script: ^ ^ ^ ^ ^ ^ c invariants: c 1 a2 =(1/9)Re (as_aa as*_bb) c 2 bs(a)2 =(1/2)Re(3as_ab as*_ab-as_aa as*_bb) c 3 ba(a)2 =(1/2)Re(3as_ab as*_ab-as_aa as*_bb) c 4 aG =(1/9)Im(as_aa Gs*_bb) c 5 bs(G)2 =(1/2)Im(3as_ab Gs*_ab-as_aa Gs*_bb) c 6 ba(G)2 =(3/2)Im(3aa_ab Ga*_ab) c 7 bs(A)2 =(w/2)Im(i as_ab eas*_ab) c ea_ab=eps_adg A_d,gb c 8 ba(A)2 =(w/2)Im(i aa_ab eaa*_ab)+i aa_ab eap_ab c eap_ab=eps_abg A_d,gd c script tensors: c 9 aG =(1/9)Im(as_aa Gs*_bb) c 10 bs(G)2 =(1/2)Im(3as_ab Gs*_ab-as_aa Gs*_bb) c 11 ba(G)2 =(3/2)Im(3aa_ab Ga*_ab) c 12 bs(A)2 =(w/2)Im(i as_ab eas*_ab) c 13 ba(A)2 =(w/2)Im(i aa_ab eaa*_ab)+i aa_ab eap_ab ldebug=.false. AMU=1822.0d0 BOHR=0.529177d0 EXCA=EXCNM*10.0d0 ECM=ef-ei CM=219474.0d0 clight=137.03599d0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA c eap_ab = eps_abc A_d,cd call calcep(eapr ,atr ) call calcep(eapi ,ati ) call calcep(eapci,atci) call calcep(eapcr,atcr) c c ea_ab = eps_adc A_d,cb ear =0.0d0 eai =0.0d0 eacr=0.0d0 eaci=0.0d0 do 2 ia=1,3 id=ia+1 if(id.gt.3)id=1 e=id+1 if(e.gt.3)e=1 do 2 b=1,3 ear( ia,b)=atr (id,e,b)-atr (e,id,b) eai( ia,b)=ati (id,e,b)-ati (e,id,b) eacr(ia,b)=atcr(id,e,b)-atcr(e,id,b) 2 eaci(ia,b)=atci(id,e,b)-atci(e,id,b) call dsa(alr,ali,alsr,alsi,alar,alai) call dsa(avr,avi,avsr,avsi,avar,avai) call dsa(avr,avi,avsr,avsi,avar,avai) call dsa(gtr,gti,gtsr,gtsi,gtar,gtai) call dsa(gtcr,gtci,gtcsr,gtcsi,gtcar,gtcai) call dsa(ear ,eai ,easr ,easi ,eaar ,eaai ) call dsa(eapr ,eapi ,eapsr ,eapsi ,eapar ,eapai ) call dsa(eacr ,eaci ,eacsr ,eacsi ,eacar ,eacai ) call dsa(eapcr,eapci,eapcsr,eapcsi,eapcar,eapcai) inv=0.0d0 c a2: (1/9) Re (als_aa als*_bb) aaar=alsr(1,1)+alsr(2,2)+alsr(3,3) aaai=alsi(1,1)+alsi(2,2)+alsi(3,3) inv(1)=(aaar*aaar+aaai*aaai)/9.0d0 c beta_s(alpha)2:(1/2)Re(3als_ab als*_ab-als_aa als*_bb) call abab(tr,ti,alsr,alsi,alsr,alsi) inv(2)=1.5d0*tr-4.5d0*inv(1) c beta_a(alpha)2:(3/2)Re(als_ab als*_ab) call abab(tr,ti,alar,alai,alar,alai) inv(3)=1.5d0*tr if(luseg)then c aG: (1/9) Im(als_aa gs*_bb) gtaar=gtsr(1,1)+gtsr(2,2)+gtsr(3,3) gtaai=gtsi(1,1)+gtsi(2,2)+gtsi(3,3) inv(4)=(aaai*gtaar-aaar*gtaai)/9.0d0 c beta_s(G)2:(1/2)Im(3als_ab Gs*_ab-als_aa Gs*_bb) call abab(tr,ti,alsr,alsi,gtsr,gtsi) inv(5)=1.5d0*ti-4.5d0*inv(4) c beta_a(G)2: (3/2)Im(ala_ab Ga*_ab) call abab(tr,ti,alar,alai,gtar,gtai) inv(6)=1.5d0*ti c aGc: gtcaar=gtcsr(1,1)+gtcsr(2,2)+gtcsr(3,3) gtcaai=gtcsi(1,1)+gtcsi(2,2)+gtcsi(3,3) inv(9)=(aaai*gtcaar-aaar*gtcaai)/9.0d0 c beta_s(Gc)2: call abab(tr,ti,alsr,alsi,gtcsr,gtcsi) inv(10)=1.5d0*ti-4.5d0*inv(9) c beta_a(Gc)2: (3/2)Im(ala_ab Gca*_ab) call abab(tr,ti,alar,alai,gtcar,gtcai) inv(11)=1.5d0*ti endif if(lusea)then c betasA2: (w/2)Im(i als_ab e_agd As*g,db) call abab(tr,ti,alsr,alsi,easr,easi) inv(7)=0.5d0*w*tr/clight c betaaA2: (w/2)Im(i[ala_ab[e_agd Aa*g,db+ e_abg Aa*d,gd) call abab(tr,ti,alar,alai,eaar,eaai) inv(8)=tr call abab(tr,ti,alar,alai,eapar,eapai) inv(8)=0.5d0*w*(inv(8)+tr)/clight c betasAc2: (w/2)Im(i als_ab e_agd Acs*g,db) call abab(tr,ti,alsr,alsi,eacsr,eacsi) inv(12)=0.5d0*w*tr/clight c betaaAc2: call abab(tr,ti,alar,alai,eacar,eacai) inv(13)=tr call abab(tr,ti,alar,alai,eapcar,eapcai) inv(13)=0.5d0*w*(inv(13)+tr)/clight endif if(ldebug.and.dabs(ECM-1691.62d0).lt.1)then c symmetric and antisymmetric tensor combinations write(6,*)ECM write(6,*)'alr:' write(6,545)alr*dsqrt(2*ECM/CM) 545 format(3f12.6) write(6,*)'gtr:' write(6,545)gtr*dsqrt(2*ECM/CM) write(6,*)'gti:' write(6,545)gti*dsqrt(2*ECM/CM) write(6,*)'gtci:' write(6,545)gtci*dsqrt(2*ECM/CM) write(6,*)'gtcr:' write(6,545)gtcr*dsqrt(2*ECM/CM) write(6,*)'Atr:' write(6,545)atr*dsqrt(2*ECM/CM)*w/clight write(6,*)'Ati:' write(6,545)ati*dsqrt(2*ECM/CM)*w/clight write(6,*)'Atcr:' write(6,545)atcr*dsqrt(2*ECM/CM)*w/clight write(6,*)'Atir:' write(6,545)atci*dsqrt(2*ECM/CM)*w/clight write(6,*)'alpha2',inv(1)*dsqrt(2*ECM/CM)*dsqrt(2*ECM/CM) write(6,*)'beta2',inv(2)*dsqrt(2*ECM/CM)*dsqrt(2*ECM/CM) write(6,*) write(6,*)'ag',inv(4)*dsqrt(2*ECM/CM)*dsqrt(2*ECM/CM) write(6,*)'agc',inv(9)*dsqrt(2*ECM/CM)*dsqrt(2*ECM/CM) write(6,*)'betaags',inv(5)*dsqrt(2*ECM/CM)*dsqrt(2*ECM/CM) write(6,*)'betaaga',inv(6)*dsqrt(2*ECM/CM)*dsqrt(2*ECM/CM) write(6,*)'betaagcs',inv(10)*dsqrt(2*ECM/CM)*dsqrt(2*ECM/CM) write(6,*) write(6,*)'betaA2s',inv(7)*2*ECM/CM write(6,*)'betaAa2',inv(8)*2*ECM/CM write(6,*)'betaA2cs',inv(12)*2*ECM/CM write(6,*)'betaA2ca',inv(13)*2*ECM/CM endif if(ldebug.and.nt.eq.37)then write(6,*)'Omega:',w/clight call wre('Eacsr:',eacsr) eacr=0.0d0 eaci=0.0d0 do 61 ia=1,3 do 61 b=1,3 do 61 c=1,3 do 61 id=1,3 eacr(ia,b)=eacr(ia,b)+eps(ia,c,id)*atcr(c,id,b) 61 eaci(ia,b)=eaci(ia,b)+eps(ia,c,id)*atci(c,id,b) call wre('Eacr:',eacr) call dsa(eacr ,eaci ,eacsr ,eacsi ,eacar ,eacai ) call wre('Eacsr:',eacsr) do 6 ia=1,3 do 6 is=1,3 6 write(6,*)ia,is,alr(ia,is),eacsr(ia,is),eacsi(ia,is) endif write(78,701)ECM 701 format(/,' Energy =',f12.4,' cm^-1') write(78,702)w*CM 702 format(' -> Omega =',f9.1,' cm^-1') do 7 ia=1,ni0 7 write(78,708)si(ia),inv(ia) 708 format(3x,a9,' =',e14.6) c spectral intensities, a1 Raman, a2 ROA: do 3 is=1,ns0 if(ldo(is))then a(1)=bf* 1 co(1,is)*(co(3,is)*inv(1)+co(4,is)*inv(2)+co(5,is)*inv(3)) a(2)=0.0d0 do 4 ii=4,ni0 4 a(2)=a(2)+bf*inv(ii)*co(2,is)*co(2+ii,is)/clight if(dabs(a(1)).gt.THRC)then call ap3(is,ns0,ECM,a,sr,wrin,wrax,npx,fwhh,lglg) YDX=0.0d0 YDY=0.0d0 ram1=a(1)*gpisvejc roa1=a(2)*gpisvejc WRITE(50+is,3000)nt,ECM,YDX,YDY,ram1,roa1 3000 FORMAT(I7,f9.2,2f3.0,g12.4,' 0 0 0 0 ',g12.4) endif endif 3 continue return end subroutine abab(tr,ti,ar,ai,br,bi) c t = a_ab b*_ab implicit none integer*4 a real*8 tr,ti,ar(3,3),ai(3,3),br(3,3),bi(3,3) tr=0.0d0 ti=0.0d0 do 1 a=1,3 tr=tr+ar(a,1)*br(a,1)+ai(a,1)*bi(a,1) 1 +ar(a,2)*br(a,2)+ai(a,2)*bi(a,2) 1 +ar(a,3)*br(a,3)+ai(a,3)*bi(a,3) 1 ti=ti-ar(a,1)*bi(a,1)+ai(a,1)*br(a,1) 1 -ar(a,2)*bi(a,2)+ai(a,2)*br(a,2) 1 -ar(a,3)*bi(a,3)+ai(a,3)*br(a,3) return end subroutine wre(s,e) character*(*) s real*8 e(3,3) write(6,'(a)')s write(6,609)e 609 format(3f12.6) return end subroutine dsa(tr,ti,tsr,tsi,tar,tai) implicit none real*8 tr(3,3),ti(3,3),tsr(3,3),tsi(3,3),tar(3,3),tai(3,3) integer*4 a,b do 1 a=1,3 do 1 b=1,3 tsr(a,b)=0.5d0*(tr(a,b)+tr(b,a)) tsi(a,b)=0.5d0*(ti(a,b)+ti(b,a)) tar(a,b)=0.5d0*(tr(a,b)-tr(b,a)) 1 tai(a,b)=0.5d0*(ti(a,b)-ti(b,a)) return end subroutine bz2(nn,n,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci, 1alr1,ali1,gr1,gi1,ar1,ai1,gcr1,gci1,acr1,aci1,avr1,avi1) implicit none integer*4 n,i,j,k,nn real*8 alr1(3,3),ali1(3,3),gr1(3,3),gi1(3,3), 1ar1(3,3,3),ai1(3,3,3),gcr1(3,3),gci1(3,3), 1acr1(3,3,3),aci1(3,3,3),avr1(3,3),avi1(3,3), 1alr(nn,3,3),ali(nn,3,3),gr(nn,3,3),gi(nn,3,3), 1ar(nn,3,3,3),ai(nn,3,3,3),gcr(nn,3,3),gci(nn,3,3), 1acr(nn,3,3,3),aci(nn,3,3,3) do 1 i=1,3 do 1 j=1,3 alr1(i,j)=alr(n,i,j) ali1(i,j)=ali(n,i,j) gr1( i,j)=gr( n,i,j) gi1( i,j)=gi( n,i,j) gcr1(i,j)=gcr(n,i,j) gci1(i,j)=gci(n,i,j) avr1(i,j)=alr(n,i,j) avi1(i,j)=ali(n,i,j) do 1 k=1,3 ar1( i,j,k)=ar( n,i,j,k) ai1( i,j,k)=ai( n,i,j,k) acr1(i,j,k)=acr(n,i,j,k) 1 aci1(i,j,k)=aci(n,i,j,k) return end subroutine initab(enm,ldo) implicit none integer*4 i,j,ns0,is parameter(ns0=18) real*8 enm logical ldo(ns0) character*11 st(ns0),s11 data st/'ICP_0 ','ICPx_90 ','ICPz_90 ','ICPs_90 ', 1 'ICPu_90 ','ICP_180 ','SCP_0 ','SCPx_90 ', 1 'SCPz_90 ','SCPs_90 ','SCPu_90 ','SCP_180 ', 1 'DCPI_0 ','DCPI_90 ','DCPI_180 ','DCPII_0 ', 1 'DCPII_90 ','DCPII_180 '/ write(s11,111)enm 111 format(f11.0) do 5 is=1,len(s11) 5 if(s11(is:is).ne.' ')goto 4 4 do 1 i=1,ns0 if(ldo(i))then do 2 j=len(st),1,-1 2 if(st(i)(j:j).ne.' ')goto 3 3 open(50+i,file=s11(is:11)//st(i)(1:j)//'.tab') write(50+i,3000)st(i) 3000 format(20x,a11,/,' n wavelength (cm) Raman ROA',/,80(1h-)) endif 1 continue return end subroutine closetab(ldo) implicit none integer*4 i,ns0 logical ldo(*) parameter(ns0=18) do 1 i=1,ns0 if(ldo(i))then write(50+i,3000) 3000 format(80(1h-)) close(50+i) endif 1 continue return end subroutine gw(s,nw) implicit none integer*4 nw,n real*8 e,eold,wcm character*80 s80 character*(*) s c howm many excitation frequencies are present? n=0 nw=0 eold=0.0d0 open(9,file=s) 1 read(9,80,end=99,err=99)s80 80 format(a80) if(s80(1:10).eq.' Energy = ')then read(s80(10:21),*)e read(9,80,end=99,err=99)s80 read(s80(14:22),*)wcm if(n.eq.0)then eold=e n=n+1 endif if(dabs(e-eold).lt.1.0d-2)then if(n.eq.1)then nw=nw+1 write(6,601)nw,wcm,1.0d7/wcm 601 format(i6,f10.1,' cm-1',f10.1,' nm') endif else n=n+1 endif eold=e endif goto 1 99 close(9) write(6,600)n,nw 600 format(i6,' modes,',i6,' excitation frequencies') return end subroutine ge(nn,s,e,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci,wau, 1nwr,bf) implicit none integer*4 n,nwr,nw,nn real*8 e(nn),alr(nn,3,3),ali(nn,3,3),gr(nn,3,3),gi(nn,3,3), 1ar(nn,3,3,3),ai(nn,3,3,3),gcr(nn,3,3),gci(nn,3,3),eold,er, 1acr(nn,3,3,3),aci(nn,3,3,3),wau,CM,wcm,bf(*) character*(*) s character*80 s80 logical lf(5),lr CM=219470.0d0 n=0 nw=0 lf=.false. eold=0.0d0 open(9,file=s) 1 read(9,80,end=99,err=99)s80 80 format(a80) if(s80(1:10).eq.' Energy = ')then read(s80(10:21),*)er read(9,80,end=99,err=99)s80 read(s80(14:22),*)wcm if(n.eq.0)then eold=er n=n+1 endif if(dabs(er-eold).lt.1.0d-2)then nw=nw+1 else n=n+1 nw=1 endif lr=nw.eq.nwr if(lr)wau=wcm/CM e(n)=er if(s80(30:39).eq.'Boltzmann:')then read(s80(40:len(s80)),*)bf(n) else bf(n)=1.0d0 endif eold=er endif if(lr)then if(s80(6:36).eq.'ELECTRIC DIPOLE-ELECTRIC DIPOLE')then call r33(nn,n,alr,ali) lf(1)=.true. endif if(s80(6:36).eq.'ELECTRIC DIPOLE-MAGNETIC DIPOLE')then call r33(nn,n,gr,gi) lf(2)=.true. endif if(s80(6:36).eq.'MAGNETIC DIPOLE-ELECTRIC DIPOLE')then call r33(nn,n,gcr,gci) lf(3)=.true. endif if(s80(4:38).eq.'ELECTRIC DIPOLE-ELECTRIC QUADRUPOLE')then call r39(nn,n,ar,ai) lf(4)=.true. endif if(s80(4:38).eq.'ELECTRIC QUADRUPOLE-ELECTRIC DIPOLE')then call r93(nn,n,acr,aci) lf(5)=.true. endif endif goto 1 99 close(9) c a,ac:first index electric dipole if( lf(1))write(6,*)' alpha' if(.not.lf(1))write(6,*)' alpha not found' if( lf(2))write(6,*)' G ' if(.not.lf(2))write(6,*)' G not found' if( lf(3))write(6,*)' Gc ' if(.not.lf(3))write(6,*)' Gc not found' if( lf(4))write(6,*)' A ' if(.not.lf(4))write(6,*)' A not found' if( lf(5))write(6,*)' Ac ' if(.not.lf(5))write(6,*)' Ac not found' if(lf(1).and.lf(2).and.lf(3).and.lf(4).and.lf(5))return end subroutine r93(nq,n,ar,ai) c quadrupole-dipole implicit none integer*4 n,nq,i,j,k real*8 ar(nq,3,3,3),ai(nq,3,3,3) character*80 s80 read(9,*) read(9,*) c electric dipole, keep it in the first index: do 1 i=1,3 do 1 j=1,3 do 1 k=1,j read(9,80)s80 80 format(a80) read(s80(7:40),*)ar(n,i,j,k),ai(n,i,j,k) ar(n,i,k,j)=ar(n,i,j,k) 1 ai(n,i,k,j)=ai(n,i,j,k) return end subroutine r39(nq,n,ar,ai) c dipole-quadrupoe implicit none integer*4 n,nq,i,j,k real*8 ar(nq,3,3,3),ai(nq,3,3,3) character*80 s80 read(9,*) read(9,*) do 1 i=1,3 do 1 k=1,3 do 1 j=1,k read(9,80)s80 read(s80(7:40),*)ar(n,i,j,k),ai(n,i,j,k) ar(n,i,k,j)=ar(n,i,j,k) 1 ai(n,i,k,j)=ai(n,i,j,k) 80 format(a80) return end subroutine r33(nq,n,alr,ali) implicit none integer*4 j,i,n,nq real*8 alr(nq,3,3),ali(nq,3,3) character*80 s80 read(9,*) read(9,*) do 1 j=1,3 do 1 i=1,3 read(9,80)s80 1 read(s80(7:40),*)alr(n,i,j),ali(n,i,j) 80 format(a80) return end function nf(s) implicit none integer*4 nf,n,nl character*(*) s character*80 s80 nl=0 n=0 open(9,file=s) 1 read(9,80,end=99,err=99)s80 80 format(a80) if(s80(1:10).eq.' Energy = ')n=n+1 nl=nl+1 goto 1 99 close(9) write(6,600)nl,n 600 format(i8,' lines,',i6,' modes') nf=n return end c ================================================================== 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 c ============================================================== FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END subroutine calcep(e,a) c eap_ab = eps_abc A_d,cd implicit none real*8 e(3,3),a(3,3,3) e(1,1)=0.0d0 e(2,2)=0.0d0 e(3,3)=0.0d0 e(1,2)= a(1,3,1)+a(2,3,2)+a(3,3,3) e(1,3)=-a(1,2,1)-a(2,2,2)-a(3,2,3) e(2,3)= a(1,1,1)+a(2,1,2)+a(3,1,3) e(2,1)=-e(1,2) e(3,1)=-e(1,3) e(3,2)=-e(2,3) return end subroutine rst(nm0,nm,ls,gs,as,ws,ic) implicit none integer*4 nm0,nm,ic,i,ix,iy,iz real*8 ls(nm0,3,3),gs(nm0,3,3),as(nm0,3,3,3),ws(nm0) open(9,file='FILE.Q.TTT',status='old') read(9,*) read(9,*)nm if(ic.eq.0)then close(9) return endif read(9,*) read(9,*) do 9 i=1,nm read(9,*)ws(i),ws(i) do 9 ix=1,3 9 read(9,*)ls(i,ix,1),(ls(i,ix,iy),iy=1,3) read(9,*) read(9,*) do 91 i=1,nm read(9,*) do 91 ix=1,3 c magnetic index iy, last 91 read(9,*)gs(i,ix,1),(gs(i,ix,iy),iy=1,3) read(9,*) read(9,*) do 92 i=1,nm read(9,*) do 92 ix=1,3 do 92 iy=1,3 c el index ix 92 read(9,*)(as(i,ix,iy,iz),iz=1,2),(as(i,ix,iy,iz),iz=1,3) close(9) write(6,*)'static polarizability read' return end subroutine adds(nn,nq,e,lr,gi,ar,gci,acr,nm,ls,gs,as,ws,ff, 1w) integer*4 nq,nm,i,j,ix,iy,iz,nn,ji,na real*8 ls(nm,3,3),gs(nm,3,3),as(nm,3,3,3),ws(nm), 1e(nn),lr(nn,3,3),gi(nn,3,3),ar(nn,3,3,3),gci(nn,3,3), 1acr(nn,3,3,3),t,ff,CM,f,d,w,clight,g integer*4 ,allocatable::asi(:) clight=137.03599d0 CM=219474.0d0 allocate(asi(nn)) asi=0 c na ... number of unassigned freuqencies na=0 do 1 i=1,nm f=ff/dsqrt(2.0d0*ws(i)/CM) g=f*clight/w ji=0 d=1.0d9 do 2 j=1,nq t=dabs(ws(i)-e(j)) if(t.lt.1.0d0.and.t.lt.d.and.asi(j).eq.0)then ji=j d=t endif 2 continue if(ji.eq.0)then write(6,600)i,ws(i) 600 format(' Frequency',i6,f12.2,' cm-1 not assigned') na=na+1 ji=nq+na e(ji)=ws(i) do 4 ix=1,3 do 4 iy=1,3 lr( ji,ix,iy )=ls(i,ix,iy )*f c gs is supposedly static limit, G'/w, so here it must be c multiplied by w: gi( ji,ix,iy )=-gs(i,ix,iy )*f*w gci(ji,ix,iy )= gs(i,ix,iy )*f*w do 4 iz=1,3 ar( ji,ix,iy,iz)=as(i,ix,iy,iz)*g 4 acr(ji,ix,iy,iz)=as(i,ix,iy,iz)*g else asi(ji)=1 do 3 ix=1,3 do 3 iy=1,3 lr( ji,ix,iy )=lr( ji,ix,iy )+ls(i,ix,iy )*f gi( ji,ix,iy )=gi( ji,ix,iy )-gs(i,ix,iy )*f*w gci(ji,ix,iy )=gci(ji,ix,iy )+gs(i,ix,iy )*f*w do 3 iz=1,3 ar( ji,ix,iy,iz)=ar( ji,ix,iy,iz)+as(i,ix,iy,iz)*g 3 acr(ji,ix,iy,iz)=acr(ji,ix,iy,iz)+as(i,ix,iy,iz)*g endif 1 continue nq=nq+na return end