program mroa write(6,600) 600 format(' Simulation of rotationally resolved magnetic circular', 1 /, ' dichroism and Raman optical activity of NO',/,/, 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 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,c,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,c,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 common/ipar/jmin,jmax real*8,allocatable::d(:),rs(:),ram(:),roa(:) character*1 branch,fstate common/spar/branch,fstate real*8 na,hprime,cvl,r,mb,m0,Debye common/rpar/na,hprime,cvl,r,mb,m0,Debye logical lwr,lzmat,lno,lno2 common/lpar/lwr,lzmat,lno,lno2 integer*4, allocatable::ntw(:) real*8, allocatable::alphad(:,:,:,:),Gpd(:,:,:,:),ad(:,:,:,:,:), 1tw(:,:),al(:,:,:),asi(:,:,:),asr(:,:,:),aii(:,:,:),air(:,:,:) 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),al(nm,3,3),asr(nm,3,5),asi(nm,3,5),ntw(1), 2 aii(nm,3,3),air(nm,3,3)) 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 6 ix=1,3 do 6 jx=1,3 6 al(im,ix,jx)=alphad(1,im,ix,jx) 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,debye) lsym=dabs(AA-BB).lt.rtol 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,')') if(lno2)then call rrno2(Jmin,Jmax,lsym,AA,BB,CC,AE,BE,CE,nm,air,aii, 1 wmin,wmax,npo,tw) else call rrspec(Jmin,Jmax,lsym,AA,BB,CC,AE,BE,CE,nm,air,aii, 1 wmin,wmax,npo,tw) 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 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*|(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 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 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 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 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) 100 deallocate(H,Cjk,et,ek) c J 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,c,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,c,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth real*8 na,hprime,cvl,r,mb,m0,Debye common/rpar/na,hprime,cvl,r,mb,m0,Debye 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,c,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,c,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,c,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,c,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 subroutine Initfac implicit none integer*4 nfac0 parameter (nfac0=100) real*8 f0,factorial common/npar/f0,factorial(nfac0) real*8 fac integer*4 i do 1 i=0 , nfac0 1 factorial(i)=fac(i) return end subroutine Setparameters implicit none real*8 d0,hw,t,c,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth common/nopar/d0,hw,t,c,bcm,becm,acm,aecm,bz,l, 1tcm,dsqrtpi,pi,y,ye,bandwidth real*8 na,hprime,cvl,r,mb,m0,Debye common/rpar/na,hprime,cvl,r,mb,m0,Debye character*1 branch,fstate common/spar/branch,fstate integer*4 jmax,jmin common/ipar/jmin,jmax real*8 wmin,wmax integer*4 npo common/opts/wmin,wmax,npo character*10 s10 logical lwr,lzmat,lno,lno2 common/lpar/lwr,lzmat,lno,lno2 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} m0=1.2566370d-6 c {Js^2C^-2m^-1} {vacuum permeativity} Debye=3.335445d-30 c {C.m} pi=4.0d0*datan(1.0d0) dsqrtpi=dsqrt(pi) jmax=10 jmin=0 t=298.0d0 bz=2.0d0 bandwidth=0.3d0 inquire(file='NO.PAR',exist=lno) inquire(file='NO2.PAR',exist=lno2) if(lno)then 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,*)c 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 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 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.'BAND')read(9,*)bandwidth if(s10(1:2).eq.'BZ')read(9,*)bz goto 1 2 close(9) tcm=t*r/na/hprime/2.0d0/pi/cvl/100.d0 write(6,6000)npo,wmin,wmax,jmin,jmax,bandwidth,t,bz,lwr 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) return end Function rthrj(x1,x2,x3,y1,y2,y3) c {real three - J symbols} implicit none integer*4 nfac0 parameter (nfac0=100) 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 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 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,debye) 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 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(4,*)'Dalton output' write(4,*)nat do 343 l=1,nat 343 if(iz(l).gt.0)write(4,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 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 ge(i,j)=-4.0d0*xp(i,j)/t(i,i)*mp 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 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 do 9 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) 9 gq(i,j)=qtot(i,j)/t(i,i)*mp 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