program ch3 c toluen ch3 rotation IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) parameter (N0=20,nf0=150,nw0=10,nh0=120) DIMENSION wc(nw0),ws(nw0), 2APTu(nw0,N0*3,3),AATu(nw0,N0*3,3), 3ALu(nw0,N0*3,3,3),Gu(nw0,N0*3,3,3),Au(nw0,N0*3,3,3,3), 4W(3*N0),hr(nh0,nh0),qr(nh0),sr(nh0,nh0), 5EM(nh0),imr(nh0),mr(nh0),anr(nh0),hri(3*N0,nh0,nh0), 6sri(3*N0,nh0,nh0),qri(3*N0,nh0),ht(nh0,nh0),st(nh0,nh0), 7qt(nh0),wu(nh0,3*N0) character*9 s9 integer*4 state(3*N0) logical ldiag,ldia2,lbath parameter (nt0=100000) integer*4 s1(3*N0,nt0),s2(3*N0,nt0),i44(4,nt0) dimension et(nt0),rt(nt0) common/trb/s1,s2,et,rt,itb,i44,lbath Z0=0.0d0 pi=4.0d0*atan(1.0d0) open(9,file='SSD.OPT') read(9,*)ncos if(ncos.gt.nw0)call report('too many cosinus freqs') read(9,*)(wc(i),i=1,ncos) read(9,*)nsin if(nsin+ncos.gt.nw0)call report('too many cos and sinus freqs') read(9,*)(ws(i),i=1,nsin) read(9,*)fis,dfi fis=fis*pi/180.0d0 dfi=dfi*pi/180.0d0 read(9,*)nf if(nf.gt.nf0)call report('too many angles') fie=fis+dble(nf-1)*dfi write(6,6767)fis,fie,nf 6767 format(' Scan ',f8.1,' ... ',f8.1,',',i3,' points.') close(9) c load in polarizability normal mode Fourier coefficients do 103 in=1,ncos+nsin write(s9,499)in 499 format(i9) open(15,file=s9//'.u.ten') read(15,*) NAT if(NAT.gt.N0)call report('too many atoms') DO 106 L=1,NAT*3 106 read(15,*) (APTu(in,L,II),II=1,3) DO 2244 L=1,NAT*3 2244 read(15,*) (AATu(in,L,II),II=1,3) CLOSE(15) open(2,file=s9//'.u.ttt') read(2,*) read(2,*) read(2,*) read(2,*) DO 15 II=1,3 READ(2,*) DO 15 IND=1,NAT*3 15 READ(2,*)LDUM,IXDUM,(ALu(in,IND,II,J),J=1,3) READ(2,*) READ(2,*) DO 13 II=1,3 READ(2,*) DO 13 IND=1,NAT*3 13 READ(2,*)LDUM,IXDUM,(Gu(in,IND,II,J),J=1,3) READ(2,*) READ(2,*) DO 14 II=1,3 DO 14 J=1,3 READ(2,*) DO 14 IND=1,NAT*3 14 READ(2,*)LDUM,IXDUM,(Au(in,IND,II,J,K),K=1,3) close(2) an1=0.0d0 an2=0.0d0 do 151 IND=1,NAT*3 do 151 II=1,3 an1=an1+APTu(in,IND,II)**2 DO 151 J=1,3 151 an2=an2+ALu(in,IND,II,J)**2 if(in.le.ncos)then write(6,6009)in,wc(in),an1,an2 6009 format(i3,' cos(',f2.0,' x f), norm1=',f6.3,' norm2=',f6.3) else write(6,6019)in,ws(in-ncos),an1,an2 6019 format(i3,' sin(',f2.0,' x f), norm1=',f6.3,' norm2=',f6.3) endif 103 continue rc=5.537d0 v6=4.7d0 v3=0.0 mmax=5 T=298.0d0 Etol=50.0d0 plim=1.0d-3 nmin=0 nmax=0 c ldiag - diagonalize rotational Hamiltonian, c otherwise use free rotor ldiag=.false. c ldia2 - diagonalize each vibr state i - incorporate c dependence of harmonic frequency on rotation ldia2=.false. lbath=.false. open(9,file='CH3.OPT') 111 read(9,900,end=599,err=599)s9 900 format(a9) if(s9(1:4).eq.'ITEN')read(9,*)rc if(s9(1:4).eq.'MMAX')read(9,*)mmax if(s9(1:4).eq.'NMAX')read(9,*)nmax if(s9(1:4).eq.'NMIN')read(9,*)nmin if(s9(1:4).eq.'TEMP')read(9,*)T if(s9(1:4).eq.'ETOL')read(9,*)Etol if(s9(1:4).eq.'DIAG')read(9,*)ldiag if(s9(1:2).eq.'V6')read(9,*)v6 if(s9(1:2).eq.'V3')read(9,*)v3 if(s9(1:4).eq.'DIA2')read(9,*)ldia2 if(s9(1:4).eq.'BATH')read(9,*)lbath if(s9(1:4).eq.'PLIM')read(9,*)plim goto 111 599 close(9) c if(lbath)write(6,*)'Last six freqs at F.INP used as a bath.' c c first state: only cosinus, M=0 (rotational quantum number) i=1 anr(i)=1.0d0/dsqrt(2.0d0*pi) imr(i)=1 mr(i)=0 c c rest of the states: sinus and cosinus do 1041 m=1,mmax do 1041 im=1,2 i=i+1 if(i.gt.nh0)call report('Hr dimension too big') mr(i)=m imr(i)=im 1041 anr(i)=1.0d0/dsqrt(pi) c nr=i do 1042 i=1,nr do 1042 j=1,nr 1042 hr(i,j)=Z0 if(ldia2)then do 1043 i=1,nr do 1043 j=1,nr do 1043 iq=1,3*NAT 1043 hri(iq,i,j)=Z0 open(22,file='u.finp') read(22,*) do 109 iq=1,3*NAT 109 read(22,*)ldumm,(wu(in,iq),in=1,ncos+nsin) write(6,*)'u.finp read' close(22) endif C c get harmonic frequencies call READS(3*N0,W) do 1091 iq=1,3*NAT write(6,1099)iq,wu(1,iq),w(iq) 1099 format(i4,2f10.2) 1091 wu(1,iq)=w(iq) write(6,*)'Average replaced by frequencies from F.INP' do 104 i=1,nr an1=anr(i) m=mr(i) im=imr(i) do 104 j=1,nr an2=anr(j) mp=mr(j) imp=imr(j) if(ldia2)then do 1044 iq=1,3*NAT hri(iq,i,j)=hri(iq,i,j) 1 +0.5d0*v6*an1*an2*(ame(im,1,imp,m,0,mp)-ame(im,1,imp,m,6,mp)) 2 +0.5d0*v3*an1*an2*(ame(im,1,imp,m,0,mp)-ame(im,1,imp,m,3,mp)) 3 +rc* dble(mp)**2*ame(im,1,imp,m,0,mp)*an1*an2 do 1044 in=2,ncos+nsin if(in.le.ncos)then ip=nint(wc(in)) iu=1 else ip=nint(ws(in-ncos)) iu=2 endif 1044 hri(iq,i,j)=hri(iq,i,j)+wu(in,iq)*ame(im,iu,imp,m,ip,mp) endif c V ~ V6/2(1-cos(6fi))+RC.d^2/dfi^2 c sin, cos are not true eigenfuncction of the rotation, c but we neglect it now 104 hr(i,j)=hr(i,j) 1+0.5d0*v6*an1*an2*(ame(im,1,imp,m,0,mp)-ame(im,1,imp,m,6,mp)) 2+0.5d0*v3*an1*an2*(ame(im,1,imp,m,0,mp)-ame(im,1,imp,m,3,mp)) 3+rc* dble(mp)**2*ame(im,1,imp,m,0,mp)*an1*an2 do 108 i=1,nr 108 write(6,610)(hr(i,j),j=1,nr) 610 format(10f8.2) if(ldiag)then CALL TRED12(nr,hr,sr,qr,2,IERR,EM,nh0) if(IERR.ne.0)call report('Diagonalization error') write(6,*)'Rotational Energies:' write(6,610)(qr(i),i=1,nr) if(ldia2)then do 1047 iq=1,3*NAT write(6,*)' Mode ',iq do 1046 i=1,nr do 1046 j=1,nr 1046 ht(i,j)=hri(iq,i,j) CALL TRED12(nr,ht,st,qt,2,IERR,EM,nh0) if(IERR.ne.0)call report('Diagonalization error') write(6,*)'Rotational Energies:' write(6,610)(qt(i),i=1,nr) do 1047 i=1,nr qri(iq,i)=qt(i) do 1047 j=1,nr 1047 sri(iq,i,j)=st(i,j) endif else do 107 i=1,nr 107 qr(i)=hr(i,i) endif open(20,file='DOG.TAB') open(30,file='ROA.TAB') WRITE (20,501) 501 FORMAT (/,' FREQ. ', 1' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]') write(20,886) WRITE(30,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(30,886) C c c calculate partition function: EkT=T*8.31441d0/4184.0d0 cmtokcal=627.5d0/219470.0d0 Qo=0.0d0 do 26 ir=1,nr 26 Qo=Qo+exp(-qr(ir)*cmtokcal/EkT) Qv=0.0d0 c vibrational ground state: E=0.0d0 Qv=Qv+1.0d0 c first vibrational excitation do 21 i1=1,3*NAT if(w(i1).gt.Etol.or.(lbath.and.i1.gt.3*NAT-6))then E=E+w(i1) Qv=Qv+exp(-E*cmtokcal/EkT) c second vibrational excitation do 22 i2=i1,3*NAT if(w(i2).gt.Etol.or.(lbath.and.i2.gt.3*NAT-6))then E=E+w(i2) Qv=Qv+exp(-E*cmtokcal/EkT) c third vibrational excitation do 23 i3=i2,3*NAT if(w(i3).gt.Etol.or.(lbath.and.i3.gt.3*NAT-6))then E=E+w(i3) Qv=Qv+exp(-E*cmtokcal/EkT) c fourth vibrational excitation do 24 i4=i3,3*NAT if(w(i4).gt.Etol.or.(lbath.and.i4.gt.3*NAT-6))then E=E+w(i4) Qv=Qv+exp(-E*cmtokcal/EkT) c fifth vibrational excitation do 25 i5=i4,3*NAT if(w(i5).gt.Etol.or.(lbath.and.i5.gt.3*NAT-6))then E=E+w(i5) Qv=Qv+exp(-E*cmtokcal/EkT) E=E-w(i5) endif 25 continue E=E-w(i4) endif 24 continue E=E-w(i3) endif 23 continue E=E-w(i2) endif 22 continue E=E-w(i1) endif 21 continue Q=Qo*Qv write(6,6500)' approx. 1',Q,Qo,Qv 6500 format(' Partition functions, ',A10,': ',3f14.3) Qar=dsqrt(pi*EkT/rc/cmtokcal)*exp(-qr(i)*cmtokcal/EkT) Qav=1.0d0 do 3 iq=1,3*NAT if(w(iq).gt.Etol.or.(lbath.and.iq.gt.3*NAT-6))then Ekcal=w(iq)*cmtokcal Qav=Qav/(1.0d0-exp(-Ekcal/EkT)) endif 3 continue write(6,6500)' approx. 2',Qar*Qav,Qar,Qav c c generate all initial states, as above: c and calculate spectra iskip=0 itran=0 itb=0 E=0.0d0 Qc=0.0d0 do 4 iq=1,3*NAT 4 state(iq)=0 do 30 ir=1,nr m=mr(ir) im=imr(ir) an=anr(ir) call addqs(E,W,EkT,state,Q,Qc,m,im,APTu,AATu,wc,ws,plim, 1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri, 2nmin,nmax,0 ,mr,imr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic) if(ic.ne.0)goto 999 do 31 i1=1,3*NAT if(w(i1).gt.Etol.or.(lbath.and.i1.gt.3*NAT-6))then call addqs(E,W,EkT,state,Q,Qc,m,im,APTu,AATu,wc,ws,plim, 1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri, 2nmin,nmax,i1,mr,imr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic) if(ic.ne.0)goto 999 do 32 i2=i1,3*NAT if(w(i2).gt.Etol.or.(lbath.and.i2.gt.3*NAT-6))then call addqs(E,W,EkT,state,Q,Qc,m,im,APTu,AATu,wc,ws,plim, 1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri, 2nmin,nmax,i2,mr,imr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic) if(ic.ne.0)goto 999 do 33 i3=i2,3*NAT if(w(i3).gt.Etol.or.(lbath.and.i3.gt.3*NAT-6))then call addqs(E,W,EkT,state,Q,Qc,m,im,APTu,AATu,wc,ws,plim, 1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri, 2nmin,nmax,i3,mr,imr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic) if(ic.ne.0)goto 999 do 34 i4=i3,3*NAT if(w(i4).gt.Etol.or.(lbath.and.i4.gt.3*NAT-6))then call addqs(E,W,EkT,state,Q,Qc,m,im,APTu,AATu,wc,ws,plim, 1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri, 2nmin,nmax,i4,mr,imr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic) if(ic.ne.0)goto 999 do 35 i5=i4,3*NAT if(w(i5).gt.Etol.or.(lbath.and.i5.gt.3*NAT-6))then call addqs(E,W,EkT,state,Q,Qc,m,im,APTu,AATu,wc,ws,plim, 1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri, 2nmin,nmax,i5,mr,imr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic) if(ic.ne.0)goto 999 state(i5)=state(i5)-1 E=E-w(i5) endif 35 continue state(i4)=state(i4)-1 E=E-w(i4) endif 34 continue state(i3)=state(i3)-1 E=E-w(i3) endif 33 continue state(i2)=state(i2)-1 E=E-w(i2) endif 32 continue state(i1)=state(i1)-1 E=E-w(i1) endif 31 continue 30 continue 999 write(6,6500)' check',Qc write(20,886) write(30,886) 886 format(80(1h-)) close(20) close(30) write(6,889)iskip,itran,itb 889 format(' Transitions from ',i10,' high states skipped',/, 1I10,' transitions',I10,' recorded') write(6,*)'ordering' ntb=itb do 36 itb=1,ntb do 36 ip=itb+1,ntb if(et(ip).gt.et(itb))then t=et(ip) et(ip)=et(itb) et(itb)=t t=rt(ip) rt(ip)=rt(itb) rt(itb)=t do 37 iq=1,3*nat it=s1(iq,ip) s1(iq,ip)=s1(iq,itb) s1(iq,itb)=it it=s2(iq,ip) s2(iq,ip)=s2(iq,itb) 37 s2(iq,itb)=it do 39 ii=1,4 it=i44(ii,ip) i44(ii,ip)=i44(ii,itb) 39 i44(ii,itb)=it endif 36 continue open(9,file='T.LST') do 38 itb=1,ntb write(9,603)et(itb),rt(itb) 603 format(f8.2,f10.3,$) do 41 i=1,3*nat 41 if(s1(i,itb).gt.0)write(9,601)i,s1(i,itb) write(9,602) 602 format('->',$) do 40 i=1,3*nat 40 if((w(i).gt.Etol.or.(lbath.and.i.gt.3*NAT-6)) 1.and.s2(i,itb).gt.0)write(9,601)i,s2(i,itb) 601 format(i3,'_',i1,' ',$) 38 write(9,604)(i44(ii,itb),ii=1,4) 604 format(4i4) close(9) write(6,*)'Transition list T.LST written.' stop end SUBROUTINE REPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(6,*)RS WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') STOP END SUBROUTINE READS(N3,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION W(*) SFAC=0.0234280d0 OPEN(4,FILE='F.INP') read(4,*)nint,N NAT=N/3 N7=N-NINT+1 do 1 i=1,NAT 1 read(4,*) read(4,*) DO 2 I=1,NAT DO 2 J=1,NINT 2 read(4,*) read(4,*)nint READ(4,4000)(W(I),I=1,NINT) 4000 FORMAT(6F11.6) DO 7 I=NINT+1,N 7 W(I)=0.0d0 CLOSE(5) WRITE(6,*)' F.INP read' RETURN end subroutine addqs(E,W,EkT,state,Q,Qc,m,im,APTu,AATu,wc,ws,plim, 1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri, 2nmin,nmax,ia,mr,imr,qr,sr,nh0,anr,an1,ldiag,ir,nr,ic) c c add excitation and get all allowed transitions c implicit none integer*4 m,im,i,ii,imp,in,inp,ix,j,jj,jx,k,kx,l,mp,nat, 1ncos,nsin,state(*),nw0,iskip,itran,ip,nmin,nmax,ia,nh0, 1mr(*),imr(*),ir,iy,jy,nr,irp,nt0,N0,ic parameter (nt0=100000,N0=20) real*8 E,W(*),EkT,Q,P,Ekcal,APTt(3),AATt(3),wc(*),ws(*),Qc, 3ALt(3,3),Gt(3,3),At(3,3,3),APTu(nw0,N0*3,3),AATu(nw0,N0*3,3), 3ALu(nw0,N0*3,3,3),Gu(nw0,N0*3,3,3),Au(nw0,N0*3,3,3,3),ame,amu, 1beta2,betaa,betag,bohr,cid1,cid2,d90x,d90z,dx,dz, 1ep,eps,fc,fd,ff,fr,gpisvejc,PI,ram1,ram2,roa1, 1roa2,roa3,s90x,s90z,sa1,sag0,sag1,sal1,spa3,spal,sumd,wr, 1ydx,ydy,rdo,sal0,p1,Etol,Ep0,an1,qr(*),sr(nh0,nh0),anr(*), 1an2,ci,cj,Eg,sri(3*N0,nh0,nh0),qri(3*N0,nh0),plim logical ldiag,ldia2 integer*4 s1(3*N0,nt0),s2(3*N0,nt0),itb,iq,i44(4,nt0) real*8 et(nt0),rt(nt0) logical lbath common/trb/s1,s2,et,rt,itb,i44,lbath c ic=0 if(ia.gt.0)then state(ia)=state(ia)+1 E=E+W(ia) endif Eg=E+qr(ir) c pi=4.0d0*atan(1.0d0) Ekcal=Eg/219470.0d0*627.5d0 FD=dsqrt(388.9219d0) FR=2586.285D-8 FC=131.14D08 P=exp(-Ekcal/EkT)/Q Qc=Qc+exp(-Ekcal/EkT) c if((P.lt.plim.and..not.lbath).or.P.lt.1.0d-4)then iskip=iskip+1 return endif c write(6,600) 600 format(1h<,$) do 2 i=1,3*nat 2 if(w(i).gt.Etol.or.(lbath.and.i.gt.3*NAT-6))write(6,601)state(i) 601 format(i1,$) write(6,602)m,im,Eg,P 602 format(1h|,1h<,2i1,1h|,f8.2,' cm-1, p = ',f6.3) c c look at excited states |state> --> |state*> c harmonic transitions only (n->n+1) do 1 i=1,3*nat if(w(i).gt.Etol.or.(lbath.and.i.gt.3*NAT-6))then if(i.ge.nmin.or.nmin.eq.0)then if(i.le.nmax.or.nmax.eq.0)then state(i)=state(i)+1 Ep0=E+w(i) do 10 irp=1,nr mp=mr(irp) an2=anr(irp) Ep=Ep0+qr(irp) if(ldia2)Ep=Ep0+qri(i,irp) imp=imr(irp) c c transition tensors do 4 ix=1,3 APTt(ix)=0.0d0 AATt(ix)=0.0d0 do 4 jx=1,3 ALt(ix,jx)=0.0d0 Gt(ix,jx)=0.0d0 do 4 kx=1,3 4 At(ix,jx,kx)=0.0d0 if(ldiag)then do 61 iy=1,nr ci=sr(iy,ir) if(ldia2)ci=sri(i,iy,ir) an1=anr(iy) do 61 jy=1,nr cj=ci*sr(jy,irp) if(ldia2)cj=ci*sri(i,jy,irp) an2=anr(jy) if(dabs(cj).gt.1.0d-6)then do 62 in=1,ncos+nsin if(in.le.ncos)then ip=nint(wc(in)) inp=1 else ip=nint(ws(in-ncos)) inp=2 endif ff=ame(imr(iy),inp,imr(jy),mr(iy),ip,mr(jy))*cj*an1*an2 do 62 ix=1,3 APTt(ix)=APTt(ix)+APTu(in,i,ix)*ff AATt(ix)=AATt(ix)+AATu(in,i,ix)*ff do 62 jx=1,3 ALt(ix,jx)=ALt(ix,jx)+ALu(in,i,ix,jx)*ff Gt(ix,jx)=Gt(ix,jx)+Gu(in,i,ix,jx)*ff do 62 kx=1,3 62 At(ix,jx,kx)=At(ix,jx,kx)+Au(in,i,ix,jx,kx)*ff endif 61 continue else do 6 in=1,ncos+nsin if(in.le.ncos)then ip=nint(wc(in)) inp=1 else ip=nint(ws(in-ncos)) inp=2 endif ff=ame(im,inp,imp,m,ip,mp)*an1*an2 do 6 ix=1,3 APTt(ix)=APTt(ix)+APTu(in,i,ix)*ff AATt(ix)=AATt(ix)+AATu(in,i,ix)*ff do 6 jx=1,3 ALt(ix,jx)=ALt(ix,jx)+ALu(in,i,ix,jx)*ff Gt(ix,jx)=Gt(ix,jx)+Gu(in,i,ix,jx)*ff do 6 kx=1,3 6 At(ix,jx,kx)=At(ix,jx,kx)+Au(in,i,ix,jx,kx)*ff endif fr=w(i) WR=dsqrt(dabs(fr)) do 7 ix=1,3 if(WR.gt.0.1d0)APTt(ix)=APTt(ix)/WR*FD 7 AATt(ix)=AATt(ix)*FC*WR SUMD=0.0d0 RDO=0.0d0 do 8 ix=1,3 SUMD=SUMD+APTt(ix)**2 8 RDO=RDO+APTt(ix)*AATt(ix) SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 SPAL=0.0d0 DO 9 II=1,3 SPAL=SPAL+ALt(II,II) DO 9 JJ=1,3 SAL0=SAL0+ALt(II,II)*ALt(JJ,JJ) SAL1=SAL1+ALt(II,JJ)*ALt(II,JJ) SAG0=SAG0+ALt(II,II)* Gt(JJ,JJ) SAG1=SAG1+ALt(II,JJ)* Gt(II,JJ) DO 9 K=1,3 DO 9 L=1,3 9 SA1=SA1+EPS(II,K,L)*ALt(II,JJ)*At(K,L,JJ)/3.0d0 S90X=7.0d0*SAL1+SAL0 D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1) DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c IL+/-IR, 90o, z S90Z=2.0d0*(3.0d0*SAL1-SAL0) D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1) DZ=D90Z IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z c polarized backscattering components (YDY+YDX=S180): YDY=3.0d0*SAL1-SAL0 YDX=4.0d0*SAL1+2.0d0*SAL0 P1=YDY/YDX AMU=1822.0d0 BOHR=0.529177d0 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) gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/4880.0d0 c alphag=SAG0/9.0d0*gpisvejc betaa=SA1*3.0d0/2.0d0*gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc c DCP 180 c ram3=24.0d0*beta2 roa3=8.0d0*(12.0d0*betag+4.0d0*beta2) S90X=S90X*(AMU*BOHR**4) D90X=D90x*gpisvejc DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c ICP(90)=90Z roa2=8.0d0*(3.0d0*betag-betaa) ram2=12.0d0*beta2 CID2=0.0d0 if(ram2.gt.1.0d-9)CID2=roa2/ram2 c ICP(180): roa1=96.0d0*(betag+betaa/3.0d0) ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) CID1=0.0d0 if(ram1.gt.1.0d-9)CID1=roa1/ram1 ram1=ram1*P roa1=roa1*P roa3=roa3*P YDX=YDX*P YDY=YDY*P if(ram1.gt.0.00001d0.or.SUMD.gt.1.0d-6)then itran=itran+1 WRITE(30,3300)itran,Ep-Eg,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,roa3 3300 FORMAT(I8,f9.2,6f13.4,f6.3,2f11.2) SUMD=SUMD*P RDO=RDO*P write(20,523)itran,Ep-Eg,SUMD,RDO,RDO,dsqrt(2.0d0*pi)*fr*SUMD 523 format(i4,5G15.6) itb=itb+1 if(itb.gt.nt0)then write(6,*)'Too many transitions' ic=1 return endif do 11 iq=1,3*nat s1(iq,itb)=state(iq) 11 s2(iq,itb)=state(iq) s1(i,itb)=s1(i,itb)-1 et(itb)=Ep-Eg rt(itb)=ram1 i44(1,itb)=m i44(2,itb)=im i44(3,itb)=mp i44(4,itb)=imp endif 10 continue state(i)=state(i)-1 endif endif endif 1 continue return end function ame(ii,jj,kk,i1,j1,k1) c int(cos/sin(if)*cos/sin(jf)*cos/sin(kf),f=0..2pi) implicit none real*8 ame,pi,sum integer*4 ii,jj,kk,i1,j1,k1 sum=0.0d0 pi=4.0d0*atan(1.0d0) if(ii.eq.1.and.jj.eq.1.and.kk.eq.2)goto 99 if(ii.eq.1.and.jj.eq.2.and.kk.eq.1)goto 99 if(ii.eq.2.and.jj.eq.1.and.kk.eq.1)goto 99 if(ii.eq.2.and.jj.eq.2.and.kk.eq.2)goto 99 if(ii.eq.1.and.jj.eq.1.and.kk.eq.1)then if(i1+j1+k1.eq.0)sum=sum+pi/2.0d0 if(i1+j1-k1.eq.0)sum=sum+pi/2.0d0 if(i1-j1+k1.eq.0)sum=sum+pi/2.0d0 if(i1-j1-k1.eq.0)sum=sum+pi/2.0d0 goto 99 endif if(ii.eq.1.and.jj.eq.2.and.kk.eq.2)then if(i1+j1+k1.eq.0)sum=sum-pi/2.0d0 if(i1+j1-k1.eq.0)sum=sum+pi/2.0d0 if(i1-j1+k1.eq.0)sum=sum+pi/2.0d0 if(i1-j1-k1.eq.0)sum=sum-pi/2.0d0 goto 99 endif if(ii.eq.2.and.jj.eq.1.and.kk.eq.2)then if(i1+j1+k1.eq.0)sum=sum-pi/2.0d0 if(i1+j1-k1.eq.0)sum=sum+pi/2.0d0 if(i1-j1+k1.eq.0)sum=sum-pi/2.0d0 if(i1-j1-k1.eq.0)sum=sum+pi/2.0d0 goto 99 endif if(ii.eq.2.and.jj.eq.2.and.kk.eq.1)then if(i1+j1+k1.eq.0)sum=sum-pi/2.0d0 if(i1+j1-k1.eq.0)sum=sum-pi/2.0d0 if(i1-j1+k1.eq.0)sum=sum+pi/2.0d0 if(i1-j1-k1.eq.0)sum=sum+pi/2.0d0 goto 99 endif write(6,*)ii,jj,kk,i1,j1,k1 call report('ame-unknown combination') 99 ame=sum return end FUNCTION EPS(I,J,K) REAL*8 EPS 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 TRED12(N,A,Z,D,IEIG,IERR,E,MX3) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E,MX3) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E,MX3) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(MX3,N),D(MX3),E(MX3) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END