PROGRAM PESSP c spectra generation IMPLICIT none real*8 qmin,qmax,dq,ee,pp,t,tol,EXCNM integer*4 np,i,n,ipl,ix,j,k real*8,allocatable::e(:),p(:),uo(:,:),psi(:,:) real*8,allocatable::al(:,:,:) integer*4,allocatable::ip(:) logical lex write(6,6000) 6000 format(' Generation of spectra from peswf output ',/,/, 1' Input : ENNTH.LST E(Q) ',/, 2' DIPOLE.TXT u(Q) ',/, 3' ALPHA.TXT a(Q) ',/, 3' N.PRN wavefunctions, N=1,2,... ',/, 3' SURF.OPT options ',/, 3' ENERGIES.TXT energies ',/,/, 3' Output: DOG.TAB IR ',/, 3' ROA.TAB Raman ',/,/) np=0 open(10,file='ENNTH.LST',status='old') 1 read(10,*,end=901,err=901)qmax,ee np=np+1 if(np.eq.1)qmin=qmax goto 1 901 close(10) if(np.lt.2)call report('not enough points in ENNTH.LST') dq=(qmax-qmin)/dble(np-1) write(6,6001)np,qmin,qmax,dq 6001 format(i10,' points in ENNTH.LST',/, 1 10x,' Qmin: ',f12.5,' Qmax: ',f12.5,' dQ: ',f12.5) n=0 do 2 i=1,2 j=0 open(10,file='ENERGIES.TXT',status='old') 3 read(10,*,end=902,err=902)ee,ee,pp j=j+1 if(i.eq.2)then e(n-j+1)=ee p(n-j+1)=pp endif goto 3 902 close(10) if(i.eq.1)then n=j write(6,*)n,' energies' allocate(e(n),p(n)) endif 2 continue write(6,*)' Energy probability:' do 4 i=1,n 4 write(6,6002)i,e(i),p(i) 6002 format(i4,f12.2,f12.6) allocate(ip(n)) call getpar(t,ipl,ip,tol,EXCNM) write(6,6003)t 6003 format(' temperatue (K):',f12.2) write(6,6004)ipl 6004 format(' number of states involved:',i5,':') write(6,6005)(ip(i),i=1,ipl) 6005 format(6i10) allocate(psi(np,ipl)) call getpsi(psi,np,ipl,ip) call testpsi(np,psi,ipl,dq) allocate(uo(np,3)) open(10,file='DIPOLE.TXT',status='old') do 5 i=1,np 5 read(10,*)uo(i,1),(uo(i,ix),ix=1,3) close(10) call dodog(e,p,np,psi,uo,ipl,tol,dq) inquire(file='AL.TXT',exist=lex) if(lex)then allocate(al(np,3,3)) open(10,file='AL.TXT',status='old') do 6 i=1,np read(10,*)((al(i,j,k),j=k,3),k=1,3) do 6 k=1,3 do 6 j=k,3 6 al(i,k,j)=al(i,j,k) close(10) call doroa(t,e,p,np,psi,al,ipl,tol,dq,EXCNM) endif end c ------------------------------------------------------------ subroutine testpsi(np,psi,ipl,dq) implicit none integer*4 np,ipl,i,j real*8 psi(np,ipl),dq,a write(6,*)' Wavefunction norms:' do 1 i=1,ipl a=0.0d0 do 2 j=1,np 2 a=a+psi(j,i)**2 a=a*dq 1 write(6,6000)i,a 6000 format(i5,f15.10) return end c ------------------------------------------------------------ subroutine doroa(t,e,p,np,psi,al,ipl,tol,dq,EXCNM) implicit none real*8 e(*),p(*),psi(np,ipl),al(np,3,3),tol, 1w,dq,aij(3,3),SAL0,SAL1,AMU,BOHR,SPAL,SpA3, 1beta2,EXCNM,YDY,YDX,ram1,P1,z,t integer*4 i,np,it,j,ix,ii,ipl,jx,JJ z=0.0d0 AMU=1822.0d0 BOHR=0.529177d0 open(18,file='ROA.TAB') WRITE(18,3304) 3304 FORMAT( 1'MODE FREQ Ramx (180) Ramy Ramtot CID90', 2' CID-X CID180 DEP ROA180 RamanDCPI',/, 3' cm-1 A^4/AMU x10^4', 4' x10^4 x10^4 A^5/AMUx1e4 au.10^4') write(18,3002) 3002 FORMAT(80(1H-)) it=0 do 1 i=1,ipl do 1 j=i+1,ipl w=e(j)-e(i) do 2 ix=1,3 do 2 jx=1,3 aij(ix,jx)=0.0d0 do 21 ii=1,np 21 aij(ix,jx)=aij(ix,jx)+psi(ii,i)*al(ii,ix,jx)*psi(ii,j) 2 aij(ix,jx)=aij(ix,jx)*dq SAL0=0.0d0 SAL1=0.0d0 SPAL=0.0d0 DO 3 II=1,3 SPAL=SPAL+aij(II,II) DO 3 JJ=1,3 SAL0=SAL0+aij(II,II)*aij(JJ,JJ) 3 SAL1=SAL1+aij(II,JJ)*aij(II,JJ) YDX=4.0d0*SAL1+2.0d0*SAL0 YDY=3.0d0*SAL1-SAL0 P1=YDY/YDX SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4)) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) YDY=6.0d0*YDY*(AMU*BOHR**4) YDX=6.0d0*YDX*(AMU*BOHR**4) c multiply by w*(1-exp(-hw/kT)) because of the drawing program c by 2.0 to get same results as gaussian (probably HO element) ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2)*p(i) 1*w/219470.0d0*2.0d0*(1.0d0-exp(-(1.44d0*w/t))) if(dabs(ram1).gt.tol)then it=it+1 WRITE(18,3000)it,w,YDX,YDY,ram1,z,z,z,P1,z,z,i,j 3000 FORMAT(I5,f9.2,6g12.3,f6.3,2g12.4,i4,' -> ',i4) endif 1 continue WRITE(18,1819)EXCNM 1819 FORMAT(60(1h-),/,f12.2,' nm') close(18) return end c ------------------------------------------------------------ subroutine dodog(e,p,np,psi,uo,ipl,tol,dq) implicit none real*8 e(*),p(*),uij(3),psi(np,ipl),uo(np,3),D,tol, 1pl,po,w,dq integer*4 i,np,it,j,ix,ii,ipl pl=4.0d0*atan(1.0d0) po=dsqrt(2.0d0*pl) open(18,file='DOG.TAB') WRITE (18,501) 501 FORMAT (/,' FREQ. ', 1' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2/,'-----------------------------------------------------------') it=0 do 1 i=1,ipl do 1 j=i+1,ipl w=e(j)-e(i) do 2 ix=1,3 uij(ix)=0.0d0 do 21 ii=1,np 21 uij(ix)=uij(ix)+psi(ii,i)*uo(ii,ix)*psi(ii,j) 2 uij(ix)=uij(ix)*dq D=(p(i)-p(j))*(uij(1)**2+uij(2)**2+uij(3)**2) if(dabs(D).gt.tol)then it=it+1 write(18,523)it,w,D,0.0d0,0.0d0,po*w*D,i,j endif 1 continue 523 FORMAT(I5,f12.2,4G15.6,i6,' -> ',i6) WRITE(18,1819) 1819 FORMAT(60(1h-)) close(18) return end c ------------------------------------------------------------ subroutine report(s) character*(*) s write(6,*)s stop end c ------------------------------------------------------------ subroutine getpar(t,ipl,ip,tol,EXCNM) real*8 t,tol,EXCNM integer*4 ipl,ip(*),ii character*4 key tol=1.0d-8 EXCNM=532.0d0 t=293.15d0 ipl=0 open(8,file='SURF.OPT') 1 read(8,800,end=88,err=88)key 800 format(a4) if(key.eq.'STAT')read(8,*)ipl,(ip(ii),ii=1,ipl) if(key.eq.'TEMP')read(8,*)t if(key.eq.'TOLE')read(8,*)tol if(key.eq.'EXCN')read(8,*)EXCNM goto 1 88 close(8) return end c ------------------------------------------------------------ subroutine getpsi(psi,np,ipl,ip) implicit none integer*4 np,ipl,ii,ist,ip(*),is,ix real*8 psi(np,ipl) character*5 s5 do 7 ii=1,ipl ist=ip(ii) write(s5,990)ist 990 format(i5) do 65 is=1,5 65 if(s5(is:is).ne.' ')goto 651 651 open(81,file=s5(is:5)//'.PRN') do 8 ix=1,np 8 read(81,*)psi(ix,ii),psi(ix,ii) 7 close(81) return end c ------------------------------------------------------------