program mroax write(6,600) 600 format(' Simulation of rotationally resolved magnetic circular', 1 /, ' dichroism and Raman optical activity of X2 molecules',/, 1 /, ' Petr Bour, 1992-2014 Chicago/Prague',/,/, 2 ' Input: NO.PAR',/, 2 ' MROA.OPT',/, 2 ' Output: d.prn',/, 1 ' r.prn',/, 1 ' ram.prn',/, 1 ' roa.prn',/, 1 ' R.TAB',/, 1 ' C.TAB',/) call Setparameters call Initfac call computespectrum end c ============================================================ subroutine Setparameters implicit none real*8 d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth,wexc,wii,gamm common/nopar/d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth real*8 na,hprime,cvl,r,mb,Debye,thresh,mn common/rpar/na,hprime,cvl,r,mb,Debye,thresh,mn integer*4 jmax,jmin,amap,i,n common/ipar/jmin,jmax,amap(3) real*8 wmin,wmax integer*4 npo common/opts/wmin,wmax,npo character*1 branch,fstate common/spar/branch,fstate integer*4 isame,vib1,vib2,nv0,ise parameter (nv0=30) real*8 I1,I2,gj,gje,g1,g2,ff,fodd,feven,Bv,Dv,fv,gjs,BBs common/x2c/Bv(nv0),Dv(nv0),fv(nv0),wexc,wii,gamm,BBs, 1I1,I2,gj,gje,gjs,g1,g2,ff,fodd,feven,isame,vib1,vib2,ise character*10 s10,s10_a integer*4 iee logical lwr,lzmat,lno,lno2,ldc,lx2 common/lpar/lwr,lzmat,lno,lno2,ldc,lx2,iee character*1 abc(3),xyz(3) data abc/'a','b','c'/ data xyz/'x','y','z'/ na =6.022045d23 c {mol^-1} {Avogadro} hprime=1.05457266d-34 c {J*s} {Planck/2pi} cvl=2.99792458d8 c {m/s} {light velicity} r=8.31451d0 c {j/K/mol} {ideal gas const.} mb=0.46686437d0 c {cm-1/T} {Bohr magneton} mn=mb/1837.0d0 c {cm-1/T} {nuclear magneton} c m0=1.2566370d-6 c {Js^2C^-2m^-1} {vacuum permeativity} Debye=3.335445d-30 c {C.m} c threshold * maximal intesities recorded in no2 spectrum to NO2.TAB thresh=0.01d0 c map of axes x y z -> a b c amap(1)=1 amap(2)=2 amap(3)=3 c calculate degree of circularity: ldc=.false. c iee:0 no weighing of excited states in rrx2 c iee:1 Boltzmann weighing of excited states in rrx2 c iee:2 f-weighting of excited states in rrx2 c iee:3 f-weighting and SO coupling of excited states in rrx2 pi=4.0d0*datan(1.0d0) dsqrtpi=dsqrt(pi) jmax=10 jmin=0 t=298.0d0 bz=2.0d0 wexc=-10.0d0 bandwidth=0.3d0 inquire(file='NO.PAR',exist=lno) inquire(file='NO2.PAR',exist=lno2) inquire(file='X2.PAR',exist=lx2) if(lno.or.lx2)then if(lno)then if(lx2)call report('NO and X2 not compatible') open(9,file='NO.PAR') READ(9,*)jmax READ(9,*)jmin READ(9,*)branch READ(9,*)d0 READ(9,*)hw READ(9,*)t READ(9,*)t READ(9,*) READ(9,*)bcm READ(9,*)becm READ(9,*)acm READ(9,*)aecm READ(9,*)bz READ(9,*)l READ(9,*)fstate READ(9,*) READ(9,*)bandwidth else wii=-10.0d0 gamm=10.0d0 I1=0.0d0 I2=0.0d0 gJ=0.0d0 gje=0.0d0 gjs=0.0d0 g1=0.0d0 g2=0.0d0 ff=0.0d0 feven=0.0d0 fodd=0.0d0 isame=-1 BBs=0.0d0 n=0 do 3 i=1,nv0 3 fv(i)=1.0d0 open(9,file='X2.PAR') 10 read(9,90,end=22,err=22)s10 if(s10(1:2).eq.'Bv')then read(9,*)n if(n.gt.nv0)call report('Too many rotational constants') read(9,*)(Bv(i),i=1,n) endif if(s10(1:2).eq.'Dv')then read(9,*)n if(n.gt.nv0)call report('Too many rotational constants') read(9,*)(Dv(i),i=1,n) endif if(s10(1:2).eq.'fv')then read(9,*)n if(n.gt.nv0)call report('Too many scale factors') read(9,*)(fv(i),i=1,n) endif if(s10(1:2).eq.'IN' )READ(9,*)I1,I2 if(s10(1:2).eq.'gJ' )READ(9,*)gJ if(s10(1:3).eq.'gje' )READ(9,*)gje if(s10(1:3).eq.'gjs' )READ(9,*)gjs if(s10(1:4).eq.'g1g2')READ(9,*)g1,g2 if(s10(1:4).eq.'MULT')READ(9,*)ff if(s10(1:4).eq.'MUEV')READ(9,*)feven if(s10(1:4).eq.'MUOD')READ(9,*)fodd if(s10(1:4).eq.'SAME')READ(9,*)isame if(s10(1:3).eq.'WII') READ(9,*)wii if(s10(1:3).eq.'BBs') READ(9,*)BBs backspace 9 read(9,90)s10_a if(s10_a.eq.s10)call report(' X2.PAR: unknown option '//s10) goto 10 22 close(9) if(I1.eq.0.0d0)call report('I1 not defined') if(I2.eq.0.0d0)call report('I2 not defined') if(gJ.eq.0.0d0)call report('gJ not defined') if(gje.eq.0.0d0)call report('gJe not defined') if(gjs.eq.0.0d0)call report('gJs not defined') if(g1.eq.0.0d0)call report('g1 not defined') if(g2.eq.0.0d0)call report('g2 not defined') if(ff.eq.0.0d0)call report('ff not defined') if(feven.eq.0.0d0)call report('feven not defined') if(fodd.eq.0.0d0)call report('fodd not defined') if(isame.eq.-1)call report('isame not defined') if(wii.eq.-10.0d0)call report('wii not defined') if(BBs.eq.0.0d0)call report('BBs not defined') endif close(9) l=l/10.0d0 d0=d0*Debye y=acm/bcm ye=aecm/becm endif lwr=.false. npo=10001 lzmat=.true. wmin=0 wmax=4000 vib1=0 vib2=0 ise=0 open(9,file='MROA.OPT') 1 read(9,90,end=2,err=2)s10 90 format(a10) if(s10(1:4).eq.'NPOI')read(9,*)npo if(s10(1:4).eq.'WMIN')read(9,*)wmin if(s10(1:4).eq.'WMAX')read(9,*)wmax if(s10(1:3).eq.'LWR')read(9,*)lwr if(s10(1:4).eq.'ZMAT')read(9,*)lzmat if(s10(1:4).eq.'JMAX')read(9,*)jmax if(s10(1:4).eq.'JMIN')read(9,*)jmin if(s10(1:4).eq.'TEMP')read(9,*)t if(s10(1:4).eq.'LDOC')read(9,*)ldc if(s10(1:4).eq.'VIB1')read(9,*)vib1 if(s10(1:4).eq.'VIB2')read(9,*)vib2 if(s10(1:4).eq.'THRE')read(9,*)thresh if(s10(1:4).eq.'BAND')read(9,*)bandwidth if(s10(1:4).eq.'AMAP')read(9,*)amap(1),amap(2),amap(3) if(s10(1:2).eq.'BZ')read(9,*)bz if(s10(1:3).eq.'ISE')read(9,*)ise if(s10(1:3).eq.'IEE')read(9,*)iee if(s10(1:4).eq.'WEXC')read(9,*)wexc if(s10(1:4).eq.'GAMM')read(9,*)gamm backspace 9 read(9,90)s10_a if(s10_a.eq.s10)call report(' MROA.OPT: unknown option '//s10) goto 1 2 close(9) if(wexc.eq.-10.0d0)call report('wexc not defined') tcm=t*r/na/hprime/2.0d0/pi/cvl/100.d0 write(6,6000)npo,wmin,wmax,jmin,jmax,bandwidth,t,bz,lwr,ldc, 1xyz(1),abc(amap(1)),xyz(2),abc(amap(2)),xyz(3),abc(amap(3)) 6000 format(' Number of points ',i10,/, 1 ' wmin ',f10.2,/, 2 ' wmax ',f10.2,/, 2 ' jmin ',i10/, 2 ' jmax ',i10/, 2 ' bandwidth ',f10.2,/, 2 ' temperature (K) ',f10.2,/, 2 ' Bz (Tesla) ',f10.2,/, 3 ' LWR ',l10,/, 3 ' LDOC ',l10,/, 3 ' axis ',a1,' is ',a1,/, 3 ' axis ',a1,' is ',a1,/, 3 ' axis ',a1,' is ',a1,/) return end c ============================================================ subroutine computespectrum implicit none integer*4 i,j,k,j1,im,ix,jx,nm,jp,ir real*8 aj,mj,ajg,aje,bjg,bje,x,fg0,xe,tj,mjp,pr,c22, 1f,mf,fe,MFE,egroundj,eexcj,enp,enm,dipp,dipm,delw, 1egg,cgcop,cgcom,so,soa,rthrj,gfactor,rcg,sig,nt,c2, 1jj,tjp,tjtjp,ileft,iright,rr,ri,li,lr,gF,kk,gFE, 1rami,roai,clight, 1g(3,3),abc(3),ta,AA,BB,CC,AE,BE,CE,u(3),rtol, 1asci(3,3),ascr(3,3),ilr,ili,irr,iri real*8 z0,half,one,th,two,three data z0,half,one,th,two,three 1/0.0d0,0.5d0,1.0d0,1.5d0,2.0d0,3.0d0/ character*4 uni real*8 d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth real*8 wmin,wmax integer*4 npo common/opts/wmin,wmax,npo integer*4 jmax,jmin,amap common/ipar/jmin,jmax,amap(3) real*8,allocatable::d(:),rs(:),ram(:),roa(:) character*1 branch,fstate common/spar/branch,fstate real*8 na,hprime,cvl,r,mb,Debye,thresh,mn common/rpar/na,hprime,cvl,r,mb,Debye,thresh,mn integer*4 iee logical lwr,lzmat,lno,lno2,ldc,lx2 common/lpar/lwr,lzmat,lno,lno2,ldc,lx2,iee integer*4, allocatable::ntw(:) real*8, allocatable::alphad(:,:,:,:),Gpd(:,:,:,:),ad(:,:,:,:,:), 1tw(:,:),al(:,:,:),asi(:,:,:),asr(:,:,:),aii(:,:,:),air(:,:,:), 1azz(:),app(:) logical lram,lg,lsym inquire(file='FILE.Q.TTT',exist=lram) if(lram)then c normal mode polarizability derivatives c read number of transition frequencies only: open(15,file='FILE.Q.TTT') READ(15,*) READ(15,*)nm close(15) allocate(alphad(1,nm,3,3),gpd(1,nm,3,3),ad(1,nm,3,3,3), 1 tw(1,nm+1),al(nm,3,3),asr(nm,3,5),asi(nm,3,5),ntw(1+1), 2 aii(nm,3,3),air(nm,3,3),azz(nm),app(nm)) ntw(1)=nm c read derivatives/transition polarizabilities/, molecular frame CALL READTTTQ(nm,1,alphad,Gpd,ad,'FILE.Q.TTT',tw,1,ntw,nm) do 6 im=1,nm do 61 ix=1,3 do 61 jx=1,3 61 al(im,amap(ix),amap(jx))=alphad(1,im,ix,jx) azz(im)=al(im,3,3) app(im)=0.5d0*(al(im,1,1)+al(im,2,2)) 6 continue c transform to spherical components call cartosphere(al,asr,asi,nm) c transform to spherical components - but dipoles individually call cartospherei(al,air,aii,nm) write(6,6109)nm 6109 format(i3,' polarizability derivatives found') if(lwr)open(78,file='R.TAB') endif inquire(file='G.OUT',exist=lg) if(lg)then rtol=1.0d0-6 clight=2.99792458d10 ta=1.0d9/clight uni='CM-1' AA=z0 BB=z0 CC=z0 call lookatout(lzmat,g,abc,ta,AA,BB,CC,AE,BE,CE,u) lsym=dabs(AA-BB).lt.rtol write(6,*)u write(6,60000)AA,BB,CC,uni,AE,BE,CE,uni, 1 u(1),u(2),u(3),lsym,rtol 60000 format(' A B C ground :', 3f12.6,a4,/, 1 ' A B C excited :', 3f12.6,a4,/, 1 ' ux uy uz :',3f12.6,/, 1 ' Symmetric top :',l12,' (tol:',g12.6,')') endif if(lno2)then if(lno.or.lx2)call report('only one of no no2 x2 possible') call rrno2(Jmin,Jmax,lsym,AA,BB,CC,AE,BE,CE,nm,air,aii, 1 wmin,wmax,npo,tw) else if(lno)then call rrspec(Jmin,Jmax,lsym,AA,BB,CC,AE,BE,CE,nm,air,aii, 1 wmin,wmax,npo,tw) else call rrx2(Jmin,Jmax,nm,azz,app,wmin,wmax,npo,tw,iee) endif endif if(.not.lno)call report('NO.PAR not found, program finished') if(lwr)open(77,file='C.TAB') nt=dble(npo)/(wmax-wmin) delw=(wmax-wmin)/dble(npo-1) allocate(d(npo),rs(npo),ram(npo),roa(npo)) do 5 i=1,npo d(i)=z0 ram(i)=z0 roa(i)=z0 5 rs(i)=z0 c write(6,6009)jmin,jmax,fstate,branch,bz 6009 format(' jmin: ',i5,' jmax: ',i5,' fstate: ',4x,a1, 1 ' branch: ',4x,a1,/, 1 ' field : ',f10.3) i=0 ir=0 do 1 j=jmin , jmax c J=1/2,3/2,...: aj=half+j tj=two*aj+one do 3 k=1 , 2 kk=dble(k) c ********************************************** c * * c * k=1 ... ground state ... approx 1/2 * c * * c * k=2 ... higher state ... approx 3/2 * c * * c ********************************************** if ((fstate.ne.'m').or.(k.ne.2))then if ((fstate.ne.'p').or.(k.ne.1))then call setxa(k,x,ajg,bjg,aj,y) egroundj=bcm *((aj-half)*(aj+1.5d0)+(dble(k)-1.5d0)*x +one) if(lram)then c RR MROA for NO (nitric oxide): write(6,*)j,k do 10 jp=-2,2 c loop over final rotational states J' = J-2, J,J+2: jj=aj+dble(jp) if(jj.ge.z0) then tjp=two*jj+one tjtjp=tj*tjp call setxa(k,xe,aje,bje,jj,ye) eexcj=becm*((jj-half)*(jj+1.5d0)+(kk-1.5d0)*xe+one) F=dabs(aj-one)-one c F = J-1,J,J+1: 502 F=F+one gF=gfactor(F,aj,y,kk)*bz*mb MF=-F-one c MF = -F .... F 501 MF=MF+one c energy of state |F MF>: c |E MF> = |JE MJ>: egg=egroundj+gF*MF pr=exp(-egg/tcm) fe=dabs(jj-one)-one c FE = J'-1 ... J'+1 503 fe=fe+one gFE=gfactor(FE,jj,ye,kk)*bz*mb MFE=mf-three c MFE = MF-2 ... MF+2, |MFE|<=FE 504 MFE=MFE+one if(dabs(MFE).le.fe)then c energy of state |FE MFE>: do 800 im=1,nm enm=tw(1,im)+eexcj+gFE*MFE-egg do 81 ix=1,3 do 81 jx=1,3 asci(ix,jx)=aii(im,ix,jx) 81 ascr(ix,jx)=air(im,ix,jx) ilr =z0 ili =z0 irr=z0 iri =z0 mj=mf-two c MJ = MF-1 ... MF+1, |MJ| <= J 505 mj=mj+one if(dabs(mj).le.aj)then c MJ' = MFE-1 ... MFE+1, |MJ'| <= J', because of nucl spin, c MJ'=MFE-MF+MJ mjp=MFE-MF+mj c2=rcg(jj,mjp,one,MF-mj,FE,MFE)*rcg(aj, mj,one,MF-mj,F , MF) c22=c2*ajg*aje c 1/2 1/2 + call ux(ri,rr,lr,li,aj,half,mj,jj,half,mjp,ascr,asci) call incr(ilr,ili,irr,iri,lr,li,rr,ri,c22) c 1/2 1/2 - call ux(ri,rr,lr,li,aj,-half,mj,jj,-half,mjp,ascr,asci) call incr(ilr,ili,irr,iri,lr,li,rr,ri,c22) c22=c2*bjg*bje c 3/2 3/2 + call ux(ri,rr,lr,li,aj,th,mj,jj,th,mjp,ascr,asci) call incr(ilr,ili,irr,iri,lr,li,rr,ri,c22) c 3/2 3/2 - call ux(ri,rr,lr,li,aj,-th,mj,jj,-th,mjp,ascr,asci) call incr(ilr,ili,irr,iri,lr,li,rr,ri,c22) endif if(mj.lt.mf+one)goto 505 ileft =ilr**2+ili**2 iright=irr**2+iri**2 rami=(iright+ileft)*pr*tjtjp roai=(iright-ileft)*pr*tjtjp if(rami.gt.1.0d-9)then ir=ir+1 if(lwr)write(78,770)ir,enm,rami,roai,aj,k, 1 f,mf,jj,fe,MFE call spread1(enm,rami,ram,npo,wmin,nt,delw) call spread1(enm,roai,roa,npo,wmin,nt,delw) endif 800 continue endif c MFE0 10 continue endif do 4 j1=-1 , 1 if (branch.ne.'q'.or.j1.eq.0) then if (branch.ne.'p'.or.j1.eq.-1) then if (branch.ne.'r'.or.j1.eq. 1) then jj=aj+dble(j1) if(jj.ge.z0)then call setxa(k,xe,aje,bje,jj,ye) eexcj=hw+becm*((jj-half)*(jj+th)+(kk-th)*xe+one) fg0=(ajg*aje*rthrj(jj,one,aj,-half,z0, half) 1 - bjg*bje*rthrj(jj,one,aj,-th,z0, th)) 2 *dsqrt((two*aj+one)*(two*jj+one)) f= dabs(aj-one)-one 999 f=f+one gF=gfactor(F,aj,y,kk)*bz*mb fe=abs(jj-one)-one 998 fe=fe+one gFE=gfactor(FE,jj,ye,kk)*bz*mb mf=-f-one 997 mf=mf+one c {LLLLLLLLcLLeft polarized light: } c { F,MF (j,k,p) -> Fe,MFe (j,k,p) } i=i+1 MFE=mf+one dipp=z0 egg=egroundj-gf*mf enp=eexcj-gfe*MFE-egg so=enp mj=-aj-one 801 mj=mj+one cgcop=rcg(aj,mj,one,mf-mj,f,mf)* 1 rcg(jj,mj+one,one,MFE-(mj+one),fe,MFE) dipp=dipp+rthrj(jj,one,aj,-(mj+one),one,mj) 1 *sig(jj-(mj+one))*cgcop if(NINT(2*mj).ne.NINT(2*aj) )goto 801 soa=enp*dipp**2*exp(-egg/tcm)*fg0**2 if(dabs(soa).lt.1.0d-10)then i=i-1 else call spread(one,so,soa,d,rs,npo,wmin,nt,delw) if(lwr)write(77,770)i,enp,soa/enp/two,soa/enp,aj,k, 1 f,mf,jj,fe,MFE 770 format(i8,f10.3,2f12.5,f5.1,'(',i1,')', 1 2f5.1,' -> ',3f5.1) endif c < G | u+ | E > absorption i=i+1 MFE=mf-one egg=egroundj-gfactor(F,aj,y,dble(k))*mf*bz*mb enm=eexcj-gfactor(fe,jj,ye,dble(k))*MFE*bz*mb-egg so=enm dipm=0 mj=-aj-one 701 mj=mj+one cgcom=rcg(aj,mj,one,mf-mj,f,mf) 1 *rcg(jj,mj-one,one,MFE-(mj-one),fe,MFE) dipm=dipm+rthrj(jj,one,aj,-(mj-one), -one,mj) 1 *sig(jj-(mj-one))*cgcom if( NINT(2*mj).ne.NINT(2*aj) )goto 701 soa=enm*dipm**2*exp(-egg/tcm)*fg0**2 if(dabs(soa).lt.1.0d-10)then i=i-1 else call spread(-one,so,soa,d,rs,npo,wmin,nt,delw) if(lwr)write(77,770)i,enm,dabs(soa)/enm/two, 1 -dabs(soa)/enm,aj,k,f,mf,jj,fe,MFE endif if(NINT(2*mf).ne.NINT(2*f) )goto 997 if(NINT(2*fe).ne.NINT(2*abs(jj+one)))goto 998 if(NINT(2*f ).ne.NINT(2*abs(aj+one)))goto 999 endif endif endif endif 4 continue c j1 endif c p endif c m 3 continue c k 1 continue c j if(lwr)then write(77,771) 771 format('---- ') close(77) if(lram)then write(78,771) close(78) endif endif call wrs('d.prn',d,wmin,wmax,npo) call wrs('r.prn',rs,wmin,wmax,npo) if(lram)then call wrs('ram.prn',ram,wmin,wmax,npo) call wrs('roa.prn',roa,wmin,wmax,npo) endif write(6,6000)i,npo 6000 format(i10,' transitions spread over ',i10,' spectral points') stop end c ============================================================ subroutine rrx2(Jmin,Jmax,nm,azz,app,wmin,wmax,npo,tw,iee) implicit none integer*4 Jmin,Jmax,J,JE,id,nm,im,I,iv,ije,npo,F,M,ME,MF,MFE,FE, 1ndim,isame,vib1,vib2,nv0,idx,ise,Js,iee real*8 Te,qp,aj,aje,amf,amfe,af,ai,gf,gfe,J2,JE2,F2,FE2,z0, 1azz(*),app(*),eee,nt,tf,tw(1,nm),a1z,a1p, 1egg,ileft,iright,rcg16,cg2,rami,roai,p,delw,wmax,wmin,ilr,roamax, 1irr,rr,lr,gi,afe,am,ame,ai2,mu,evib,CID,rammax,BBs,drr,dlr, 1I1,I2,gj,gje,g1,g2,ff,fodd,feven,Bv,Dv,DD,fv,d0,hw,t,BB,BE,acm, 1aecm,bz,l,tcm,dsqrtpi,pi,y,ye,bandwidth,gjs,ggg,wexc,wii,gamm, 1dl,dr,doci real*8,allocatable::ram(:),roa(:),rrJ(:),lrJ(:),DOCl(:),DOCr(:), 1cg2s(:),docs(:) parameter (nv0=30) common/x2c/Bv(nv0),Dv(nv0),fv(nv0),wexc,wii,gamm,BBs, 1I1,I2,gj,gje,gjs,g1,g2,ff,fodd,feven,isame,vib1,vib2,ise data z0/0.0d0/ common/nopar/d0,hw,t,BB,BE,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth real*8 na,hprime,cvl,r,mb,Debye,thresh,mn common/rpar/na,hprime,cvl,r,mb,Debye,thresh,mn logical allowed,bosons,Jeven,Ieven allocate(ram(npo),roa(npo),docs(npo)) write(6,611) 611 format(/' Simulation of MROA spectrum of a X2 molecule',/,/, 1 ' Input parameters',/,/, 2' v(cm) Bv(10-2cm-1) Dv(10-8cm-1) sf:',/) write(6,609)0,0.0d0,1.0d2*Bv(1),1.0D8*Dv(1) do 10 im=1,nm 10 write(6,609)im,tw(1,im),1.0d2*Bv(im+1),1.0d8*Dv(im+1),fv(im) 609 format(i4,f10.2,3f10.5) write(6,610)g1,g2,gJ,gje,gjs,Jmin,Jmax,I1,I2,npo,t,wmin,wmax, 1bz,mn,bandwidth,ff,fodd,feven,isame,iee,wii,wexc,gamm, 2vib1,vib2 610 format(/,' g1 ',f10.3,' g2 ',f10.3,/, 1 ' gJ ',f10.3,' gJe ',f10.3,/, 1 ' g* ',f10.3,/, 1 ' Jmin ',i10 ,' Jmax ',i10 ,/, 1 ' I1 ',f10.3,' I1 ',f10.3,/, 1 ' npoints ',i10 ,' T [K] ',f10.3,/, 1 ' wmin [cm-1] ',f10.3,' wmax [cm-] ',f10.3,/, 1 ' Bz [tesla] ',f10.3,' uN[cm-1 tesla-1]',f10.4,/, 1 ' bandwidth [cm-1]',f10.4,' mult. factor ',f10.4,/, 1 ' mult. fac. Jodd ',f10.4,' mult. fac. Jeven',f10.4,/, 1 ' same(0) nuclei ',i10 ,' IEE exc',i10,/, 1 ' w_res ',f10.3,/, 1 ' w_exc (cm-1) ',f10.3,' Gamma (cm-1) ',f10.3,/, 1 ' vib1 ',i10 ,' vib2 ',i10/,/) if(isame.eq.1)then write(6,*)'different nuclei' else bosons=mod(nint(I1*2.0d0),2).eq.0 if(bosons)then write(6,*)'same nuclei - bosons' else write(6,*)'same nuclei - fermions' endif endif if(ise.eq.0.or.isame.eq.1)then write(6,*)'all J* summed over' else write(6,*)'only allowed J* summed over' endif if(iee.eq.0)write(6,*)'non-resonance ux0' if(iee.eq.1)write(6,*)'resonance Boltzmann ux0' if(iee.eq.2)write(6,*)'resonance f ux0' if(iee.eq.3)write(6,*)'resonance f-SO ux0' if(iee.eq.4)write(6,*)'resonance f-SO ux0, Sz = 0 only' write(6,*) do 9 id=1,npo ram(id)=z0 docs(id)=z0 9 roa(id)=z0 gi=(g1+g2)/2.0d0 Te=t c kelvin to cm-1 tf=0.695d0 delw=(wmax-wmin)/dble(npo-1) nt=dble(npo)/(wmax-wmin) c approximate partition function: qp=0 do 101 iv=0,nm if(iv.eq.0)then evib=0.0d0 else evib=tw(1,iv) endif 101 qp=qp+exp(-evib/(Te*tf)) qp=qp*dsqrt((Te*tf)/Bv(1)) c vibrational ground state do 1101 iv=vib1,vib2 if(iv.eq.0)then evib=0.0d0 else evib=tw(1,iv) endif BB=Bv(iv+1) DD=Dv(iv+1) c rotational quantum numbers do 1101 J=Jmin,Jmax Jeven=mod(J,2).eq.0 if(Jeven)then mu=feven*ff else mu=fodd*ff endif aj=dble(J) J2=aj*(aj+1.0d0) c precalculate for each J ux0 elements c ux0(M=-J...J, Je=J-2,J,2, Me = M-2...M+2 * nvib * nJ* c (2J+1)* 3 * 5 *(nm-iv) *3 c J* = |J-1| ... J+1 ndim=(2*J+1)*15*(nm-iv)*3 allocate(rrJ(ndim),lrJ(ndim),cg2s(2*J+1),DOCl(ndim),DOCr(ndim)) do 1102 im=iv+1,nm a1z=azz(im-iv)*fv(im-iv) a1p=app(im-iv)*fv(im-iv) do 1102 M=-J,J am=dble(M) do 1102 JE=J-2,J+2,2 if(JE.ge.0)then ije=(JE-J+2+1)/2+1 aje=dble(JE) do 11021 ME=M-2,M+2 ame=dble(ME) if(iee.eq.0)then call ux0(rr,lr,aj,am,aje,ame,a1z,a1p,Jeven,isame,ise) else if(iee.eq.1)then do 11022 Js=abs(J-1),J+1 call ux1(rr,lr,aj,am,aje,ame,a1p,tf,te,qp,bz,gjs,Js,JE) idx=ME-M+2+1+5*(ije-1)+15*(M+J)+15*(2*J+1)*(im-iv-1) 1 +15*(2*J+1)*(im-iv-1)*(Js-abs(J-1)) rrJ(idx)=rr lrJ(idx)=lr 11022 continue else c approximate ground-state energy: egg=evib+J2*(BB-DD*J2) if(iee.eq.2)then call ux2(rr,lr,aj,am,aje,ame,a1p,wexc, 1 wii,gjs,mn,mb,bz,egg,evib,BBs,DD,gamm) else if(iee.eq.3)then call ux3(rr,lr,aj,am,aje,ame,a1p,wexc, 1 wii,gjs,mn,mb,bz,egg,evib,BBs,DD,gamm) else if(iee.eq.4)then call ux4(dl,dr,rr,lr,aj,am,aje,ame,a1p,wexc, 1 wii,gjs,mn,mb,bz,egg,evib,BBs,DD,gamm) else call report('unknown type of IEE') endif endif endif endif endif if(iee.ne.1)then idx=ME-M+2+1+5*(ije-1)+15*(M+J)+15*(2*J+1)*(im-iv-1) rrJ(idx)=rr lrJ(idx)=lr DOCl(idx)=dl DOCr(idx)=dr endif 11021 continue endif 1102 continue c total nuclear spin do 100 I=nint(abs(I1-I2)),nint(I1+I2) ai=dabs(dble(i)) ai2=ai*(ai+1.0d0) Ieven=mod(I,2).eq.0 if(isame.eq.1)then allowed=.true. else if(bosons)then allowed=(Jeven.and.Ieven).or.((.not.Jeven).and.(.not.Ieven)) else allowed=((.not.Jeven).and.Ieven).or.(Jeven.and.(.not.Ieven)) endif endif if(allowed)then c total momentum: do 1001 F=abs(J-I),J+I if(F.ge.0)then af=dble(F) F2=dble(F*(F+1)) c ground state g-factor x uN x B gf=ggg(gJ,gI,F2,J2,ai2)*mn*bz do 1501 MF=-F,F amf=dble(MF) c ground state energy egg=evib+J2*(BB-DD*J2)-gf*amf c Boltzmann factor: p=exp(-egg/tf/Te)/qp if(MF.eq.-F)write(6,6005)iv,F,J,I,egg,p 6005 format(' |v F J I >: ',4i4,' e:',f10.4,' cm-1 p:',f14.9) c excited state rotational quantum number: do 1501 JE=J-2,J+2,2 if(JE.ge.0)then ije=(JE-J+2+1)/2+1 aje=dble(JE) JE2=aje*(aje+1.0d0) c total momentum (IE=I): do 10011 FE=abs(JE-I),JE+I if(FE.ge.0)then afe=dble(FE) FE2=afe*(afe+1.0d0) c excited state g-factor x mn x bz gfe=ggg(gje,gI,FE2,JE2,ai2)*mn*bz do 10012 MFE=max(-FE,MF-2),min(FE,MF+2) amfe=dble(MFE) c pre-calculate cg2s: do 1103 M=-J,J am=dble(M) ME=MFE-MF+M if(abs(ME).le.JE)then ame=dble(ME) cg2s(M+J+1)= rcg16(aj ,am ,ai ,amf -am ,af ,amf ) 1 * rcg16(aje,ame,ai ,amfe-ame,afe,amfe) endif 1103 continue do 8 im=iv+1,nm eee=tw(1,im)+JE2*(Bv(im+1)-Dv(im+1)*JE2)-gfe*amfe-egg c Stokes scattering: if(eee.gt.0.0d0)then c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee if(iee.eq.1)then do 11024 Js=abs(J-1),J+1 ilr=0.0d0 irr=0.0d0 dlr=0.0d0 drr=0.0d0 do 10021 M=-J,J am=dble(M) c Clebsh-Gordon ME=MFE-MF+M if(abs(ME).le.JE)then c Clebsh-Gordon cg2=cg2s(M+J+1) idx=ME-M+2+1+5*(ije-1)+15*(M+J)+15*(2*J+1)*(im-iv-1) 1 +15*(2*J+1)*(im-iv-1)*(Js-abs(J-1)) ilr =ilr+cg2*lrJ(idx) irr =irr+cg2*rrJ(idx) endif c |ME| <= JE 10021 continue ileft =J2*aje*ilr**2 iright=J2*aje*irr**2 rami=(iright+ileft)*p*mu roai=(iright-ileft)*p*mu if(rami.gt.1.0d-19)then call spread1(eee,rami,ram,npo,wmin,nt,delw) call spread1(eee,roai,roa,npo,wmin,nt,delw) endif 11024 continue else c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee ilr=0.0d0 irr=0.0d0 dlr=0.0d0 drr=0.0d0 do 1002 M=-J,J am=dble(M) c Clebsh-Gordon ME=MFE-MF+M if(abs(ME).le.JE)then c Clebsh-Gordon cg2=cg2s(M+J+1) idx=ME-M+2+1+5*(ije-1)+15*(M+J)+15*(2*J+1)*(im-iv-1) c ROA: ilr =ilr+cg2*lrJ(idx) irr =irr+cg2*rrJ(idx) c DOC: dlr =dlr+cg2*DOCl(idx) drr =drr+cg2*DOCr(idx) c write(6,*)idx,cg2,riJ(idx) endif c |ME| <= JE 1002 continue ileft =J2*aje*ilr**2 iright=J2*aje*irr**2 rami=(iright+ileft)*p*mu roai=(iright-ileft)*p*mu if(rami.gt.1.0d-19)then call spread1(eee,rami,ram,npo,wmin,nt,delw) call spread1(eee,roai,roa,npo,wmin,nt,delw) endif ileft =J2*aje*dlr**2 iright=J2*aje*drr**2 doci=(iright-ileft)*p*mu if(dabs(doci).gt.1.0d-19) 1 call spread1(eee,doci,docs,npo,wmin,nt,delw) endif c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee endif c Stokes 8 continue 10012 continue c MFE endif c FE>0 10011 continue c FE endif c JE>=0 1501 continue c JE=J-2,J+2,2 c MF=-F...F endif c F>0 1001 continue endif c if((mod(J,2).eq.0.and.mod(I,2).ne.0).or. c 1 (mod(J,2).ne.0.and.mod(I,2).eq.0))then 100 continue deallocate(rrJ,lrJ,cg2s,DOCr,DOCl) 1101 continue call wrs('ramJ.prn',ram ,wmin,wmax,npo) call wrs('docJ.prn',docs,wmin,wmax,npo) call wrs('roaJ.prn',roa ,wmin,wmax,npo) CID=0.0d0 rammax=0.0d0 roamax=0.0d0 do 901 i=1,npo if( ram(i) .gt. rammax )rammax=ram(i) 901 if(dabs(roa(i)).gt.dabs(roamax))roamax=roa(i) if(rammax.ne.0.0d0)CID=roamax/rammax write(6,6901)CID*1.0d4 6901 format(' Largest CID x 10^4: ',f12.6) return end c ============================================================ subroutine ux4(dl,dr,srr,slr,J,M,JP,MP,app, 1wexc,wii,gjs,mn,mb,bz,egg,evib,BB,DD,gamm) c uxl/r= c c sum(J*,K*) (-1)^(K*+M+M') (2J*+1) [( J* 1 J) - (J* 1 J ) ] c (M'+/-1 -1 -M) (M'+/-1 1 -M ) c c (J' 1 J*) sum(q q') (J* 1 J ) ( J' 1 J*) al(q,q') c (M' +/-1 -M'-/+1) (K* q 0 ) ( 0 qp -K*) c c (wii^2-wexc^2) * wii c * ------------------------------------ c (wii^2 - wexc^2)^2+4*wii^2*Gamma^2 c c if ise=1 and isame=1, sum only over allowed Js=J=J c as ux3, but only Sz =0 allowed implicit none real*8 J,M,JP,MP,Js,js2,sumr,z0,one,two,wexc,wii,gamm, 1app,slr,srr,rthrj16,r1m,r1p,r2m,r2p,wbb,s2,Fs,fs2, 2es,gs,ff,ggg,gel,gjs,mn,mb,bz,egg,evib,BB,DD,ess,gf,MFs, 3Ms,rcg16,ip,im,op,om,pf,dr,dl data z0,one,two/0.0d0,1.0d0,2.0d0/ c g-factor for electron gel=-2.00231930436153d0 s2=one*(one+1) srr=z0 slr=z0 dr=z0 dl=z0 Js=J-two c J* = J-1,J,J+1 1 Js=Js+one if(Js.ge.z0.and.dabs(Js-JP).le.one)then js2=Js*(Js+one) r1m=rthrj16(Js,one, J,-one, one, z0) r1p=rthrj16(Js,one, J, one,-one, z0) r2m=rthrj16(JP,one,Js, z0,-one, one) r2p=rthrj16(JP,one,Js, z0, one,-one) c [(J* 1 J)*(JP 1 J*) + (J* 1 J)*(JP 1 J*)] * (alpha_xx+alpha_yy)/2 c [(-1 1 0) ( 0 -1 1 ) ( 1 -1 0) ( 0 1 -1 )] sumr=(r1m*r2m+r1p*r2p)*app c approximate g-factor of the excited(J*) state: gs=(gjs*mn*(js2-one)-mb)/js2 c approximate energy of the excited(*) state: ess=wii+evib+(js2-one)*(BB-DD*js2) c construct excited states F* = | J* - S | ... | J* + S|, S=1 Fs=Js-two 2 Fs=Fs+one if(Fs.ge.z0.and.dabs(Js-JP).le.one)then fs2=Fs*(Fs+one) c approximate g-factor of the excited(F*) state: gf=ggg(gs,gel*mb,fs2,js2,s2)*bz c MFS = M-2 ... M+2 MFs=M-two-one 3 MFs=MFs+1 if(abs(MFs).le.Fs.and.abs(MFs-MP).eq.one)then Ms=MFs es=ess-MFs*gf c ff=(es-egg)/(es-egg+wexc)/(es-egg-wexc) wbb=(es-egg+wexc)*(es-egg-wexc) ff=wbb*(es-egg)/(wbb**2+(two*(es-egg)*gamm)**2) if(Ms.eq.MP+one)then c CPR M* = M'+1 c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) * ff c [( M'+1 -1 -M ) ( M'+1 1 -M )] ( M' 1 -M'-1 ) im=rthrj16(Js,one, J,MP+one,-one, -M) ip=rthrj16(Js,one, J,MP+one, one, -M) op=rthrj16(JP,one,Js,MP , one, -MP-one) pf=sumr*js2*rcg16(Js,MFs,one,z0,Fs,MFs )**2*ff srr=srr+(im-ip)*op*pf c DOC: input right output right dr =dr +( -ip)*op*pf endif if(Ms.eq.MP-one)then c CPL M* = M'-1 c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) * ff c [( M'-1 -1 -M ) ( M'-1 1 -M )] ( M' -1 -M'+1 ) im=rthrj16(Js,one,J,MP-one,-one, -M) ip=rthrj16(Js,one,J,MP-one, one, -M) om=rthrj16(JP,one,Js,MP,-one, -MP+one) pf=sumr*js2*rcg16(Js,MFs,one,z0,Fs,MFs )**2*ff slr=slr+(im-ip)*om*pf c DOC: input right output left dl =dl +( -ip)*om*pf endif endif if(MFs.lt.M+two)goto 3 endif if(Fs.lt.Js+one)goto 2 endif if(Js.lt.J+one)goto 1 return end c =========================================================== subroutine ux3(srr,slr,J,M,JP,MP,app, 1wexc,wii,gjs,mn,mb,bz,egg,evib,BB,DD,gamm) c uxl/r= c c sum(J*,K*) (-1)^(K*+M+M') (2J*+1) [( J* 1 J) - (J* 1 J ) ] c (M'+/-1 -1 -M) (M'+/-1 1 -M ) c c (J' 1 J*) sum(q q') (J* 1 J ) ( J' 1 J*) al(q,q') c (M' +/-1 -M'-/+1) (K* q 0 ) ( 0 qp -K*) c c (wii^2-wexc^2) * wii c * ------------------------------------ c (wii^2 - wexc^2)^2+4*wii^2*Gamma^2 c c if ise=1 and isame=1, sum only over allowed Js=J=J implicit none real*8 J,M,JP,MP,Js,js2,sumr,z0,one,two,wexc,wii,gamm, 1app,slr,srr,fr,fl,rthrj16,r1m,r1p,r2m,r2p,wbb,s2,Fs,fs2, 2es,gs,ff,ggg,gel,gjs,mn,mb,bz,egg,evib,BB,DD,ess,gf,MFs, 3aispin,Ms,rcg16 integer*4 ispin data z0,one,two/0.0d0,1.0d0,2.0d0/ c g-factor for electron gel=-2.00231930436153d0 s2=one*(one+1) srr=z0 slr=z0 Js=J-two c J* = J-1,J,J+1 1 Js=Js+one if(Js.ge.z0.and.dabs(Js-JP).le.one)then js2=Js*(Js+one) r1m=rthrj16(Js,one, J,-one, one, z0) r1p=rthrj16(Js,one, J, one,-one, z0) r2m=rthrj16(JP,one,Js, z0,-one, one) r2p=rthrj16(JP,one,Js, z0, one,-one) c [(J* 1 J)*(JP 1 J*) + (J* 1 J)*(JP 1 J*)] * (alpha_xx+alpha_yy)/2 c [(-1 1 0) ( 0 -1 1 ) ( 1 -1 0) ( 0 1 -1 )] sumr=(r1m*r2m+r1p*r2p)*app c approximate g-factor of the excited(J*) state: gs=(gjs*mn*(js2-one)-mb)/js2 c approximate energy of the excited(*) state: ess=wii+evib+(js2-one)*(BB-DD*js2) c construct excited states F* = | J* - S | ... | J* + S|, S=1 Fs=Js-two 2 Fs=Fs+one if(Fs.ge.z0.and.dabs(Js-JP).le.one)then fs2=Fs*(Fs+one) c approximate g-factor of the excited(F*) state: gf=ggg(gs,gel*mb,fs2,js2,s2)*bz c MFS = M-2 ... M+2 MFs=M-two-one 3 MFs=MFs+1 if(abs(MFs).le.Fs.and.abs(MFs-MP).le.two)then es=ess-MFs*gf c ff=(es-egg)/(es-egg+wexc)/(es-egg-wexc) wbb=(es-egg+wexc)*(es-egg-wexc) ff=wbb*(es-egg)/(wbb**2+(two*(es-egg)*gamm)**2) do 4 ispin=-1,1 aispin=dble(ispin) Ms=MFs-aispin if(Ms.eq.MP+one)then c CPR M* = M'+1 c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) * ff c [( M'+1 -1 -M ) ( M'+1 1 -M )] ( M' 1 -M'-1 ) fr=(rthrj16(Js,one,J,MP+one,-one, -M) 1 - rthrj16(Js,one,J,MP+one, one, -M)) 1 * rthrj16(JP,one,Js,MP, one, -MP-one)*ff srr=srr+fr*sumr*js2*rcg16(Js,MFs-aispin,one,aispin,Fs,MFs )**2 endif if(Ms.eq.MP-one)then c CPL M* = M'-1 c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) * ff c [( M'-1 -1 -M ) ( M'-1 1 -M )] ( M' -1 -M'+1 ) fl=(rthrj16(Js,one,J,MP-one,-one, -M) 1 - rthrj16(Js,one,J,MP-one, one, -M)) 1 * rthrj16(JP,one,Js,MP,-one, -MP+one)*ff slr=slr+fl*sumr*js2*rcg16(Js,MFs-aispin,one,aispin,Fs,MFs )**2 endif 4 continue endif if(MFs.lt.M+two)goto 3 endif if(Fs.lt.Js+one)goto 2 endif if(Js.lt.J+one)goto 1 return end c ============================================================= subroutine ux2(srr,slr,J,M,JP,MP,app, 1wexc,wii,gjs,mn,mb,bz,egg,evib,BB,DD,gamm) c uxl/r= c c sum(J*,K*) (-1)^(K*+M+M') (2J*+1) [( J* 1 J) - (J* 1 J ) ] c (M'+/-1 -1 -M) (M'+/-1 1 -M ) c c (J' 1 J*) sum(q q') (J* 1 J ) ( J' 1 J*) al(q,q') c (M' +/-1 -M'-/+1) (K* q 0 ) ( 0 qp -K*) c c (wii^2-wexc^2) * wii c * ------------------------------------ c (wii^2 - wexc^2)^2+4*wii^2*Gamma^2 c c if ise=1 and isame=1, sum only over allowed Js=J=J implicit none real*8 J,M,JP,MP,Js,js2,sumr,z0,one,two,wexc,wii,gamm, 1app,slr,srr,fr,fl,rthrj16,r1m,r1p,r2m,r2p,wbb, 2es,gs,ff,gjs,mn,mb,bz,egg,evib,BB,DD,ess data z0,one,two/0.0d0,1.0d0,2.0d0/ c g-factor for electron srr=z0 slr=z0 Js=J-two c J* = J-1,J,J+1 1 Js=Js+one if(Js.ge.z0.and.dabs(Js-JP).le.one)then js2=Js*(Js+one) c approximate energy of the excited(*) state: ess=wii+evib+(js2-one)*(BB-DD*js2) c approximate g-factor of the excited(*) state: gs=(gjs*mn*(js2-one)-mb)/js2*bz r1m=rthrj16(Js,one, J,-one, one, z0) r1p=rthrj16(Js,one, J, one,-one, z0) r2m=rthrj16(JP,one,Js, z0,-one, one) r2p=rthrj16(JP,one,Js, z0, one,-one) sumr=(r1m*r2m+r1p*r2p)*app c [(J* 1 J)*(JP 1 J*) + (J* 1 J)*(JP 1 J*)] * (alpha_xx+alpha_yy)/2 c [(-1 1 0) ( 0 -1 1 ) ( 1 -1 0) ( 0 1 -1 )] c CPR M* = M'+1 c write(6,*)' J J* ess wii,gs,egg,wexc',J,Js,ess,wii,gs,egg,wexc es=ess-(MP+one)*gs c ff=(es-egg)/(es-egg+wexc)/(es-egg-wexc) wbb=(es-egg+wexc)*(es-egg-wexc) ff=wbb*(es-egg)/(wbb**2+(two*(es-egg)*gamm)**2) c write(6,*)' fp',ff c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) * ff c [( M'+1 -1 -M ) ( M'+1 1 -M )] ( M' 1 -M'-1 ) fr=(rthrj16(Js,one,J,MP+one,-one, -M) 1 - rthrj16(Js,one,J,MP+one, one, -M)) 1 * rthrj16(JP,one,Js,MP, one, -MP-one)*ff c CPL M* = M'-1 es=ess-(MP-one)*gs wbb=(es-egg+wexc)*(es-egg-wexc) ff=wbb*(es-egg)/(wbb**2+(two*(es-egg)*gamm)**2) c write(6,*)' fm',ff c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) * ff c [( M'-1 -1 -M ) ( M'-1 1 -M )] ( M' -1 -M'+1 ) fl=(rthrj16(Js,one,J,MP-one,-one, -M) 1 - rthrj16(Js,one,J,MP-one, one, -M)) 1 * rthrj16(JP,one,Js,MP,-one, -MP+one)*ff slr=slr+fl*sumr*js2 srr=srr+fr*sumr*js2 endif if(Js.lt.J+one)goto 1 return end c ============================================================ subroutine ux0(srr,slr,J,M,JP,MP,azz,app,Jeven,isame, 1ise) c uxl/r= c c sum(J*,K*) (-1)^(K*+M+M') (2J*+1) [( J* 1 J) - (J* 1 J ) ] c (M'+/-1 -1 -M) (M'+/-1 1 -M ) c c (J' 1 J*) sum(q q') (J* 1 J ) ( J' 1 J*) al(q,q') c (M' +/-1 -M'-/+1) (K* q 0 ) ( 0 qp -K*) c c if ise=1 and isame=1, sum only over allowed Js=J=J implicit none real*8 J,M,JP,MP,Js,js2,sumr,z0,one,two,azz, 1app,slr,srr,fr,fl,rthrj16,r1m,r10,r1p,r2m,r20,r2p integer*4 isame,ise logical allowed,Jseven,Jeven data z0,one,two/0.0d0,1.0d0,2.0d0/ srr=z0 slr=z0 Js=J-two if(Js.lt.-one)Js=Js+one c J* = J-1,J,J+1 1 Js=Js+one if(dabs(Js-JP).le.one)then if(ise.eq.0)then allowed=.true. else if(isame.eq.1)then allowed=.true. else Jseven=mod(nint(Js),2).eq.0 allowed=Jseven.eqv.Jeven endif endif if(allowed)then js2=Js*(Js+one) r1m=-rthrj16(Js,one, J,-one, one, z0) r10= rthrj16(Js,one, J, z0, z0, z0) r1p=-rthrj16(Js,one, J, one,-one, z0) r2m=-rthrj16(JP,one,Js, z0,-one, one) r20= rthrj16(JP,one,Js, z0, z0, z0) r2p=-rthrj16(JP,one,Js, z0, one,-one) sumr=r10*r20*azz+(r1m*r2m+r1p*r2p)*app c ( J* 1 J ) * ( JP 1 J* ) * (-1)^(K*+K+M+MP+1) *alpha(q,qp): c ( K* q 0 ) ( 0 qp -K* ) c CPR c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) c [( M'+1 -1 -M ) ( M'+1 1 -M )] ( M' 1 -M'-1 ) fr=(rthrj16(Js,one,J,MP+one,-one, -M) 1 - rthrj16(Js,one,J,MP+one, one, -M)) 1 * rthrj16(JP,one,Js,MP, one, -MP-one) c c CPL c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) c [( M'-1 -1 -M ) ( M'-1 1 -M )] ( M' -1 -M'+1 ) fl=(rthrj16(Js,one,J,MP-one,-one, -M) 1 - rthrj16(Js,one,J,MP-one, one, -M)) 1 * rthrj16(JP,one,Js,MP,-one, -MP+one) slr=slr+fl*sumr*js2 srr=srr+fr*sumr*js2 endif endif if(Js.lt.J+one)goto 1 return end c ============================================================ subroutine ux1(srr,slr,J,M,JP,MP,app,tf,te,qp,bz,gjs,Jsi,JPi) c uxl/r= c c sum(J*,K*) (-1)^(K*+M+M') (2J*+1) [( J* 1 J) - (J* 1 J ) ] c (M'+/-1 -1 -M) (M'+/-1 1 -M ) c c (J' 1 J*) sum(q q') (J* 1 J ) ( J' 1 J*) al(q,q') c (M' +/-1 -M'-/+1) (K* q 0 ) ( 0 qp -K*) c c if ise=1 and isame=1, sum only over allowed Js=J=J implicit none real*8 J,M,JP,MP,Js,js2,sumr,z0,one, 1app,slr,srr,fr,fl,rthrj16,r1m,r1p,r2m,r2p, 1ggg,tf,te,qp,fl2,gfs,bz,gjs,es,p,gel integer*4 L,Jsi,JPi data z0,one/0.0d0,1.0d0/ real*8 na,hprime,cvl,r,mb,Debye,thresh,mn common/rpar/na,hprime,cvl,r,mb,Debye,thresh,mn c g-factor for electron gel=-2.00231930436153d0 c excited state electronic angular momentum = 1 L=1 fl2=dble(L*(L+1)) srr=z0 slr=z0 Js=dble(Jsi) js2=Js*(Js+one) if(abs(Jsi-JPi).le.1)then r1m=-rthrj16(Js,one, J,-one, one, z0) r1p=-rthrj16(Js,one, J, one,-one, z0) r2m=-rthrj16(JP,one,Js, z0,-one, one) r2p=-rthrj16(JP,one,Js, z0, one,-one) sumr=(r1m*r2m+r1p*r2p)*app c ( J* 1 J ) * ( JP 1 J* ) * (-1)^(K*+M+MP+1) *alpha(q,qp): c ( K* q 0 ) ( 0 qp -K* ) c because K* = omega = Lz = +/- 1, alpha(0,0) is not there c gfs=ggg(gjs*mn,gel*mb,js2,js2,fl2)*bz es=-(MP+one)*gfs p=dsqrt(exp(-es/tf/Te)/qp) c CPR c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) c [( M'+1 -1 -M ) ( M'+1 1 -M )] ( M' 1 -M'-1 ) fr=(rthrj16(Js,one,J,MP+one,-one, -M) 1 - rthrj16(Js,one,J,MP+one, one, -M)) 1 * rthrj16(JP,one,Js,MP, one, -MP-one) srr=srr+fr*sumr*js2*p c state |Fs K Ms> es=-(MP-one)*gfs p=dsqrt(exp(-es/tf/Te)/qp) c c CPL c input: output: c [( J* 1 J ) -( J* 1 J )] * ( J' 1 J* ) c [( M'-1 -1 -M ) ( M'-1 1 -M )] ( M' -1 -M'+1 ) fl=(rthrj16(Js,one,J,MP-one,-one, -M) 1 - rthrj16(Js,one,J,MP-one, one, -M)) 1 * rthrj16(JP,one,Js,MP,-one, -MP+one) slr=slr+fl*sumr*js2*p endif return end c ============================================================ subroutine rrspec(Jmin,Jmax,lsym,AA,BB,CC,AE,BE,CE,nm,air,aii, 1wmin,wmax,npo,tw) implicit none integer*4 Jmin,Jmax,J,Ndim,IERR,JE,Ndime,id,ide,nm,im,ix,jx, 1jd,jde,jdeend,jdestart,jdend,jdstart,npo,K,KE real*8 AA,BB,CC,AE,BE,CE,Te,qp,wo,aj,aje,mf,MFE,f,MJ,nt, 1gf,gfe,S2,J2,JE2,F2,FE2,aii(nm,3,3),air(nm,3,3),eee, 1ascr(3,3),asci(3,3),ak,ake,tf,tw(1,nm),c1,co,egg,ileft,iright, 1rcg,MJE,cg1,cg2,rami,roai,FE,p,delw,wmax,wmin,ilr, 1ili,irr,iri,rr,ri,lr,li real*8,allocatable::H(:,:),Cjk(:,:),et(:),ek(:) real*8,allocatable::He(:,:),Cjke(:,:),ete(:),eke(:) real*8,allocatable::ram(:) ,roa(:) logical lsym real*8 z0,half,one,two,three data z0,half,one,two,three 1/0.0d0,0.5d0,1.0d0,2.0d0,3.0d0/ real*8 d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth real*8 na,hprime,cvl,r,mb,Debye,thresh,mn common/rpar/na,hprime,cvl,r,mb,Debye,thresh,mn IERR=0 allocate(ram(npo),roa(npo)) write(6,*)nm,' transition frequencies (cm):' do 10 im=1,nm 10 write(6,609)im,tw(1,im) 609 format(i4,f10.2) do 9 id=1,npo ram(id)=z0 9 roa(id)=z0 Te=t c kelvin to cm-1 tf=0.695d0 delw=(wmax-wmin)/dble(npo-1) nt=dble(npo)/(wmax-wmin) c approximate partition function: qp=dsqrt(pi*(three*Te*tf)**3/(AA-BB+CC)**3) S2=half*(half+1) do 100 J=Jmin,Jmax Ndim=2*J+1 aj=dble(J) J2=aj*(aj+1) allocate(H(Ndim,Ndim),Cjk(Ndim,Ndim),et(Ndim),ek(Ndim)) call tz(Cjk,Ndim) call tz(H,Ndim) if(lsym)then call ddiag(J,ek,Cjk,AA,BB,CC,Ndim) else call seth(h,J,AA,BB,CC) call TRED12(Ndim,H,Cjk,ek,2,IERR,et) endif c rotational, not total quantum number: do 1001 JE=max(0,J-1),J+1 write(6,*)J,JE aje=dble(JE) Ndime=2*JE+1 JE2=aje*(aje+1) allocate(He(Ndime,Ndime),Cjke(Ndime,Ndime),ete(Ndime),eke(Ndime)) call tz(Cjke,Ndime) call tz(He,Ndime) if(lsym)then call ddiag(JE,eke,Cjke,AE,BE,CE,Ndime) else call seth(He,JE,AE,BE,CE) call TRED12(Ndime,He,Cjke,eke,2,IERR,ete) endif c ground state id = Cjk(jd,id)*|JKM>(jd): do 103 id=1,Ndim c excited state ide = Cjke(jde,ide)*|JEKEME>(jde): do 103 ide=1,Ndime c transition energy |G> -> |E>: wo=eke(ide)-ek(id) c transition moments, c = Cjke_jde Cjk_jd , etc if(wo.gt.0.0d0)then if(lsym)then jdstart=id jdend=id jdestart=ide jdeend=ide else jdstart=1 jdend=Ndim jdestart=1 jdeend=Ndime endif c total moment-ground F=J-S,J,J+S F=aj-half-half 1 F=F+half if(F.gt.z0)then F2=F*(F+1) c ground state g-factor gf=(F2+S2-J2)/F2*bz*mb mf=-f-one c MF=-F...F 2 mf=mf+one c ground state energy: egg=ek(id)+gf*mf p=exp(-egg/tf/Te)/qp c total moment-excited FE=JE+S c FE=JE-S,JE,JE+S fe=aje-half-half 3 fe=fe+half FE2=FE*(FE+1) c excited state g-factor gfe=(FE2+S2-JE2)/FE2*bz*mb c MFE=-FE...FE,|MF-MFE|<=2 MFE=-fe-one 5 MFE=MFE+one if(dabs(MFE-mf).le.two)then c excited state energy for each vibration: do 8 im=1,nm eee=tw(1,im)+ek(ide)+gfe*MFE-egg do 82 ix=1,3 do 82 jx=1,3 asci(ix,jx)=aii(im,ix,jx) 82 ascr(ix,jx)=air(im,ix,jx) ilr=0.0d0 ili=0.0d0 irr=0.0d0 iri=0.0d0 do 81 jd=jdstart,jdend K=-J+jd-1 ak=dble(K) c1=Cjk(jd,id) do 81 jde=jdestart,jdeend KE=-JE+jde-1 if(abs(K-KE).le.2)then ake=dble(KE) co=Cjke(jde,ide)*c1 c MJ=MF-SZ,MF+SZ MJ=MF-half-half 6 MJ=MJ+half c Clebsh-Gordon cg1=rcg(aj,mj,half,mf-mj,f,mf)*co c MJE=MFE-SZ,MF+SZ c MJE=MFE-MF+MJ becasue of nuclear spin MJE=MFE-MF+MJ c Clebsh-Gordon cg2=cg1*rcg(aje,mje,half,MFE-mje,fe,MFE) if(cg2.gt.z0)then call ux(ri,rr,lr,li,aj,ak,MJ,aje,ake,MJE,ascr,asci) ilr =ilr +cg2*lr ili =ili +cg2*li irr=irr+cg2*rr iri =iri +cg2*ri endif if(MJ.lt.MF+half)goto 6 endif 81 continue ileft =J2*JE*(ilr**2 +ili**2 ) iright=J2*JE*(irr**2+iri**2) rami=(iright+ileft)*p roai=(iright-ileft)*p if(rami.gt.1.0d-9)then call spread1(eee,rami,ram,npo,wmin,nt,delw) call spread1(eee,roai,roa,npo,wmin,nt,delw) endif 8 continue endif c MFE=MF +/-2 if(MFE.lt.MFE)goto 5 if(FE.lt.aje+half)goto 3 if(MF.lt.F)goto 2 endif c F>0 if(F.lt.aj+half)goto 1 endif c wo>0 103 continue c id,ide 1001 deallocate(He,Cjke,ete,eke) c JE 100 deallocate(H,Cjk,et,ek) c J call wrs('ramJ.prn',ram,wmin,wmax,npo) call wrs('roaJ.prn',roa,wmin,wmax,npo) return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(N,N),Z(N,N),D(N),E(N) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 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 DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) JP1 = J + 1 IF (L .LT. JP1) GO TO 220 DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE HH = F / (H + H) DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) 290 D(I) = H 300 CONTINUE 320 D(1) = 0.0D0 E(1) = 0.0D0 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 DO 360 J = 1, L G = 0.0D0 DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE 500 CONTINUE 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Z(N,N),D(N),E(N) 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 IERR = 0 IF (N .EQ. 1) GO TO 1001 DO 100 I = 2, N 100 E(I-1) = E(I) F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 110 CONTINUE 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 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) DO 140 I = L1, N 140 D(I) = D(I) - H F = F + H P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L 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)) 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 200 CONTINUE 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 GO TO 1001 1000 IERR = L 1001 RETURN END c ================================================================== subroutine trt(g,a) implicit none real*8 g(3,3),a(3,3),t(3,3) integer*4 i,j,ii,jj do 1 i=1,3 do 1 j=1,3 t(i,j)=0.0d0 do 1 ii=1,3 1 t(i,j)=t(i,j)+g(ii,j)*a(ii,i) do 2 i=1,3 do 2 j=1,3 g(i,j)=0.0d0 do 2 jj=1,3 2 g(i,j)=g(i,j)+t(i,jj)*a(jj,j) return end c ================================================================== subroutine trv(g,a) implicit none real*8 g(3),a(3,3),t(3) integer*4 i,ii do 1 i=1,3 t(i)=0.0d0 do 1 ii=1,3 1 t(i)=t(i)+g(ii)*a(ii,i) do 2 i=1,3 2 g(i)=t(i) return end c ================================================================== subroutine tz(C,N) implicit none real*8 C(N,N) integer*4 i,j,N do 1 i=1,N do 1 j=1,N 1 C(i,j)=0.0d0 return end c ================================================================== subroutine ddiag(J,ek,Cjk,AA,BB,CC,n) implicit none integer*4 J,id,K,n real*8 ek(n),Cjk(n,n),AA,BB,CC,JJP,K2 id=0 JJP=dble(J*(J+1)) do 1012 K=-J,J id=id+1 K2=dble(K*K) ek(id)=0.5d0*(AA*(JJP-K2)/2.0d0+BB*(JJP-K2)/2.0d0+CC*K2) 1012 Cjk(id,id)=1.0d0 return end c ================================================================== subroutine sethno2(h,J,A,B,C,aa,bb,cc,N) c add non-zero elements of Hamiltonian c NO2 molecule aa bb cc -- spin-rotational coupling constants c J rotational number, N total momentum (N = S + J) implicit none integer*4 J,i,K real*8 h(2*J+1,2*J+1),K2,KKP,KKQ,A,B,C,JJP,aa,bb,cc,en,N, 1NNP NNP=N*(N+1) JJP=dble(J*(J+1)) i=0 do 101 K=-J,J i=i+1 K2=dble(K*K) KKP=dble(K*(K+1)) KKQ=dble((K+1)*(K+2)) en=0.50d0*(bb+cc)+(aa-(bb+cc)/2.0d0)*K2/JJP c +/- delta(|K|,1)*(bb+cc) neglected c : H(i,i)=0.5d0*(B+C)*(JJP-K2)+A*K2 1+0.5d0*en*(NNP-JJP-0.75d0) c : c : if(K+2.le.J)then H(i,i+2)=0.25d0*(B-C)*dsqrt((JJP-KKP)*(JJP-KKQ)) H(i+2,i)=H(i,i+2) endif 101 continue return end c ================================================================== subroutine seth(h,J,A,B,C) c add non-zero elements of Hamiltonian implicit none integer*4 J,i,K real*8 h(2*J+1,2*J+1),K2,KKP,KKQ,A,B,C,JJP JJP=dble(J*(J+1)) i=0 do 101 K=-J,J i=i+1 K2=dble(K*K) KKP=dble(K*(K+1)) KKQ=dble((K+1)*(K+2)) c : H(i,i)=0.5d0*(A+B)*(JJP-K2)+C*K2 c : c : if(K+2.le.J)then H(i,i+2)=0.25d0*(A-B)*dsqrt((JJP-KKP)*(JJP-KKQ)) H(i+2,i)=H(i,i+2) endif 101 continue return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine spread(ic,e,a,d,r,npo,wmin,nt,delw) c c ic .. left/right CPL, sign of r implicit none real*8 e,a,d(*),x,delw,wmin,nt,ass,ff,ic,r(*) integer*4 npo,i,n real*8 d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth c approximate position of e: c n=nint((e-wmin)/(wmax-wmin)*dble(npo)) n=nint((e-wmin)*nt) x=wmin+dble(n-2)*delw do 1 i=n,npo x=x+delw ass=-((x-e)/bandwidth)**2 if(ass.gt.-25.0d0)then ff=exp(ass)/bandwidth/dsqrtpi d(i)=d(i)+a*ff/2.0d0 r(i)=r(i)+ic*a*ff else goto 2 endif 1 continue 2 x=wmin+dble(n-1)*delw do 3 i=n-1,1,-1 x=x-delw ass=-((x-e)/bandwidth)**2 if(ass.gt.-25.0d0)then ff=exp(ass)/bandwidth/dsqrtpi d(i)=d(i)+a*ff/2.0d0 r(i)=r(i)+ic*a*ff else goto 4 endif 3 continue 4 continue return end subroutine spread1(e,a,d,npo,wmin,nt,delw) c implicit none real*8 e,a,d(*),x,delw,wmin,nt,ass,ff integer*4 npo,i,n real*8 d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth c approximate position of e: c n=nint((e-wmin)/(wmax-wmin)*dble(npo)) n=nint((e-wmin)*nt) if(n.lt.1.or.n.gt.npo)return x=wmin+dble(n-2)*delw do 1 i=n,npo x=x+delw ass=-((x-e)/bandwidth)**2 if(ass.gt.-25.0d0)then ff=exp(ass)/bandwidth/dsqrtpi d(i)=d(i)+a*ff else goto 2 endif 1 continue 2 x=wmin+dble(n-1)*delw do 3 i=n-1,1,-1 x=x-delw ass=-((x-e)/bandwidth)**2 if(ass.gt.-25.0d0)then ff=exp(ass)/bandwidth/dsqrtpi d(i)=d(i)+a*ff else goto 4 endif 3 continue 4 continue return end subroutine wrs(f,s,wmin,wmax,npo) implicit none character*(*) f real*8 s(*),x,wmin,wmax,delw integer*4 npo,i open(9,file=f) delw=(wmax-wmin)/dble(npo-1) x=wmin-delw do 1 i=1,npo x=x+delw 1 write(9,900)x,s(i) 900 format(f15.4,g20.8) close(9) write(6,*)f,' written' return end SUBROUTINE READTTTQ(nm,n,alphad,Gpd,ad,ft,tw,ip,ntw,ntmax) IMPLICIT none INTEGER*4 nm,nq,LL,I,J,K,n,ip,ntmax,ntw(*) real*8 alphad(n,nm,3,3),Gpd(n,nm,3,3),Ad(n,nm,3,3,3),tw(n,ntmax) character*(*) ft open(15,file=ft) READ(15,*) READ(15,*)nq if(nq.ne.ntw(ip))call report('nq<>ntw') READ(15,*) READ(15,*) DO 1 LL=1,nq READ(15,*)tw(ip,LL),tw(ip,LL) DO 1 I=1,3 1 READ(15,*)alphad(ip,LL,I,1),(alphad(ip,LL,I,J),J=1,3) tw(ip,nq+1)=0.0d0 READ(15,*) READ(15,*) DO 2 LL=1,nq READ(15,*) DO 2 I=1,3 2 READ(15,*)GPD(ip,LL,I,1),(GPD(ip,LL,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 LL=1,nq READ(15,*) DO 3 I=1,3 DO 3 J=1,3 3 READ(15,*)(Ad(ip,LL,I,J,K),K=1,2),(Ad(ip,LL,I,J,K),K=1,3) c first index electric CLOSE(15) RETURN END subroutine setc(cr,ci) implicit none real*8 cr(3,5,3,3),ci(3,5,3,3) integer*4 ij,ik,ix,jx do 2 ij=1,3 c J = 0 1 2 do 2 ik=1,5 c K = -J ... +J do 2 ix=1,3 do 2 jx=1,3 ci(ij,ik,ix,jx)=0.0d0 2 cr(ij,ik,ix,jx)=0.0d0 c J=0 K=0 cr(1,1,1,1)=-1.0d0/dsqrt(3.0d0) cr(1,1,2,2)=-1.0d0/dsqrt(3.0d0) cr(1,1,3,3)=-1.0d0/dsqrt(3.0d0) c J=1 K=-1 cr(2,1,1,3)=-0.5d0 cr(2,1,3,1)= 0.5d0 ci(2,1,2,3)= 0.5d0 ci(2,1,3,2)=-0.5d0 c J=1 K=0 ci(2,2,1,2)= 1.0d0/dsqrt(2.0d0) ci(2,2,2,1)=-1.0d0/dsqrt(2.0d0) c J=1 K=1 cr(2,3,1,3)=-0.5d0 cr(2,3,3,1)= 0.5d0 ci(2,3,2,3)=-0.5d0 ci(2,3,3,2)= 0.5d0 c J=2 K=-2 cr(3,1,1,1)= 0.5d0 cr(3,1,2,2)=-0.5d0 ci(3,1,1,2)=-0.5d0 ci(3,1,2,1)=-0.5d0 c J=2 K=-1 cr(3,2,1,3)= 0.5d0 cr(3,2,3,1)= 0.5d0 ci(3,2,3,2)=-0.5d0 ci(3,2,2,3)=-0.5d0 c J=2 K=0 cr(3,3,1,1)= -1.0d0/dsqrt(6.0d0) cr(3,3,2,2)= -1.0d0/dsqrt(6.0d0) ci(3,3,3,3)= 2.0d0/dsqrt(6.0d0) c J=2 K=1 cr(3,4,1,3)=-0.5d0 cr(3,4,3,1)=-0.5d0 ci(3,4,2,3)=-0.5d0 ci(3,4,3,2)=-0.5d0 c J=2 K=2 cr(3,5,1,1)= 0.5d0 cr(3,5,2,2)=-0.5d0 ci(3,5,1,2)= 0.5d0 ci(3,5,2,1)= 0.5d0 return end Function Fac(n) implicit none integer*4 n,ii real*8 Fac,bb bb=1.0d0 do 1 ii=2 , n 1 bb=bb*dble(ii) fac=bb return end Function Fac16(n) implicit none integer*4 n,ii real*16 Fac16,bb bb=1.0q0 do 1 ii=2 , n 1 bb=bb*real(ii,16) fac16=bb return end subroutine Initfac implicit none integer*4 nfac0,i parameter (nfac0=170) real*8 f0,factorial,fac real*16 f016,fc16,fac16 common/npar/f0,factorial(nfac0) common/npar16/f016,fc16(2*nfac0) do 1 i=0 , nfac0 1 factorial(i)=fac(i) do 2 i=0 ,2*nfac0 2 fc16(i)=fac16(i) return end Function rcg(j1,m1,j2,m2,j3,m3) c {Clebsh-Gordon coefficients} implicit none real*8 rcg,j1,m1,j2,m2,j3,m3,rthrj if (mod(NINT(j1-j2+m3) , 2).eq. 0) then rcg= dsqrt(2.0d0*j3+1.0d0)*rthrj(j1,j2,j3,m1,m2,-m3) else rcg=-dsqrt(2.0d0*j3+1.0d0)*rthrj(j1,j2,j3,m1,m2,-m3) endif return end Function rcg16(j1,m1,j2,m2,j3,m3) c {Clebsh-Gordon coefficients} implicit none real*8 rcg16,j1,m1,j2,m2,j3,m3,rthrj16 if (mod(NINT(j1-j2+m3) , 2).eq. 0) then rcg16= dsqrt(2.0d0*j3+1.0d0)*rthrj16(j1,j2,j3,m1,m2,-m3) else rcg16=-dsqrt(2.0d0*j3+1.0d0)*rthrj16(j1,j2,j3,m1,m2,-m3) endif return end Function gfactor(F,J,lambda,k) c {g-factor according , H.E.Radford, Phys.Rev. (1961)122,1,114} c {or W.Herrmann, W.Rohrbeck and W.Urban, Appl.Phys.22(1980),71-75} c {-gF expression derived by Petr Bour} implicit none real*8 gfactor,F,J,lambda,k,sqy,x,sign,gj sign=3.0d0-2.0d0*k sqy=(j+1.5d0)*(j-0.5d0) x=dsqrt(4.0d0*(j+0.5d0)**2+lambda*(lambda-4.0d0)) gj=(1.5d0+sign*(2.0d0*sqy-1.5d0*lambda+3.0d0)/x)/j/(j+1.0d0) gfactor=gJ*(F*(F+1.0d0)+J*(J+1.0d0)-2.0d0)/2.0d0/F/(F+1.0d0) return end Function Sig(a) implicit none real*8 Sig,a if(mod(NINT(a) , 2).eq. 0) then sig=1.0d0 else sig=-1.0d0 endif return end subroutine cartosphere(ac,asr,asi,n) c second order tensor, Cartesian to spherical implicit none integer*4 n,i,ij,ik,J,ix,jx real*8 ac(n,3,3),asr(n,3,5),cr(3,5,3,3),ci(3,5,3,3),asi(n,3,5) call setc(cr,ci) do 1 i=1,n do 1 ij=1,3 J=ij-1 do 1 ik=1,2*J+1 c K=-J+ik-1 asr(i,ij,ik)=0.0d0 asi(i,ij,ik)=0.0d0 do 1 ix=1,3 do 1 jx=1,3 asr(i,ij,ik)=asr(i,ij,ik)+cr(ij,ik,ix,jx)*ac(i,ix,jx) 1 asi(i,ij,ik)=asi(i,ij,ik)+ci(ij,ik,ix,jx)*ac(i,ix,jx) return end subroutine cartospherei(ac,asr,asi,n) c Cartesian to spherical, individual indices implicit none integer*4 n,i real*8 ac(n,3,3),asr(n,3,3),asi(n,3,3), 1xx,xy,xz,yx,yy,yz,zx,zy,zz c spherical -1 0 1 = 1 2 3 do 1 i=1,n xx=ac(i,1,1) xy=ac(i,1,2) xz=ac(i,1,3) yx=ac(i,2,1) yy=ac(i,2,2) yz=ac(i,2,3) zx=ac(i,3,1) zy=ac(i,3,2) zz=ac(i,3,3) asr(i,1,1)= (xx-yy)/2.0d0 asi(i,1,1)= -(yx+xy)/2.0d0 asr(i,1,3)=-(xx+yy)/2.0d0 asi(i,1,3)= +(yx-xy)/2.0d0 asr(i,3,1)=-(xx+yy)/2.0d0 asi(i,3,1)= -(yx-xy)/2.0d0 asr(i,3,3)= (xx-yy)/2.0d0 asi(i,3,3)= +(yx+xy)/2.0d0 asr(i,1,2)= xz /dsqrt(2.0d0) asi(i,1,2)= -yz /dsqrt(2.0d0) asr(i,2,1)= zx /dsqrt(2.0d0) asi(i,2,1)= -zy /dsqrt(2.0d0) asr(i,3,2)= -xz /dsqrt(2.0d0) asi(i,3,2)= -yz /dsqrt(2.0d0) asr(i,2,3)= -zx /dsqrt(2.0d0) asi(i,2,3)= -zy /dsqrt(2.0d0) asr(i,2,2)= zz 1 asi(i,2,2)= 0.0d0 return end c ============================================================ subroutine lookatout(lzmat,gt,abc,ta,AA,BB,CC,AE,BE,CE,ud) implicit none integer*4 n0,MENDELEV parameter (n0=3000,MENDELEV=89) integer*4 n3,nat,l,i,iz(n0),KA,ig98,j,ia,k,IERR,NG character*80 FN,s80 logical lex,lzmat real*8 bohr,r(3,n0),ama(n0),xp(3,3),qt(3,3),qtot(3,3),debye, 1abc(3),c,g(3,3),cm(3),mt,t(3,3),q(3,3),a(3,3),z(3),u,mp, 1gn(3,3),ge(3,3),gq(3,3),amu,e(3),ta,AA,BB,CC,AE,BE,CE,ud(3), 1qtr(3,3),gt(3,3) CHARACTER*2 atsy(MENDELEV),s2 data atsy/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne', 3'Na','Mg','Al','Si',' P',' S','Cl','Ar', 4' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn', 4 'Ga','Ge','As','Se','Br','Kr', 5'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd', 5 'In','Sn','Sb','Te',' I','Xe', 6'Cs','Ba','La', 6 'Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho', 6 'Er','Tm','Yb','Lu', 6'Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg', 6 'Tl','Pb','Bi','Po','At','Rn', 7'Fr','Ra','Ac'/ real*8 amas(MENDELEV) data amas/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ bohr=0.52917721092d0 mp=1.00727646688d0 c=2.99792458d10 u=219474.63068d0 amu=1822.88839d0 c debye-in au, different definition from the rest of the program: debye=2.541746913d0 do 92 i=1,3 do 92 j=1,3 qt(i,j)=0.0d0 92 qtot(i,j)=0.0d0 inquire(file='G.OUT',exist=lex) if(.not.lex)return OPEN(2,FILE='G.OUT') OPEN(4,FILE='FILE.X') NG=0 1 READ(2,2000,END=1000)FN 2000 FORMAT(A80) c c dalton geometry: IF(FN(3:30).EQ.'Cartesian Coordinates (a.u.)')then NG=NG+1 READ(2,2000,END=1000)FN READ(2,2000,END=1000)FN READ(2,2000,END=1000)FN READ(FN(31:35),*)n3 nat=n3/3 if(nat.gt.n0)call report('too many atoms') do 13 l=1,nat read(2,1313)s2,(r(i,l),i=1,3) do 132 i=1,3 132 r(i,l)=r(i,l)*bohr 1313 format(1x,a2,9x,3(8x,f15.10)) iz(l)=0 do 131 i=1,MENDELEV 131 if(s2.eq.atsy(i))iz(l)=i if(iz(l).lt.1)then write(6,*)'Unknown atom',s2 ama(l)=0.0d0 else ama(l)=amas(iz(l)) endif 13 continue write(4,*)'Dalton output' write(4,*)nat do 34 l=1,nat 34 if(iz(l).gt.0)write(4,4000)iz(l),(r(i,l),i=1,3),(0,i=1,7),0.0d0 endif c Gaussian geometry: IF((lzmat.and.(FN(19:39).EQ.'Z-Matrix orientation:'.OR. 1 FN(26:46).EQ.'Z-Matrix orientation:'.OR. 1 FN(20:37).EQ.'Input orientation:'.OR. 1 FN(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (FN(20:40).EQ.'Standard orientation:'.OR. 2 FN(26:46).EQ.'Standard orientation:')))THEN NG=NG+1 ig98=0 if(FN(26:46).EQ.'Z-Matrix orientation:'.OR. 1 FN(27:44).EQ.'Input orientation:'.OR. 1 FN(26:46).EQ.'Standard orientation:')ig98=1 WRITE(4,2000) DO 4 I=1,4 4 READ(2,*) l=0 5 READ(2,2000)FN IF(FN(2:4).NE.'---')THEN l=l+1 if(l.gt.n0)call report('too many atoms') BACKSPACE 2 if(ig98.eq.0)then READ(2,*)IA,KA,(r(i,l),i=1,3) else READ(2,*)IA,KA,IA,(r(i,l),i=1,3) endif iz(l)=KA IF(KA.EQ.-1)l=l-1 GOTO 5 ENDIF nat=l write(6,*)'Gaussian geometry' write(6,*)nat do 343 l=1,nat 343 if(iz(l).gt.0)write(6,4000)iz(l),(r(i,l),i=1,3),(0,i=1,7),0.0d0 4000 format(I3,3F12.6,7(1x,i1),f4.1) ENDIF c Dalton: IF(FN(30:53).eq.'Dipole moment components')then write(6,2000)FN read(2,*) read(2,*) read(2,*) read(2,*) read(2,206)ud(1) read(2,206)ud(2) read(2,206)ud(3) 206 format(7x,f16.8) endif c Gaussian: IF(FN(2:42).eq.'Dipole moment (field-independent basis, D')then write(6,2000)FN read(2,207)(ud(i),i=1,3) 207 format(3(10x,f16.4)) ud(1)=ud(1)/debye ud(2)=ud(2)/debye ud(3)=ud(3)/debye endif c Dalton: IF(FN(22:61).eq.'Paramagnetic magnetizability tensor (au)')then write(6,2000)FN read(2,*) read(2,*) read(2,*) read(2,*) do 1010 i=1,3 1010 read(2,202)(xp(i,j),j=1,3) 202 format(7x,3F20.12) endif c Gaussian: IF(FN(2:40).eq.'Paramagnetic susceptibility tensor (au)')then write(6,2000)FN do 10 i=1,3 10 read(2,203)(xp(i,j),j=1,3) 203 format(3(6x,f15.4)) endif IF(FN(27:55).eq.'Total quadrupole moments (au)'.or. 1FN(2:43).eq.'Quadrupole moment (field-independent basis')then IF(FN(27:55).eq.'Total quadrupole moments (au)')then c Dalton: write(6,2000)FN read(2,*) read(2,*) read(2,*) read(2,*) do 121 i=1,3 121 read(2,2031)(qtr(i,j),j=1,3) 2031 format(17x,3f12.6) do 1211 i=1,3 do 1211 j=1,3 qt(i,j)=0.0d0 qtot(i,j)=0.0d0 c gaussian traceless quadrupole moment = (2/3) x Dalton quadrupole moment: 1211 qtr(i,j)=qtr(i,j)*2.0d0/3.0d0 else c Gaussian, Debye*Ang: write(6,2000)FN read(2,204)qt(1,1),qt(2,2),qt(3,3) read(2,204)qt(1,2),qt(1,3),qt(2,3) 204 format(3(6x,f20.4)) qt(2,1)=qt(1,2) qt(3,1)=qt(1,3) qt(3,2)=qt(2,3) do 122 i=1,3 do 122 j=1,3 122 qt(i,j)=qt(i,j)/bohr/debye c make inertia-like: do 12 i=1,3 do 12 j=1,3 if(i.eq.j)then k=i+1 if(k.gt.3)k=1 l=k+1 if(l.gt.3)l=1 qtot(i,j)=qt(k,k)+qt(l,l) else qtot(i,j)=-qt(i,j) endif 12 continue c traceless: qtr(1,1)=(2.0d0*qt(1,1)-qt(2,2)-qt(3,3))/3.0d0 qtr(2,2)=(2.0d0*qt(2,2)-qt(3,3)-qt(1,1))/3.0d0 qtr(3,3)=(2.0d0*qt(3,3)-qt(1,1)-qt(3,3))/3.0d0 qtr(1,2)=qt(2,1) qtr(2,1)=qt(1,2) qtr(1,3)=qt(3,1) qtr(3,1)=qt(1,3) qtr(2,3)=qt(3,2) qtr(3,2)=qt(2,3) endif endif c Dalton: IF(FN(2:21).eq.'Rotational constants')then READ(2,2000,END=1000)s80 IF(s80(2:21).eq.'--------------------')then read(2,*) read(2,*) read(2,*) read(2,*)(ABC(I),I=1,3) ABC(1)=ABC(1)/1000.0d0 ABC(2)=ABC(2)/1000.0d0 ABC(3)=ABC(3)/1000.0d0 else backspace 2 endif endif c Gaussian: IF(FN(2:28).eq.'Rotational constants (GHZ):')then write(6,2000)FN read(FN(30:80),*)(ABC(I),I=1,3) endif c Dalton: IF(FN(3:34).eq.'Molecular rotational g-tensor in')then write(6,2000)FN read(2,*) read(2,*) do 81 i=1,3 81 read(2,2011)(g(i,j),j=1,3) 2011 format(10x,3f16.8) endif c Gaussian: IF(FN(2:24).eq.'GIAO rotational tensor:')then write(6,2000)FN do 8 i=1,3 8 read(2,205)(g(i,j),j=1,3) 205 format(3(6x,f11.4)) endif c Gaussian: IF(FN(10:26).eq.'has atomic number')then write(6,2000)FN read(FN(6:8),*)I read(FN(35:44),*)ama(I) endif IF(FN(13:29).eq.'has atomic number')then write(6,2000)FN read(FN(9:11),*)I read(FN(42:51),*)ama(I) endif GOTO 1 1000 CLOSE(2) if(ng.eq.0)then write(*,*)'Geometry not found' else WRITE(*,*)' File FILE.X written' write(*,*)nat,' atoms' endif c c coordinates to atomic units: do 2 ia=1,nat do 2 i=1,3 2 r(i,ia)=r(i,ia)/bohr c subtract mass center: DO 21 I=1,3 CM(I)=0.0d0 mt=0.0d0 DO 31 J=1,nat mt=mt+ama(J) 31 CM(I)=CM(I)+R(I,J)*ama(J) 21 CM(I)=CM(I)/mt do 22 ia=1,nat do 22 i=1,3 22 r(i,ia)=r(i,ia)-CM(i) write(6,*)'mass center (A):' write(6,6001)(cm(i)*bohr,i=1,3) c c inertia and quadrupole moment: do 6 i=1,3 t(i,i)=0.0d0 q(i,i)=0.0d0 j=i+1 if(j.gt.3)j=1 k=j+1 if(k.gt.3)k=1 do 61 ia=1,nat q(i,i)=q(i,i)+dble(iz(ia))*(r(j,ia)*r(j,ia)+r(k,ia)*r(k,ia)) 61 t(i,i)=t(i,i)+ama(ia)*(r(j,ia)*r(j,ia)+r(k,ia)*r(k,ia)) do 6 j=1,i-1 t(i,j)=0.0d0 q(i,j)=0.0d0 do 62 ia=1,nat q(i,j)=q(i,j)-dble(iz(ia))*r(i,ia)*r(j,ia) 62 t(i,j)=t(i,j)- ama(ia)*r(i,ia)*r(j,ia) q(j,i)=q(i,j) 6 t(j,i)=t(i,j) CALL TRED12(3,T,A,Z,2,IERR,E) do 7 i=1,3 7 z(i)=u*c/1.0d9/z(i)/2.0d0 write(6,*)'Gaussian (GHz):' write(6,6000)(abc(i),i=1,3) write(6,*)'Recalc:' write(6,6000)(Z(i)/amu,i=1,3) 6000 format(3f18.6) do 15 i=1,3 do 15 j=1,3 if(t(i,i).ne.0.0d0)ge(i,j)=-4.0d0*xp(i,j)/t(i,i)*mp if(t(i,i).ne.0.0d0)gn(i,j)=q(i,j)/t(i,i)*mp 15 gt(i,j)=ge(i,j)+gn(i,j) do 14 i=1,3 do 14 j=1,3 14 if(t(i,i).ne.0.0d0)gq(i,j)=qtot(i,j)/t(i,i)*mp write(6,*)'Laboratory system:' write(6,*)'xp:' write(6,6001)((xp(i,j),i=1,3),j=1,3) write(6,*)'g-tensor,electronic:' write(6,6001)((ge(i,j),i=1,3),j=1,3) write(6,*)'g-tensor,nuclear:' write(6,6001)((gn(i,j),i=1,3),j=1,3) write(6,*)'sum:' write(6,6001)((gt(i,j),i=1,3),j=1,3) write(6,*)'g-tensor in output:' write(6,6001)((g(i,j),i=1,3),j=1,3) write(6,*)'Quadrupole moment nuclear:' write(6,6001)((q(i,j),i=1,3),j=1,3) write(6,*)'Quadrupole moment total, inertia-like:' write(6,6001)((qtot(i,j),i=1,3),j=1,3) write(6,*)'Quadrupole moment total, traceless:' write(6,6001)((qtr(i,j),i=1,3),j=1,3) write(6,*)'Inertia moment:' write(6,6001)((t(i,j),i=1,3),j=1,3) write(6,*)'g-tensor from quadrupole:' write(6,6001)((gq(i,j),i=1,3),j=1,3) write(6,*)'dipole (au):' write(6,6001)(ud(i),i=1,3) 6001 format(/,3(3f18.6,/),/) write(6,*)'Trafo:' write(6,6001)((a(i,j),i=1,3),j=1,3) call trt(g,a) call trt(t,a) call trt(q,a) call trt(xp,a) call trt(qtot,a) call trt(qtr,a) call trv(ud,a) do 9 i=1,3 write(6,*)t(i,i),mp if(t(i,i).ne.0.0d0)then do 91 j=1,3 ge(i,j)=-4.0d0*xp(i,j)/t(i,i)*mp gn(i,j)=q(i,j)/t(i,i)*mp gt(i,j)=ge(i,j)+gn(i,j) 91 gq(i,j)=qtot(i,j)/t(i,i)*mp endif 9 continue write(6,*) write(6,*)'Molecular System:' write(6,*)'*****************' write(6,*)'xp:' write(6,6001)((xp(i,j),i=1,3),j=1,3) write(6,*)'g-tensor,electronic:' write(6,6001)((ge(i,j),i=1,3),j=1,3) write(6,*)'g-tensor,nuclear:' write(6,6001)((gn(i,j),i=1,3),j=1,3) write(6,*)'sum:' write(6,6001)((gt(i,j),i=1,3),j=1,3) write(6,*)'g-tensor from output:' write(6,6001)((g(i,j),i=1,3),j=1,3) write(6,*)'Quadrupole moment nuclear:' write(6,6001)((q(i,j),i=1,3),j=1,3) write(6,*)'Quadrupole moment total, inertia-like:' write(6,6001)((qtot(i,j),i=1,3),j=1,3) write(6,*)'Quadrupole moment total, traceless:' write(6,6001)((qtr(i,j),i=1,3),j=1,3) write(6,*)'Inertia moment:' write(6,6001)((t(i,j),i=1,3),j=1,3) write(6,*)'g-tensor from quadrupole:' write(6,6001)((gq(i,j),i=1,3),j=1,3) write(6,*)'dipole (au):' write(6,6001)(ud(i),i=1,3) if(AA+BB+CC.eq.0.0d0)then write(6,*)' Rotational moments from G.OUT used' AA=abc(1)*ta BB=abc(2)*ta CC=abc(3)*ta AE=abc(1)*ta BE=abc(2)*ta CE=abc(3)*ta endif return end c ============================================================ subroutine setsign(Ks,sr,si) implicit none c (-1)^Ks: real*8 Ks,sr,si,z0,one,two,four data z0,one,two,four 1/0.0d0,1.0d0,2.0d0,4.0d0/ sr=z0 si=z0 if(mod(nint(Ks),2).eq.0)then sr=one else if(mod(nint(two*Ks),2).eq.0)then sr=-one else if(mod(nint(four*Ks)+2,8).eq.0)then si=-one else si=one endif endif endif return end c ============================================================ subroutine setxa(k,xe,aje,bje,jj,ye) implicit none integer*4 k real*8 four,two,half,xe,aje,bje,jj,ye data half,two,four/0.5d0,2.0d0,4.0d0/ xe=dsqrt(four*(jj+half)**2+ye*(ye-four)) if(k.eq.2) then aje=-dsqrt((xe-(ye-two))/two/xe) bje= dsqrt((xe+(ye-two))/two/xe) else aje= dsqrt((xe+(ye-two))/two/xe) bje= dsqrt((xe-(ye-two))/two/xe) endif return end c ============================================================ subroutine incr(ilr,ili,irr,iri,lr,li,rr,ri,c22) implicit none real*8 ilr,ili,irr,iri,lr,li,rr,ri,c22 ilr =ilr +lr*c22 ili =ili +li*c22 irr=irr+rr*c22 iri =iri +ri*c22 return end Function rthrj(x1,x2,x3,y1,y2,y3) c {real three - J symbols} implicit none integer*4 nfac0 parameter (nfac0=170) real*8 f0,factorial common/npar/f0,factorial(nfac0) real*8 rthrj,x1,x2,x3,y1,y2,y3,ar,u,sum,msign integer*4 j1d,j2d,j3d,m1d,m2d,m3d,s,ni,i1,i2,i3,i4,i5,i6 j1d=NINT(x1*2.0d0) j2d=NINT(x2*2.0d0) j3d=NINT(x3*2.0d0) m1d=NINT(y1*2.0d0) m2d=NINT(y2*2.0d0) m3d=NINT(y3*2.0d0) if(abs(m1d).gt.j1d.or.abs(m2d).gt.j2d.or.abs(m3d).gt.j3d.or. 1m1d+m2d .ne.-m3d.or.j1d+j2d.lt.j3d.or.j2d+j3d.lt.j1d.or. 2j3d+j1d.lt.j2d .or.mod(j1d-m1d,2).ne.0 .or. 3mod(j2d-m2d,2) .ne.0.or.mod(j3d-m3d,2).ne. 0.or. 4mod(j1d+j2d+j3d,2).ne.0) then ar=0.0d0 else if (mod((j1d-j2d-m3d)/2,2).eq.0) then msign=1.0d0 else msign=-1.0d0 endif s=NINT(x1+x2+x3) ar=factorial(s-j3d)*factorial(s-j2d)*factorial(s-j1d) 1 /factorial(s+1) ar=ar*factorial((j1d+m1d) / 2)*factorial((j1d-m1d) / 2) ar=ar*factorial((j2d+m2d) / 2)*factorial((j2d-m2d) / 2) ar=ar*factorial((j3d+m3d) / 2)*factorial((j3d-m3d) / 2) ar=msign*dsqrt(ar) sum=0.0d0 do 1 ni=0 , s i1=ni i2=(j1d+j2d-j3d) / 2-ni i3=(j1d-m1d) / 2-ni i4=(j2d+m2d) / 2-ni i5=(j3d-j2d+m1d) / 2+ni i6=(j3d-j1d-m2d) / 2+ni if( mod(ni, 2) .ne. 0) then msign=-1.0d0 else msign=1.0d0 endif if(((((i2.ge.0).and.(i3.ge.0)).and.(i4.ge.0)).and.(i5.ge.0)) 1 .and.(i6.ge.0)) then u=factorial(i1)*factorial(i2) u=u*factorial(i3)*factorial(i4) u=u*factorial(i5)*factorial(i6) sum=sum+msign/u endif 1 continue ar=ar*sum if(abs(ar).lt.1.0d-38) ar=0.0d0 endif rthrj=ar return end Function rthrj16(x1,x2,x3,y1,y2,y3) c {real three - J symbols} implicit none integer*4 nfac0 parameter (nfac0=170) real*16 f016,fc16,ar,msign,u,sum common/npar16/f016,fc16(2*nfac0) real*8 rthrj16,x1,x2,x3,y1,y2,y3 integer*4 j1d,j2d,j3d,m1d,m2d,m3d,s,ni,i1,i2,i3,i4,i5,i6 j1d=NINT(x1*2.0d0) j2d=NINT(x2*2.0d0) j3d=NINT(x3*2.0d0) m1d=NINT(y1*2.0d0) m2d=NINT(y2*2.0d0) m3d=NINT(y3*2.0d0) if(abs(m1d).gt.j1d.or.abs(m2d).gt.j2d.or.abs(m3d).gt.j3d.or. 1m1d+m2d .ne.-m3d.or.j1d+j2d.lt.j3d.or.j2d+j3d.lt.j1d.or. 2j3d+j1d.lt.j2d .or.mod(j1d-m1d,2).ne.0 .or. 3mod(j2d-m2d,2) .ne.0.or.mod(j3d-m3d,2).ne. 0.or. 4mod(j1d+j2d+j3d,2).ne.0) then ar=0.0q0 else if (mod((j1d-j2d-m3d)/2,2).eq.0) then msign=1.0q0 else msign=-1.0q0 endif s=NINT(x1+x2+x3) ar=fc16(s-j3d)*fc16(s-j2d)*fc16(s-j1d)/fc16(s+1) 1 *fc16((j1d+m1d) / 2)*fc16((j1d-m1d) / 2) 2 *fc16((j2d+m2d) / 2)*fc16((j2d-m2d) / 2) 3 *fc16((j3d+m3d) / 2)*fc16((j3d-m3d) / 2) ar=msign*sqrt(ar) sum=0.0q0 do 1 ni=0 , s i1=ni i2=(j1d+j2d-j3d) / 2-ni i3=(j1d-m1d) / 2-ni i4=(j2d+m2d) / 2-ni i5=(j3d-j2d+m1d) / 2+ni i6=(j3d-j1d-m2d) / 2+ni if( mod(ni, 2) .ne. 0) then msign=-1.0q0 else msign=1.0q0 endif if(((((i2.ge.0).and.(i3.ge.0)).and.(i4.ge.0)).and.(i5.ge.0)) 1 .and.(i6.ge.0)) then u=fc16(i1)*fc16(i2)*fc16(i3)*fc16(i4)*fc16(i5)*fc16(i6) sum=sum+msign/u endif 1 continue ar=ar*sum if(abs(ar).lt.1.0q-38) ar=0.0q0 endif rthrj16=dble(ar) return end c ============================================================ subroutine ux(sri,srr,sli,slr,J,K,M,JP,KP,MP,ascr,asci) c uxl/r= c c sum(J*,K*) (-1)^(K*+K+M+M') (2J*+1) [( J* 1 J) - (J* 1 J ) ] c (M'+/-1 -1 -M) (M'+/-1 1 -M ) c c (J' 1 J*) sum(q q') (J* 1 J ) ( J' 1 J*) al(q,q') c (M' +/-1 -M'-/+1) (K* q -K ) ( K' qp -K*) implicit none real*8 sr,J,K,M,JP,KP,MP,Js,Ks,js2,q,qp,sumr,sumi,si real*8 z0,one,two,asci(3,3),ascr(3,3),f,sli,slr real*8 sri,srr,fr,fl,rthrj integer*4 i,l data z0,one,two/0.0d0,1.0d0,2.0d0/ sri=z0 srr=z0 sli=z0 slr=z0 c J*=J-1...J+2, J*>0,|J*-J'|<=1 Js=dabs(J-one)-one c J* = J-1,1,J+1 1 Js=Js+one if(dabs(Js-JP).le.one)then js2=Js*(Js+one) sumr=z0 sumi=z0 Ks=K-two c K* = K-1,K,K+1,|K*-KP|<=1,|K*|0,|J*-J'|<=1 Js=dabs(J-one)-one c J* = J-1,1,J+1 1 Js=Js+one if(dabs(Js-JP).le.one)then js2=Js*(Js+one) sumr=z0 sumi=z0 Ks=K-two c K* = K-1,K,K+1,|K*-KP|<=1,|K*|(jd): do 103 id=1,Ndim c excited state ide = Cjke(jde,ide)*|JEKEME>(jde): do 103 ide=1,Ndime c transition moments, c = Cjke_jde Cjk_jd , etc if(lsym)then jdstart=id jdend=id jdestart=ide jdeend=ide else jdstart=1 jdend=Ndim jdestart=1 jdeend=Ndime endif c total moment F=J+S+I,-ground F=J-S-I,J,J+S+I c J-rotational c S-spin c I-nuclear c F=N+I=N-1...N+1 F=N-two 101 F=F+one if(F.gt.z0)then F2=F*(F+one) c ground state g-factor gF=gN*(F2+N2-two)/two/F2*bz*mb mf=-f-one c MF=-F...F 2 mf=mf+one c ground state energy: egg=ek(id)+gf*mf p=exp(-egg/tf/t)/qp FE=NE-two 106 FE=FE+one if(FE.gt.z0)then FE2=FE*(FE+one) c excited state g-factor gFE=gNE*(FE2+NE2-two)/two/FE2*bz*mb c MFE=-FE...FE,|MF-MFE|<=2 MFE=-fe-one 5 MFE=MFE+one if(dabs(MFE-mf).le.two)then c excited state energy for each vibration: eee=ww(im)+ek(ide)+gFE*MFE-egg ilr=0.0d0 ili=0.0d0 irr=0.0d0 iri=0.0d0 ilrdoc=0.0d0 ilidoc=0.0d0 irrdoc=0.0d0 iridoc=0.0d0 do 81 jd=jdstart,jdend K=-J+jd-1 ak=dble(K) c1=Cjk(jd,id) do 81 jde=jdestart,jdeend KE=-JE+jde-1 co=Cjke(jde,ide)*c1 if(abs(K-KE).le.2)then ake=dble(KE) c MN=MF-1..MF+1 MN=MF-two 107 MN=MN+one if(dabs(MN).le.N)then c cg1=rcg(N,MN,one,MF-MN,F,MF)*co c cg2=rcg(NE,MFE-MF+MN,one,MF-MN,FE,MFE)*cg1 c MJ=MN-SZ,MF+SZ MJ=MN-half-one 6 MJ=MJ+one if(dabs(MJ).le.J)then c cg3=rcg(aj,mj,half,MN-mj,N,MN)*cg2 c cg4=rcg(ajE,MJ+MFE-MF,half,MN-MJ,NE,MFE-MF+MN)*cg3 if(abs(cg2).gt.z0)then call ux(ri,rr,lr,li,aj,ak,MJ,aje,ake,MJ+MFE-MF,ascr,asci) ilr=ilr+cg4*lr ili=ili+cg4*li irr=irr+cg4*rr iri=iri+cg4*ri if(ldc)then call uxdoc(ridoc,rrdoc,lrdoc,lidoc,aj,ak,MJ, 1 aje,ake,MJ+MFE-MF,ascr,asci) ilrdoc=ilrdoc+cg4*lrdoc ilidoc=ilidoc+cg4*lidoc irrdoc=irrdoc+cg4*rrdoc iridoc=iridoc+cg4*ridoc endif endif endif if(MJ.lt.MN+half)goto 6 endif if(MN.lt.MF+one)goto 107 endif 81 continue ileft =J2*JE*(ilr**2+ili**2 ) iright=J2*JE*(irr**2+iri**2) rami=(iright+ileft)*p*1.0d9/eee roai=(iright-ileft)*p*1.0d9/eee if(ldc)then ileftdoc =J2*JE*(ilrdoc**2+ilidoc**2) irightdoc=J2*JE*(irrdoc**2+iridoc**2) ramidoc=irightdoc+ileftdoc if(ramidoc.ne.z0)then doci=(irightdoc-ileftdoc)/ramidoc else doci=z0 endif endif if(lwr)then if(imax.lt.rami)imax=rami if(rami.gt.thresh*imax)then ipeak=ipeak+1 write(12,1212)ipeak,eee,rami,roai,J,N,F,MF,JE,NE,FE,MFE 1212 format(i12,3f16.7,i3,3f4.0,' -> ',i3,3f4.0) endif endif if(rami.gt.1.0d-9)then call spread1(eee,rami,ram,npo,wmin,nt,delw) call spread1(eee,roai,roa,npo,wmin,nt,delw) if(ldc)call spread1(eee,doci,doc,npo,wmin,nt,delw) endif endif if(MFE.lt.FE)goto 5 endif if(FE.lt.NE+one)goto 106 if(MF.lt.F)goto 2 endif c F>0 if(F.lt.N+one)goto 101 103 continue c id,ide endif if(NE.lt.aje+half)goto 102 1001 deallocate(He,Cjke,ete,eke) c JE endif if(N.lt.aj+half)goto 1 call wrs('ramJ.prn',ram,wmin,wmax,npo) call wrs('roaJ.prn',roa,wmin,wmax,npo) if(ldc)call wrs('docJ.prn',doc,wmin,wmax,npo) if(lwr)close(12) 100 deallocate(H,Cjk,et,ek) c J return end c ============================================================ function ggg(gJ,gI,F2,J2,I2) implicit none real*8 ggg,gJ,gI,F2,J2,I2 if(F2.eq.0)then ggg=0.0d0 else ggg=(gJ*(F2+J2-I2)+gI*(F2+I2-J2))/(F2+F2) endif return end c ============================================================