program rcd c extracts from G.OUT molecular parameters c and calculates rotationa spectra and circular dichroism implicit none integer*4 i,IERR,j,k,id,jd,Jmax,Jmin,Ndim,Ndime,NPOINT, 1JE,KE,ide,jde,jdstart,jdend,jdestart,jdeend,io,lwr,units, 1JEMIN,JEMAX,f,ib,ic,JPP,KPP,M,ME,nbuf,nba real*8 ABC(6),g(3,3),p,qp,c,Te,gcm, 1DD,wo,WMAX,WMIN,WEXC,u(3),zd,zr,dw,tf,kcm,pi,rtol,ak,j2,zq, 1spi,cjje,coselJ,uN,debye,jkp,jkm,mfac,mf,z0,uer,uei,mer,mei, 1ur(2),ui(2),uq(3,3),sgr(3,3),sgi(3,3),cp,CM,dsum,rsum,qsum, 1ckr,cki,er,ei,uetr,ueti,erp,eip,w, 1p1r,p1i,ckbr,ckbi,ckcr,ckci,pG,pE,ar(3),ai(3),br,bi, 1mr(2),mi(2),qr(2),qi(2),tolzd,cop,rx,ry,dx,dy,qx,qy real*8,allocatable:: H(:,:),Cjk(:,:),et(:),ek(:),He(:,:), 1Cjke(:,:),ete(:),eke(:),sd(:),sr(:),zrb(:,:) integer*4,allocatable::irb(:,:) LOGICAL lzmat,lg,lsym,lquad,ltab,lnum character*80 fil write(6,6000) 6000 format(/,' Rotational circular dichroism intensities',/,/, 1 ' Input: RCD.PAR options',/,/, 1 ' G.OUT Dalton/Gaussian output',/, 1 ' Output: RCD.TAB line spectrum',/, 1 ' LARGE.TAB line spectrum, selection',/, 1 ' DETAIL/N.TAB details',/, 1 ' D/N.PRN absorption',/, 1 ' R/N.PRN RCD',/,/) call setpar(c,kcm,lzmat,Jmin,Jmax,LG,lwr,ABC,dw,Te,WEXC,WMAX, 1WMIN,DD,NPOINT,u,tf,rtol,JEMIN,JEMAX,uN,debye,pi,spi,uq,lquad, 1ltab,nbuf,fil,lnum,units) CM=219474.63d0 tolzd=1.0d-10 z0=0.0d0 gcm=0.0333564095198152d0 call lookatout(fil,lzmat,g,ABC,u,uq) c write some selected strengths: call wrinv(g,u,uN,debye) lsym=dabs(ABC(1)-ABC(2)).lt.rtol write(6,60000)(ABC(i),i=1,3),(ABC(i),i=4,6), 1JMIN,JMAX,WEXC,WMIN,WMAX, 1NPOINT,LG,Te,LWR,nbuf,u(1),u(2),u(3),((uq(i,j),j=1,3),i=1,3), 1lquad,lnum,lsym,rtol 60000 format(' A B C ground :', 3f12.6,' GHz',/, 1 ' A B C excited :', 3f12.6,' GHz',/, 1 ' Jmin Jmax :',2i12,/, 1 ' WEXC WMIN WMAX :',3F12.6,' GHz',/, 1 ' Npoints LG :',i12,l12,/, 1 ' Temperature :',f12.6,' K',/, 1 ' Large output :',i12,/, 1 ' Buffer size :',i12,/, 1 ' ux uy uz :',3f12.6,/, 1 ' Qxx Qxy Qxz :',3f12.6,/, 1 ' Qyx Qyy Qyz :',3f12.6,/, 1 ' Qzx Qzy Qzz :',3f12.6,/, 1 ' Quadrupolar part:',l12,/, 1 ' Nymerical path :',l12,/, 1 ' Symmetric top :',l12,' (tol:',e8.2,')',/) allocate(sd(NPOINT),sr(NPOINT),zrb(3,nbuf),irb(6,nbuf)) sr=z0 sd=z0 io=0 nba=0 zrb(1,1)=z0 irb(1,1)=0 c g - first index mechanic, second magnetic c spherical sg-tensor, first index mechanic spherical, c second magnetic Cartesian do 1 f=1,3 c - = x/2 - iy/2: sgr(1,f)= g(1,f)/2.0d0 sgi(1,f)=-g(2,f)/2.0d0 c 0 = z: sgr(2,f)= g(3,f) sgi(2,f)= z0 c + = x/2 + iy/2: sgr(3,f)= g(1,f)/2.0d0 1 sgi(3,f)= g(2,f)/2.0d0 c call initab(lwr,lnum,ltab,units) c calculate partititon function: call partf(Jmin,Jmax,qp,ABC,Te,tf,pi,lsym) do 2 J=Jmin,Jmax Ndim=2*J+1 write(6,609)J 609 format(i3,$) j2=dble(Ndim) 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,ABC(1),ABC(2),ABC(3),Ndim) else call seth(h,J,ABC(1),ABC(2),ABC(3)) call TRED12(Ndim,H,Cjk,ek,2,IERR,et) endif if(lwr.gt.1)then write(6,6002)J 6002 format(' J = ',I4,', energies/GHz:') call wre(ek,Ndim,Te,tf,qp,z0) call wrm(H,Ndim) endif do 3 JE=max(0,J-1),J+1 Ndime=2*JE+1 cjje=coselJ(J,JE) mf=mfac(J,JE) do 710 f=1,2 do 710 M=-J,J ME=M+1 call coselM(J ,M ,f,JE,ME,er ,ei ) 710 call coselM(JE,ME,f,J ,M ,erp,eip) 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,ABC(4),ABC(5),ABC(6),Ndime) else call seth(He,JE,ABC(4),ABC(5),ABC(6)) call TRED12(Ndime,He,Cjke,eke,2,IERR,ete) if(lwr.gt.1)then write(6,601)JE 601 format(' JE = ',I4,' energies/GHz:') call wre(eke,Ndime,tf,Te,qp,WEXC) call wrm(He,Ndime) endif endif c ground state id = Cjk(jd,id)*|JKM>(jd): do 4 id=1,Ndim pG=exp(-ek(id)/tf/Te)/qp c excited state ide = Cjke(jde,ide)*|JEKEME>(jde): do 4 ide=1,Ndime c transition energy |G> -> |E> in GHz: wo=WEXC+eke(ide)-ek(id) c transition energy |G> -> |E> in au: w=wo*gcm/CM c transition moments, c = Cjke_jde Cjk_jd , etc if(wo.gt.z0)then pE=exp(-eke(ide)/tf/Te)/qp p=pG-pE if(lsym)then jdstart=id jdend=id jdestart=ide jdeend=ide else jdstart=1 jdend=Ndim jdestart=1 jdeend=Ndime endif c |G=id> ->|E=ide> if(ltab)write(67,603)J,id,JE,ide 603 format(' |J i =',2i3,'> -> |JE ie =',2i3,'>') c total reduced el dipole uer=z0 uei=z0 c total reduced el dipole uetr=z0 ueti=z0 c total reduced mag dipole mer=z0 mei=z0 do 51 jd=jdstart,jdend K=-J+jd-1 ak=dble(k) c J+|J K>-sqrt(J(J+1)-K(K+1)) |J K+1> c J-|J K>-sqrt(J(J+1)-K(K-1)) |J K+1> jkp=dsqrt(dble(J*(J+1)-K*(K+1))) jkm=dsqrt(dble(J*(J+1)-K*(K-1))) do 51 jde=jdestart,jdeend KE=-JE+jde-1 cp=Cjke(jde,ide)*Cjk(jd,id)*cjje if(abs(K-KE).le.1)then do 91 f=1,3 call coselK(J ,K ,f,JE,KE,ckr,cki) c =sum(f) CC (JJ') uf (J K f JE KE): uer =uer +cp*u(f)*ckr uei =uei +cp*u(f)*cki call coselK(JE ,KE ,f,J,K,ckr,cki) c =sum(f) CC (JJ') uf (JE KE f J K): uetr=uetr+cp*u(f)*ckr 91 ueti=ueti+cp*u(f)*cki endif if(abs(K-KE).le.2)then c A+ = sum(f) (JE KE f J K+1) gf- = A(1): c A0 = sum(f) (JE KE f J K ) gf0 = A(2): c A- = sum(f) (JE KE f J K-1) gf+ = A(3): ar=z0 ai=z0 do 71 ic=1,3 do 71 f=1,3 call coselK(JE,KE ,f,J,K+2-ic,ckr,cki) ar(ic)=ar(ic)+ckr*sgr(ic,f)-cki*sgi(ic,f) 71 ai(ic)=ai(ic)+ckr*sgi(ic,f)+cki*sgr(ic,f) do 101 f=1,3 mer=mer+uN*cp*(ar(1)*jkp+ar(3)*jkm+ar(2)*ak) 101 mei=mei+uN*cp*(ai(1)*jkp+ai(3)*jkm+ai(2)*ak) endif 51 continue c jd,jde = K = (JMaJEME) c mag dipole = (JEMEaJM) a=X,Y: do 61 f=1,2 call coselM(J,M,f,JE,ME,er,ei) ur(f)=uer*er-uei*ei ui(f)=uer*ei+uei*er call coselM(JE,ME,f,J,M,er,ei) mr(f)=mer*er-mei*ei 61 mi(f)=mer*ei+mei*er qr=z0 qi=z0 if(lquad)then c quadrupole ,A=XZ,YZ: do 5 jd=jdstart,jdend K=-J+jd-1 ak=dble(k) do 5 jde=jdestart,jdeend KE=-JE+jde-1 if(abs(K-KE).le.2)then do 11 JPP=max(0,J-1),J+1 c cc (JEJPP)(JPP,J): cop=Cjke(jde,ide)*Cjk(jd,id)*coselJ(JE,JPP)*coselJ(JPP,J) c B = sum(b,c,KPP) Qbc (JEKEbJPPKPP)(JPPKPPcJK) br=z0 bi=z0 do 111 KPP=K-1,K+1 do 111 ib=1,3 call coselK(JE,KE,ib,JPP,KPP,ckbr,ckbi) do 111 ic=1,3 call coselK(JPP,KPP,ic,J,K,ckcr,ckci) br=br+uq(ib,ic)*(ckbr*ckcr-ckbi*ckci) 111 bi=bi+uq(ib,ic)*(ckbi*ckcr+ckbr*ckci) call coselM(JPP,M,3,J,M ,er,ei) c B (JPPMZJM): p1r=br*er-bi*ei p1i=bi*er+br*ei do 11 f=1,2 call coselM(JE,ME,f ,JPP,M ,er,ei) c w cc (JEJPP)(JPP,J) B (JPPMZJM) (JEME f JPPM): qr(f)=qr(f)+w*cop*(p1r*er-p1i*ei) 11 qi(f)=qi(f)+w*cop*(p1i*er+p1r*ei) endif 5 continue c jd,jde endif dx=ur(1)*ur(1)+ui(1)*ui(1) dy=ur(2)*ur(2)+ui(2)*ui(2) zd=1.5d0*(dx+dy)/j2 if(zd.gt.tolzd)then dsum=dsum+zd rx=ur(1)*mi(1)+ui(1)*mr(1) ry=ur(2)*mi(2)+ui(2)*mr(2) zr=1.5d0*(rx+ry)/j2 rsum=rsum+zr call smooth(tolzd,qr(1),qr(2),qi(1),qi(2)) c Q = - 1.5 w Re (ux qyz - uy qxz) qx= ur(1)*qr(2)-ui(1)*qi(2) qy= - ur(2)*qr(1)+ui(2)*qi(1) zq=-1.5d0*(qx+qy)/j2 qsum=qsum+zq if(ltab)then write(67,602)J,id,M,JE,ide,ME 602 format(' |J i M = ',3i3,'> -> |JE ie ME = ',3i3,'>') write(67,672)ur(1),ui(1),ur(2),ui(2),mr(1),mi(1),mr(2),mi(2), 1 qr(1),qi(1),qr(2),qi(2),dx,dy,dx+dy,zd,rx,ry,rx+ry, 1 qx,qy,qx+qy,zr+zq 672 format(' el dipole X Y = ',4e12.4,/, 1 ' mag dipole X Y = ',4e12.4,/, 1 ' el quadru wXZ wYZ = ',4e12.4,/, 1 ' D X Y tot = ',3e12.4,/, 1 ' D ave = ',1e12.4,/, 1 ' Rm X Y tot = ',3e12.4,/, 1 ' Rq X Y tot = ',3e12.4,/, 1 ' R ave = ',1e12.4,/) endif call recordtrans(nbuf,nba,J,JE,id,ide,M,ME,wo,zd,zr,zq, 1 zrb,irb,debye,lwr,io,pG,pE,p,wmin,dw,dd,npoint,LG,sd,sr, 2 spi,pi,units) endif 6 continue c M,ME c NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN else c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA c c dx = sum (JMXJEM+1)(JEM+1XJM) c is not complex conjugate of !! c because (JMXJEM+1) is real, we can take the real part here: dx=(uer*uetr-uei*ueti)*mf*2.0d0 c 2.0d0 above for ME=M+1 and ME=M-1 dy=dx c D = 1.5 Re ( + )/(2J+1) zd=1.5d0*(dx+dy)/j2 c because (JMXJEM+1) is real, we can take the imaginary part here: rx=(uei*mer+uer*mei)*mf*2.0d0 c 2.0d0 above for ME=M+1 and ME=M-1 ry=rx c R = 1.5 Im ( + ) zr=1.5d0*(rx+ry)/j2 dsum=zd rsum=zr if(zd.gt.tolzd)then call recordtrans(nbuf,nba,J,JE,id,ide,0,0,wo,zd,zr,z0, 1 zrb,irb,debye,lwr,io,pG,pE,p,wmin,dw,dd,npoint,LG,sd,sr, 2 spi,pi,units) endif c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA endif c lnum if(ltab)write(67,673)uer,uei,uetr,ueti, 1 mer,mei,dsum,rsum,qsum 673 format(' ',60(1H*),/, 1 ' , = ',4e12.4,/, 1 ' = ',2e12.4,/, 1 ' = ',1e12.4,/, 1 ' = ',1e12.4,/, 1 ' = ',1e12.4,/) endif c wo>0 4 continue c id,ide 3 deallocate(He,Cjke,ete,eke) c JE if(mod(J,23).eq.0)write(6,*) 2 deallocate(H,Cjk,et,ek) c J write(6,*) call writeb(nbuf,nba,zrb,irb,units) if(lwr.gt.0)then write(66,660) close(660) 660 format(60(1h-)) write(6,*)'RCD.TAB' endif if(ltab)close(67) if(lnum)then call wrspec('DN.PRN',npoint,wmin,dw,sd,units) call wrspec('RN.PRN',npoint,wmin,dw,sr,units) else call wrspec('D.PRN',npoint,wmin,dw,sd,units) call wrspec('R.PRN',npoint,wmin,dw,sr,units) endif end c ================================================================== subroutine wrspec(sn,nw,wmin,dw,spec,units) implicit none character*(*) sn integer*4 nw,iw,units real*8 wmin,dw,spec(*),w,f,gcm c 1GHz = 0.033 cm-1 gcm=0.0333564095198152d0 if(units.eq.0)then f=1.0d0 elseif(units.eq.1)then f=1.0d3 elseif(units.eq.2)then f=gcm else call report('unknown units requested') endif open(8,file=sn) w=wmin-dw do 1 iw=1,nw w=w+dw 1 write(8,805)w*f,spec(iw) 805 format(f18.8,' ',g18.6) close(8) write(6,*)sn,' written' return end c ================================================================== subroutine trt(g,a) implicit none real*8 g(3,3),a(3,3),t(3,3) integer*4 i,j do 1 i=1,3 do 1 j=1,3 1 t(i,j)=g(1,j)*a(1,i)+g(2,j)*a(2,i)+g(3,j)*a(3,i) do 2 i=1,3 do 2 j=1,3 2 g(i,j)=t(i,1)*a(1,j)+t(i,2)*a(2,j)+t(i,3)*a(3,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 SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) 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) IMPLICIT INTEGER*4(I-N) 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 subroutine wrm(H,N) implicit none real*8 H(N,N) integer*4 i,j,N do 1 i=1,N 1 write(6,60)(H(i,j),j=1,N) 60 format(10f8.2) return end subroutine lookatout(fil,lzmat,gt,ABC,ud,qt) 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,ic character*80 s80 character*160 FN character*(*) fil logical lzmat real*8 bohr,r(3,n0),ama(n0),xp(3,3),qt(3,3),qtot(3,3),debye, 1ABC(6),c,g(3,3),cm(3),mt,t(3,3),q(3,3),a(3,3),z(3),u,mp,cha, 1gn(3,3),ge(3,3),gq(3,3),amu,e(3),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/ debye=2.541746913d0 bohr=0.52917721092d0 mp=1.00727646688d0 c=2.99792458d10 u=219474.63068d0 amu=1822.88839d0 cha=0.0d0 OPEN(2,FILE=fil,status='old') NG=0 1 READ(2,2000,END=1000)FN 2000 FORMAT(A160) 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 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 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 ENDIF c Gaussian: IF(FN(2:9).eq.'Charge= ')then write(6,2000)FN read(FN,208)cha 208 format(20x,f8.4) 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(27:55).eq.'Total quadrupole moments (au)')then 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 c gaussian traceless quadrupole moment = (2/3) x Dalton quadrupole moment: qtr(i,j)=qtr(i,j)*2.0d0/3.0d0 1211 qt(i,j)=qtr(i,j) c many possibilities to make diagonal elements from traceless, c choose xx=0: qt(1,1)=0.0d0 qt(2,2)= qtr(2,2)-qtr(1,1) qt(3,3)= qtr(3,3)-qtr(1,1) endif c Gaussian, Debye*Ang: IF(FN(2:43).eq.'Quadrupole moment (field-independent basis')then 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 qt=sum(i)qi ri ri 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(2,2))/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 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(ABC(1).eq.0.0d0)then IF(FN(2:28).eq.'Rotational constants (GHZ):')then c Gaussian: write(6,2000)FN read(FN(30:len(FN)),*)(ABC(I),I=1,3) ABC(4)=ABC(1) ABC(5)=ABC(2) ABC(6)=ABC(3) ELSEIF(FN(2:21).eq.'Rotational constants')then c Dalton: 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 ABC(4)=ABC(1) ABC(5)=ABC(2) ABC(6)=ABC(3) else backspace 2 endif endif 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 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 write(6,*)'mass center (A):' write(6,6001)(cm(i)*bohr,i=1,3) c subtract mass center: do 22 ia=1,nat do 22 i=1,3 22 r(i,ia)=r(i,ia)-CM(i) write(6,*)'coordinates shifted' do 23 i=1,3 23 ud(i)=ud(i)-cha*CM(i) write(6,*)'dipole shifted' do 24 i=1,3 do 24 j=1,3 24 qt(i,j)=qt(i,j)-CM(i)*ud(j)-CM(j)*ud(i) write(6,*)'quadrupole shifted' do 25 i=1,3 do 25 j=1,3 25 qt(i,j)=qt(i,j)/2.0d0 write(6,*)'quadrupole redefined as (1/2)qrr' c c inertia t and nuclear quadrupole q 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) c check that we have a right-handed system if(A(3,3)*(A(1,1)*A(2,2)-A(1,2)*A(2,1)).lt.0.0d0)then A(1,3)=-A(1,3) A(2,3)=-A(2,3) A(3,3)=-A(3,3) write(6,*)'chirality fixed' endif c make inertia-like electric quadrupole moment: 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 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 99 ic=1,2 if(ic.eq.1)then write(6,*)'Laboratory system:' else write(6,*)'Trafo:' write(6,6001)((a(i,j),i=1,3),j=1,3) call trt(xp,a) call trt(t,a) call trt(q,a) call trt(g,a) call trt(qtot,a) call trt(qtr,a) call trt(qt,a) call trv(ud,a) write(6,*) write(6,*)'Molecular System:' write(6,*)'*****************' endif 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,*)'xp:' write(6,6001)((xp(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,*)'Quadrupole moment nuclear:' write(6,6001)((q(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 (1/2) qrr:' write(6,6001)((qt(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,*)'g-tensor from quadrupole:' write(6,6001)((gq(i,j),i=1,3),j=1,3) write(6,*)'dipole (au):' 99 write(6,6001)(ud(i),i=1,3) 6001 format(/,3(3f18.6,/),/) return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine setpar(c,kcm,lzmat,Jmin,Jmax,LG,lwr,ABC,dw, 1Te,WEXC,WMAX,WMIN,DD,NPOINT,u,tf, 1rtol,JEMIN,JEMAX,uN,debye,pi,spi,uq,lquad,ltab,nbuf,fil,lnum, 1units) implicit none real*8 c,kcm,ABC(6),Te,WEXC,WMAX,WMIN,u(3), 1DD,tf,dw,rtol,uN,debye,pi,spi,uq(3,3) integer*4 Jmin,Jmax,NPOINT,lwr,JEMIN,JEMAX,i,j,nbuf,units logical lzmat,LG,lex,lquad,ltab,lnum character*80 FN,fil pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) debye=2.541746913d0 c nuclear magneton in au (e hbar/2/Mp): uN=0.5d0/1836.152663302331d0 c tolerance for symmetric top rotational cosntant differences: rtol=1.0d-6 c=2.99792458d10 c kelvin to cm-1: kcm=0.695d0 c kelvin to GHz: tf=2.083674d1 lzmat=.true. c rotational J-number limits: Jmax=7 Jmin=0 Jemin=0 Jemax=7 c gaussian/lorentzian peaks: LG=.true. c extensive output: lwr=0 c rotational constants of ground and excited states: ABC=0.0d0 c temperature: Te=293.0d0 c central frequency: WEXC=0.0d0 c frequency limits: WMAX=200.0d0 WMIN=-200.0d0 c peak width: DD=1.0d0 NPOINT=4001 c electric dipole moment: u(1)=0.0d0 u(2)=0.0d0 u(3)=1.0d0 c electric quadrupole moment: uq=0.0d0 c use the quadrupole moment: lquad=.false. c evaluate numerically all M transitions: lnum=.false. c write tables for all transitions: ltab=.false. c buffer size for the biggest transitions: nbuf=1000 c default file name: fil='G.OUT' c units ... plot spectrum in units c units = 0 ... GHz c 1 ... MHz c 2 ... cm-1 units=0 inquire(file='RCD.PAR',exist=lex) if(lex)then write(6,*)'RCD.PAR' open(55,file='RCD.PAR') 551 read(55,2000,end=555,err=555)FN 2000 format(a80) if(FN(1:2).eq.'AA')read(55,*)ABC(1) if(FN(1:2).eq.'BB')read(55,*)ABC(2) if(FN(1:2).eq.'CC')read(55,*)ABC(3) if(FN(1:2).eq.'AE')read(55,*)ABC(4) if(FN(1:2).eq.'BE')read(55,*)ABC(5) if(FN(1:2).eq.'CE')read(55,*)ABC(6) if(FN(1:2).eq.'DD')read(55,*)DD if(FN(1:4).eq.'FILE')read(55,*)fil if(FN(1:5).eq.'JEMAX')read(55,*)JEMAX if(FN(1:5).eq.'JEMIN')read(55,*)JEMIN if(FN(1:4).eq.'JMAX')read(55,*)JMAX if(FN(1:4).eq.'JMIN')read(55,*)JMIN if(FN(1:4).eq.'GAUS')read(55,*)LG if(FN(1:3).eq.'LWR')read(55,*)lwr if(FN(1:4).eq.'LNUM')read(55,*)lnum if(FN(1:4).eq.'LQUA')read(55,*)lquad if(FN(1:3).eq.'LZM')read(55,*)lzmat if(FN(1:4).eq.'LTAB')read(55,*)ltab if(FN(1:4).eq.'NBUF')read(55,*)nbuf if(FN(1:4).eq.'NPOI')read(55,*)NPOINT if(FN(1:4).eq.'QXYZ')read(55,*)((uq(i,j),j=1,3),i=1,3) if(FN(1:4).eq.'RTOL')read(55,*)rtol if(FN(1:4).eq.'TEMP')read(55,*)Te if(FN(1:4).eq.'UNIT')read(55,*)units if(FN(1:4).eq.'UXYZ')read(55,*)(u(i),i=1,3) if(FN(1:4).eq.'WEXC')read(55,*)WEXC if(FN(1:4).eq.'WMAX')read(55,*)WMAX if(FN(1:4).eq.'WMIN')read(55,*)WMIN goto 551 555 close(55) endif dw=(WMAX-WMIN)/dble(NPOINT) return end 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 subroutine wre(e,n,tf,Te,qp,WEXC) implicit none real*8 e(*),tf,Te,qp,WEXC,p integer*4 i,n do 1 i=1,n p=exp(-(e(i)+WEXC)/tf/Te)/qp*dble(n) 1 write(6,600)i,e(i),' GHz',p*100.0d0 600 format(i6,f12.6,a4,f12.5,'%') return end 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 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 subroutine wrinv(g,u,uN,debye) implicit none real*8 g(3,3),u(3),uN,debye,r,d character*1 xyz(3) integer i,j,k data xyz/'x','y','z'/ do 1 j=1,3 k=j+1 if(k.gt.3)k=1 i=k+1 if(i.gt.3)i=1 r=uN/2.0d0*debye**2*u(i)*(g(j,k)-g(k,j))*1.0d4 d=u(i)**2*debye**2 write(6,600)xyz(i),xyz(j),xyz(k),xyz(k),xyz(j),r,xyz(i),d 600 format(' u0/2*u',a1,'(g',2a1,' - g',2a1,') = r = ',f12.6,/, 1 ' u(',a1,')^2 = d = ',f12.6) r=uN/4.0d0*debye**2*u(i)*(g(j,k)+g(k,j))*1.0d4 d=u(i)**2*debye**2/2.0d0 write(6,610)xyz(i),xyz(j),xyz(k),xyz(k),xyz(j),r,xyz(i),d 610 format(' u0/4*u',a1,'(g',2a1,' + g',2a1,') = r = ',f12.6,/, 1 ' u(',a1,')^2 / 2 = d = ',f12.6) 1 write(6,*) return end function mfac(J,JE) c sum(M=-J...J) |(J M CXa J' M+1)|**2=sum |(J M CYa J' M+1)|**2 implicit none real*8 mfac integer*4 J,JE if(J.eq.JE)then mfac=dble(2*J*(J+1)*(2*J+1))/3.0d0 else if(J+1.eq.JE)then mfac=-dble(2*(J*(4*J*(J+3)+11)+3))/3.0d0 else if(J-1.eq.JE)then mfac=-dble(2*J*(2*J-1)*(2*J+1))/3.0d0 else mfac=0.0d0 endif endif endif return end subroutine coselM(J,M,f,JE,ME,jmr,jmi) implicit none integer*4 J,M,F,JE,ME,s real*8 jmr,jmi,z jmr=0.0d0 jmi=0.0d0 s=ME-M if(f.eq.3)then c : if(M.eq.ME)then if(JE.eq.J+1)jmr= 2.0d0*dsqrt(dble(JE**2-M**2)) if(JE.eq.J-1)jmr=-2.0d0*dsqrt(dble(J **2-M**2)) if(JE.eq.J) jmr=dble(M+M) endif endif z=dble(s) if(abs(s).ne.1)return c and : if(f.eq.1)then if(JE.eq.J+1)jmr=z*dsqrt(dble((JE+s*M)*(JE+1+s*M))) if(JE.eq.J-1)jmr=z*dsqrt(dble((J -s*M)*(J -1-s*M))) if(JE.eq.J) jmr=-dsqrt(dble(J*(J+1)-M*(M+s))) endif if(f.eq.2)then if(JE.eq.J+1)jmi= dsqrt(dble((JE+s*M)*(JE+1+s*M))) if(JE.eq.J-1)jmi= dsqrt(dble((J -s*M)*(J -1-s*M))) if(JE.eq.J) jmi=-z*dsqrt(dble(J*(J+1)-M*(M+s))) endif return end subroutine coselK(J,K,f,JE,KE,jkr,jki) implicit none integer*4 J,K,F,JE,KE,s real*8 jkr,jki,z jkr=0.0d0 jki=0.0d0 if(f.eq.3)then c : if(K.eq.KE)then if(JE.eq.J+1)jkr= 2.0d0*dsqrt(dble(JE**2-K**2)) if(JE.eq.J-1)jkr=-2.0d0*dsqrt(dble(J **2-K**2)) if(JE.eq.J) jkr=dble(K+K) endif endif s=KE-K z=dble(s) if(abs(s).ne.1)return c and : if(f.eq.1)then if(JE.eq.J+1)jkr=z*dsqrt(dble((JE+s*K)*(JE+1+s*K))) if(JE.eq.J-1)jkr=z*dsqrt(dble((J -s*K)*(J -1-s*K))) if(JE.eq.J) jkr=-dsqrt(dble(J*(J+1)-K*(K+s))) endif if(f.eq.2)then if(JE.eq.J+1)jki=-dsqrt(dble((JE+s*K)*(JE+1+s*K))) if(JE.eq.J-1)jki=-dsqrt(dble((J -s*K)*(J -1-s*K))) if(JE.eq.J) jki= z*dsqrt(dble(J*(J+1)-K*(K+s))) endif return end function coselJ(J,JE) implicit none integer*4 J,JE real*8 coselJ coselJ=0.0d0 if(JE.eq.J+1)coselJ=1.0d0/(dble(4*JE)*dsqrt(dble(4*JE**2-1))) if(JE.eq.J-1)coselJ=1.0d0/(dble(4*J )*dsqrt(dble(4*J **2-1))) if(JE.eq.J.and.J.ne.0)coselJ=1.0d0/dble(4*J*(J+1)) return end subroutine saveb(nbuf,nba,J,JE,jd,jde,M,ME,wo,zd,zr,zrb,irb) implicit none integer*4 J,JE,jd,jde,nba,nbuf,irb(6,nbuf),i,n,nend,k,M,ME,l real*8 zd,zr,zrb(3,nbuf),a,tol,wo tol=1.0d-10 a=dabs(zr) if(a.lt.tol)return c nba .. number of non-zero elements if(nba.eq.0)then nend=1 n=1 goto 99 endif do 1 i=1,nba if(a.gt.dabs(zrb(3,i)))then nend=min(nbuf,nba+1) do 2 k=nend,i+1,-1 zrb(1,k)=zrb(1,k-1) zrb(2,k)=zrb(2,k-1) zrb(3,k)=zrb(3,k-1) do 2 l=1,6 2 irb(l,k)=irb(l,k-1) n=i goto 99 endif 1 continue if(nba.lt.nbuf)then nend=nba+1 n=nend goto 99 else return endif 99 nba=nend zrb(1,n)=wo zrb(2,n)=zd zrb(3,n)=zr irb(1,n)=J irb(2,n)=jd irb(3,n)=JE irb(4,n)=jde irb(5,n)=M irb(6,n)=ME return end subroutine writeb(nbuf,nba,zrb,irb,units) implicit none integer*4 nba,nbuf,irb(6,nbuf),i,units,j real*8 zrb(3,nbuf) character*4 lu(3) data lu/' GHz',' MHz',' cm-1'/ open(46,file='LARGE.TAB') write(46,662)lu(units+1) 662 format(' freq D R 4R/D', 1' J j Je k M ME Probab',/, 1 'transition ',a4,' debye^2 10^-4debye^2') write(46,660) 660 format(60(1h-)) do 1 i=1,nba 1 write(46,661)i,(zrb(j,i),j=1,3),zrb(3,i)/zrb(2,i), 1(irb(j,i),j=1,6),1.0d0 661 format(i7,f11.4,3f10.3,6i3,f9.5) write(46,660) close(46) return end subroutine sprd(zd,zr,p,wo,wmin,dw,dd,npoint,LG,sd,sr,spi,pi) implicit none real*8 zd,zr,p,si,wo,wmin,dw,w,du,dd,sr(*),sd(*),ss,spi,pi integer*4 iu,is,i,npoint logical LG zd=p*zd zr=p*zr do 901 iu=1,2 if(iu.eq.1)then is=-1 si=-1.0d0 i=int((wo-wmin)/dw)-1 w=wmin+dw*dble(i-1)-dw else is=1 si=1.0d0 i=int((wo-wmin)/dw) w=wmin+dw*dble(i-1)+dw endif 900 w=w-dw*si i=i-is if(i.lt.1.or.i.gt.npoint)goto 901 du=((w-wo)/dd)**2 if(LG.and.du.lt.20.0d0)then ss=exp(-du)/dd/spi sd(i)=sd(i)+ss*zd sr(i)=sr(i)+ss*zr goto 900 else if(.not.LG.and.du.lt.1.0d5)then ss=1.0d0/dd/(du+1.0d0)/pi sd(i)=sd(i)+ss*zd sr(i)=sr(i)+ss*zr goto 900 endif endif 901 continue return end subroutine initab(lwr,lnum,ltab,units) implicit none logical lnum,ltab integer*4 units,lwr character*4 lu(3) data lu/' GHz',' MHz',' cm-1'/ if(lwr.gt.0)then if(lnum)then open(66,file='RCDN.TAB') else open(66,file='RCD.TAB') endif write(66,662)lu(units+1) 662 format(' freq',11x,'D',11x,'R,'7x,'4R/D', 1 ' J i M Je ie Me pG pE',/, 1 'transition ',a4,5x,'debye^2',5x,'debye^2', 1 38x,'%',8x,'%') write(66,660) 660 format(90(1h-)) endif if(ltab)then if(lnum)then open(67,file='DETAILN.TAB') else open(67,file='DETAIL.TAB') endif endif return end subroutine recordtrans(nbuf,nba,J,JE,id,ide,M,ME,wo,zd,zr,zq, 1zrb,irb,debye,lwr,io,pG,pE,p,wmin,dw,dd,npoint,LG,sd,sr,spi,pi, 1units) implicit none integer*4 nbuf,nba,J,JE,id,ide,M,ME,irb(6,nbuf),lwr,npoint,io, 1jt,jet,idt,idet,mt,met,units real*8 wo,zd,zr,zq,debye,zrb(3,nbuf),cid,pG,pE,p,zdd,zqrd, 1wmin,dw,dd,sd(*),sr(*),spi,pi,gcm,wu logical LG gcm=0.0333564095198152d0 c 1GHz=0.0333564095198152 cm-1 io=io+1 zdd=zd*debye**2 zqrd=(zr+zq)*debye**2 cid=4.0d0*zqrd/zdd c temporary variables so that compiler does not complain: jt=J jet=JE idt=id idet=ide mt=M met=ME c save big transitions call saveb(nbuf,nba,jt,jet,idt,idet,mt,met,wo,zdd,zqrd,zrb,irb) c write to .TAB file, in debyes^2: if(units.eq.0)then wu=wo elseif(units.eq.1)then wu=wo*1.0d3 elseif(units.eq.2)then wu=wo*gcm else call report('unknown units requested') endif if(lwr.gt.0)write(66,661)io,wu,zdd,zqrd,cid,J,id,M, 1JE,ide,ME,pG*100.0d0,pE*100.0d0 661 format(i7,f11.4,3e12.4,6i3,2f9.5) c multiply by probability and spread over spectrum: call sprd(zdd,zqrd,p,wo,wmin,dw,dd,npoint,LG,sd,sr,spi,pi) return end subroutine partf(Jmin,Jmax,Q,ABC,Te,tf,pi,lsym) implicit none integer*4 J,Jmin,Jmax,IERR,id,Ndim real*8 Q,qa,ABC(6),Te,tf,pi,j2 logical lsym real*8,allocatable:: H(:,:),Cjk(:,:),et(:),ek(:) Q=0.0d0 do 1 J=Jmin,Jmax Ndim=2*J+1 j2=dble(Ndim) 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,ABC(1),ABC(2),ABC(3),Ndim) else call seth(h,J,ABC(1),ABC(2),ABC(3)) IERR=0 call TRED12(Ndim,H,Cjk,ek,2,IERR,et) endif do 2 id=1,Ndim 2 Q=Q+j2*exp(-ek(id)/tf/Te) 1 deallocate(H,Cjk,et,ek) qa=dsqrt(pi*(3.0d0*Te*tf)**3/(ABC(1)-ABC(2)+ABC(3))**3) write(6,600)Q,qa 600 format(' Partitition function: numerical = ',G12.4,/, 1 ' approximate = ',G12.4,/,/) return end subroutine smooth(t,a,b,c,d) implicit none real*8 t,a,b,c,d if(dabs(a).lt.t)a=0.0d0 if(dabs(b).lt.t)b=0.0d0 if(dabs(c).lt.t)c=0.0d0 if(dabs(d).lt.t)d=0.0d0 return end