program rcd c extracts from G.OUT implicit none integer*4 nf0 parameter (nf0=170) c nf0 ...maximal factorial for real*8 (for real*16 nf0=1754) integer*4 i,IERR,j,k,id,jd,Jmax,Jmin,Ndim,ic,Ndime,NPOINT, 1JE,KE,M,ide,jde,jdstart,jdend,jdestart,jdeend,mm,io,lwr, 1iu,is,JEMIN,JEMAX,f,a,b,ME,ia real*8 abc(3),g(3,3),khz,p,qp,AA,c,ta,cid,kke,BB,CC,Te,AE,BE,CE, 1DD,wo,WMAX,WMIN,WEXC,u(3),u0i,u0r,u1i,u1r,umr,umi,zd,ss,zrr,zri, 1co,ithrj,dw,factorial(nf0),tf,kcm,sk,pi,jje,rtol,gr(3,3),gi(3,3), 1skm,d10,dm0,d00,ur,ui,spi,si,myr,cjje,coselJ, 1du,w,sn,uN,debye,jkp,jkm,jkem,jkep,elr,eli,av, 1ck0r,ck0i,ck1r,ck1i,ck2r,ck2i,ck3r,ck3i,ck4r,ck4i real*8,allocatable:: H(:,:),Cjk(:,:),et(:),ek(:),He(:,:), 1Cjke(:,:),ete(:),eke(:),sd(:),sr(:),t1(:),tm(:), 1gu1r(:),gu1i(:),gumr(:),gumi(:),mXYr(:,:,:),mXYi(:,:,:), 1pXYr(:,:,:),pXYi(:,:,:),cmmer(:),cmmei(:), 1cmmert(:),cmmeit(:) LOGICAL lzmat,lg,lcm,lsym character*4 uni call initfac(factorial,nf0) pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) call setpar(c,kcm,khz,lzmat,Jmin,Jmax,LG,lwr,AA,BB,CC,dw, 1AE,BE,CE,Te,WEXC,WMAX,WMIN,DD,NPOINT,u,lcm,tf,ta,uni, 1rtol,JEMIN,JEMAX,uN,debye) call lookatout(lzmat,g,abc,ta,AA,BB,CC,AE,BE,CE,u,debye) c transforms cartesian g-tensor to spherical G-tensor: call makeg(g,gr,gi) write(6,*)'Spherical tensor real:' call wrm(gr,3) write(6,*)'Spherical tensor imaginary:' call wrm(gi,3) c write some selected strengths: call wrinv(g,u,uN,debye) lsym=dabs(AA-BB).lt.rtol write(6,60000)AA,BB,CC,uni,AE,BE,CE,uni, 1JMIN,JMAX,WEXC,WMIN,WMAX,uni, 1NPOINT,LG,Te,LWR,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 ' Jmin Jmax :',2i12,/, 1 ' WEXC WMIN WMAX:',3F12.6,a4,/, 1 ' Npoints LG :',i12,l12,/, 1 ' Temperature :',f12.6,' K',/, 1 ' Large output :',i12,/, 1 ' ux uy uz :',3f12.6,/, 1 ' Symmetric top :',l12,' (tol:',g12.6,')') c spherical tensor molecular dipoles: c u1=-(ux+iuy)/sqrt(2): u1r=-u(1)/dsqrt(2.00d0) u1i=-u(2)/dsqrt(2.00d0) c u0=uz: u0r= u(3) u0i=0.0d0 c u1= (ux-iuy)/sqrt(2): umr= u(1)/dsqrt(2.00d0) umi=-u(2)/dsqrt(2.00d0) allocate(sd(NPOINT),sr(NPOINT)) do 41 i=1,NPOINT sr(i)=0.0d0 41 sd(i)=0.0d0 qp=0.0d0 io=0 c c ic: first cycle -build Qp (partition function) c second cycle - estimate probabilities if(lwr.gt.0)then open(66,file='RCD.TAB') write(66,662)uni 662 format(' freq D R 4R/D J j Je', 1 ' k Probab',/, 1 'transition ',a4,' debye^2 10^-4debye^2') write(66,660) 660 format(60(1h-)) endif do 100 ic=1,2 do 100 J=Jmin,Jmax Ndim=2*J+1 allocate(H(Ndim,Ndim),Cjk(Ndim,Ndim),et(Ndim),ek(Ndim), 1t1(Ndim),tm(Ndim),cmmer(2*2*3*Ndim),cmmei(2*2*3*Ndim), 1cmmert(2*2*3*Ndim),cmmeit(2*2*3*Ndim), 1gu1r(Ndim),gu1i(Ndim),gumr(Ndim),gumi(ndiM), 1mXyr(2,2,Ndim),mXYi(2,2,Ndim),pXYr(2,2,Ndim),pXYi(2,2,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 if(ic.eq.1)then c estimate partition function only, 2J+1... magnetic degeneracy: do 1034 id=1,Ndim 1034 qp=qp+exp(-ek(id)/tf/Te)*dble(Ndim) else if(lwr.gt.1)then write(6,6002)J 6002 format(' J = ',I4,', energies:') call wre(ek,Ndim,uni,Te,tf,qp,0.0d0) call wrm(H,Ndim) endif do 1001 JE=max(0,J-1),J+1 Ndime=2*JE+1 jje=dsqrt(dble(2*JE+1)*dble(2*J+1)) cjje=coselJ(J,JE) 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) if(lwr.gt.1)then write(6,60021)JE 60021 format(' JE = ',I4,' energies:') call wre(eke,Ndime,uni,tf,Te,qp,WEXC) call wrm(He,Ndime) endif endif c pre-calculate M-factors: do 1035 M=-J,J t1(M+J+1)=ithrj(JE,1,J,M-1, 1,-M,factorial)*skm(M) 1035 tm(M+J+1)=ithrj(JE,1,J,M+1,-1,-M,factorial)*skm(M) ia=0 do 10351 f=1,3 do 10351 a=1,2 M=-J-1 do 10351 mm=1,Ndim M=M+1 ia=ia+1 call coselM(J,M,a,JE,M-1,elr,eli) cmmer(ia)=elr cmmei(ia)=eli call coselM(JE,M-1,a,J,M,elr,eli) cmmert(ia)=elr cmmeit(ia)=eli ia=ia+1 call coselM(J,M,a,JE,M+1,elr,eli) cmmer(ia)=elr cmmei(ia)=eli call coselM(JE,M+1,a,J,M,elr,eli) cmmert(ia)=elr 10351 cmmeit(ia)=eli c ground state id = Cjk(jd,id)*|JKM>(jd): do 103 id=1,Ndim p=exp(-ek(id)/tf/Te)/qp*dble(Ndim) c excited state ide = Cjke(jde,ide)*|JEKEME>(jde): do 103 ide=1,Ndime c transition energy |G> -> |E>: wo=WEXC+eke(ide)-ek(id) c transition moments, c = Cjke_jde Cjk_jd , etc if(wo.gt.0.0d0)then do 1036 M=-J,J c real and imaginary: gu1r(M+J+1)=0.0d0 gu1i(M+J+1)=0.0d0 gumr(M+J+1)=0.0d0 gumi(M+J+1)=0.0d0 c , , do 1036 a=1,2 do 1036 f=1,2 pXYr(a,f,M+J+1)=0.0d0 pXYi(a,f,M+J+1)=0.0d0 c : mXYr(a,f,M+J+1)=0.0d0 1036 mXYi(a,f,M+J+1)=0.0d0 if(lsym)then jdstart=id jdend=id jdestart=ide jdeend=ide else jdstart=1 jdend=Ndim jdestart=1 jdeend=Ndime endif do 1032 jd=jdstart,jdend K=-J+jd-1 jkp=dsqrt(dble(J*(J+1)-K*(K-1))) jkm=dsqrt(dble(J*(J+1)-K*(K+1))) sk=jje*skm(K) do 1032 jde=jdestart,jdeend KE=-JE+jde-1 if(abs(K-KE).le.2)then jkep=dsqrt(dble(JE*(JE+1)-KE*(KE-1))) jkem=dsqrt(dble(JE*(JE+1)-KE*(KE+1))) kke=dble(KE+K) co=Cjke(jde,ide)*Cjk(jd,id) c 1)make electric dipoles: if(abs(K-KE).le.1)then c (-1)^(K)sqrt((2J+1)(2JE+1)) (JE1J//KEq-K)(JE1J//MEq-M)um1q: if(KE.eq.K)THEN d00=ithrj(JE,1,J,KE , 0,-K,factorial) ur=u0r*d00*sk*co ui=u0i*d00*sk*co else if(KE.eq.K+1)THEN dm0=ithrj(JE,1,J,KE ,-1,-K,factorial) ur=umr*dm0*sk*co ui=umi*dm0*sk*co else d10=ithrj(JE,1,J,KE , 1,-K,factorial) ur=u1r*d10*sk*co ui=u1i*d10*sk*co endif endif M=-J-1 do 11372 mm=1,Ndim M=M+1 c =CGJ CEJE etc.: c ME=M-1: gu1r(mm)=gu1r(mm)+ur*t1(mm) gu1i(mm)=gu1i(mm)+ui*t1(mm) c ME=M+1: gumr(mm)=gumr(mm)+ur*tm(mm) 11372 gumi(mm)=gumi(mm)+ui*tm(mm) ia=0 do 10372 f=1,3 do 10372 a=1,2 call coselK(J,K,f,JE,KE,ck0r,ck0i) M=-J-1 do 10372 mm=1,Ndim M=M+1 ia=ia+1 c again by cosinus directional elements Xf, elr=ck0r*cmmer(ia)-ck0i*cmmei(ia) eli=ck0i*cmmer(ia)+ck0r*cmmei(ia) pXYr(1,a,mm)=pXYr(1,a,mm)+u(f)*elr*co*cjje pXYi(1,a,mm)=pXYi(1,a,mm)+u(f)*eli*co*cjje ia=ia+1 elr=ck0r*cmmer(ia)-ck0i*cmmei(ia) eli=ck0i*cmmer(ia)+ck0r*cmmei(ia) pXYr(2,a,mm)=pXYr(2,a,mm)+u(f)*elr*co*cjje 10372 pXYi(2,a,mm)=pXYi(2,a,mm)+u(f)*eli*co*cjje ENDIF c 2)make magnetic dipoles ia=0 do 10371 f=1,3 do 10371 a=1,2 call coselK(JE,KE,f,J,K,ck0r,ck0i) call coselK(JE,KE,f,J,K+1,ck1r,ck1i) call coselK(JE,KE-1,f,J,K,ck2r,ck2i) call coselK(JE,KE,f,J,K-1,ck3r,ck3i) call coselK(JE,KE+1,f,J,K,ck4r,ck4i) M=-J-1 do 10371 mm=1,Ndim M=M+1 c Saltzman facon: c call magelx(J,K,M,g,JE,KE,M+1,elr,eli) c if(elr**2+eli**2.lt.1.00d-10) c 1 call magelx(JE,KE,M+1,g,J,K,M,elr,eli) c sXr(mm)=sXr(mm)+elr c sXi(mm)=sXi(mm)+eli c by cosinus directional elements Xf: c both MX, My, all M': b=0 do 10371 ME=M-1,M+1,2 b=b+1 ia=ia+1 elr=(ck1r*cmmert(ia)-ck1i*cmmeit(ia))*co*cjje*jkp eli=(ck1i*cmmert(ia)+ck1r*cmmeit(ia))*co*cjje*jkp mXYr(b,a,mm)=mXYr(b,a,mm)+gr(f,1)*elr-gi(f,1)*eli mXYi(b,a,mm)=mXYi(b,a,mm)+gr(f,1)*eli+gi(f,1)*elr elr=(ck2r*cmmert(ia)-ck2i*cmmeit(ia))*co*cjje*jkem eli=(ck2i*cmmert(ia)+ck2r*cmmeit(ia))*co*cjje*jkem mXYr(b,a,mm)=mXYr(b,a,mm)+gr(f,1)*elr-gi(f,1)*eli mXYi(b,a,mm)=mXYi(b,a,mm)+gr(f,1)*eli+gi(f,1)*elr elr=(ck3r*cmmert(ia)-ck3i*cmmeit(ia))*co*cjje*jkm eli=(ck3i*cmmert(ia)+ck3r*cmmeit(ia))*co*cjje*jkm mXYr(b,a,mm)=mXYr(b,a,mm)+gr(f,2)*elr-gi(f,2)*eli mXYi(b,a,mm)=mXYi(b,a,mm)+gr(f,2)*eli+gi(f,2)*elr elr=(ck4r*cmmert(ia)-ck4i*cmmeit(ia))*co*cjje*jkep eli=(ck4i*cmmert(ia)+ck4r*cmmeit(ia))*co*cjje*jkep mXYr(b,a,mm)=mXYr(b,a,mm)+gr(f,2)*elr-gi(f,2)*eli mXYi(b,a,mm)=mXYi(b,a,mm)+gr(f,2)*eli+gi(f,2)*elr elr=(ck0r*cmmert(ia)-ck0i*cmmeit(ia))*co*cjje*kke eli=(ck0i*cmmert(ia)+ck0r*cmmeit(ia))*co*cjje*kke mXYr(b,a,mm)=mXYr(b,a,mm)+gr(f,3)*elr-gi(f,3)*eli 10371 mXYi(b,a,mm)=mXYi(b,a,mm)+gr(f,3)*eli+gi(f,3)*elr endif 1032 continue c c m-average rotational strengths: zd=0.0d0 zrr=0.0d0 zri=0.0d0 myr=0.0d0 M=-J-1 do 1038 mm=1,Ndim zd=zd+gu1r(mm)*gu1r(mm)+gu1i(mm)*gu1i(mm) 1 + gumr(mm)*gumr(mm)+gumi(mm)*gumi(mm) c M->M+1,M->M-1,X and Y: myr=myr+pXYi(1,1,mm)*mXYr(1,1,mm)+pXYr(1,1,mm)*mXYi(1,1,mm) 1 +pXYi(2,1,mm)*mXYr(2,1,mm)+pXYr(2,1,mm)*mXYi(2,1,mm) 1 +pXYi(1,2,mm)*mXYr(1,2,mm)+pXYr(1,2,mm)*mXYi(1,2,mm) 1 +pXYi(2,2,mm)*mXYr(2,2,mm)+pXYr(2,2,mm)*mXYi(2,2,mm) c M->M+1,X : zrr=zrr+pXYr(2,1,mm)*mXYr(2,1,mm)-pXYi(2,1,mm)*mXYi(2,1,mm) 1038 zri=zri+pXYi(2,1,mm)*mXYr(2,1,mm)+pXYr(2,1,mm)*mXYi(2,1,mm) io=io+1 c from au back to debyes,r with prefactor 1e4: zd=zd*debye**2 myr=myr*debye**2*uN*1.0d4 zrr=zrr*debye**2*uN*1.0d4 zri=zri*debye**2*uN*1.0d4 cid=4.0d0*zri/(zd+1.0d-10) if(cid.gt.9999.9d0)cid=9999.9d0 if(cid.lt.-9999.9d0)cid=-9999.9d0 c average norm for M-averaging: av=6.0d0/dble(2*J+1) zd=av*zd/2.0d0 myr=myr*av/2.0d0 zri=av*zri zrr=av*zrr if(lwr.gt.0) 1 write(66,661)io,wo,zd,zri,myr,cid,J,id,JE,ide,p*100.0d0 661 format(i7,f11.4,4f12.5,4i3,f10.5) c multiply by probability and spread over spectrum: zd=p*zd zri=p*zri 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*zri 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*zri goto 900 endif endif 901 continue endif 103 continue 1001 deallocate(He,Cjke,ete,eke) endif 100 deallocate(H,Cjk,et,ek,cmmer,cmmei,cmmert,cmmeit, 1t1,tm,gu1r,gu1i,gumr,gumi,mXYr,mXYi,pXYr,pXYi) if(lwr.gt.0)then write(66,660) close(660) write(6,*)'RCD.TAB' endif sn=0.0d0 do 902 i=1,npoint 902 sn=sn+sd(i)**2 sn=dsqrt(sn) do 903 i=1,npoint sd(i)=sd(i)/sn 903 sr(i)=sr(i)/sn call wrspec2('D.PRN',npoint,wmin,dw,sd) call wrspec2('R.PRN',npoint,wmin,dw,sr) deallocate(sd,sr) STOP end c ================================================================== subroutine wrspec2(sn,nw,wmin,dw,spec) implicit none character*(*) sn integer*4 nw,iw real*8 wmin,dw,spec(*),w open(8,file=sn) do 1 iw=1,nw w=wmin+dble(iw-1)*dw 1 write(8,805)w,spec(iw) 805 format(f18.8,' ',g25.12) 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,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 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 Function rthrj(x1,x2,x3,y1,y2,y3) c {real three - J symbols} implicit none real*8 rthrj,x1,x2,x3,y1,y2,y3,ar,u,sum,msign,fac 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) 1.or.(m1d+m2d.ne.-m3d).or.(j1d+j2d.lt.j3d).or.(j2d+j3d.lt.j1d) 2.or.(j3d+j1d.lt.j2d).or. (mod(j1d-m1d,2).ne.0) 3.or.(mod(j2d-m2d,2).ne.0).or.(mod(j3d-m3d,2).ne.0) 4.or.(mod(j1d+j2d+j3d,2).ne.0) )then rthrj=0.0d0 else if (mod(j1d-j2d-m3d,2).eq.0)then msign=1.0d0 else msign=-1.0d0 endif s=nint(x1+x2+x3) ar=fac(s-j3d)*fac(s-j2d)*fac(s-j1d)/fac(s+1) ar=ar*fac((j1d+m1d) / 2)*fac((j1d-m1d) / 2) ar=ar*fac((j2d+m2d) / 2)*fac((j2d-m2d) / 2) ar=ar*fac((j3d+m3d) / 2)*fac((j3d-m3d) / 2) ar=msign*sqrt(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. 2 (i5.ge.0)).and.(i6.ge.0)) then u=fac(i1)*fac(i2) u=u*fac(i3)*fac(i4) u=u*fac(i5)*fac(i6) sum=sum+msign/u endif 1 continue ar=ar*sum if(abs(ar).lt.1.0d-38)then rthrj=0.0d0 else rthrj=ar endif endif return end Function ithrj(x1,x2,x3,y1,y2,y3,fc) c {integer three - J symbols} implicit none real*8 ar,sum,ithrj,fc(*) integer*4 ni,i1,i2,i3,i4,i5,i6, 1x1,x2,x3,y1,y2,y3 if((abs(y1).gt.x1).or.(abs(y2).gt.x2).or.(abs(y3).gt.x3) 1.or.(y1+y2.ne.-y3) 1.or.(x1+x2.lt.x3).or.(x2+x3.lt.x1).or.(x3+x1.lt.x2))then ithrj=0.0d0 else ar=sqrt(fc(1+x1+x2-x3)*fc(1+x1+x3-x2)*fc(1+x2+x3-x1) 1 *fc(1+x1+y1)*fc(1+x1-y1)*fc(1+x2+y2)*fc(1+x2-y2) 2 *fc(1+x3+y3)*fc(1+x3-y3)/fc(1+x1+x2+x3+1)) sum=0.0d0 do 1 ni=0,x1+x2+x3 i1=ni i2=x1+x2-x3-ni i3=x1-y1-ni i4=x2+y2-ni i5=x3-x2+y1+ni i6=x3-x1-y2+ni if(mod(ni,2).ne.0) then if(i2.ge.0.and.i3.ge.0.and.i4.ge.0.and.i5.ge.0.and.i6.ge.0) 1 sum=sum 1 -1.0d0/(fc(1+i1)*fc(1+i2)*fc(1+i3)*fc(1+i4)*fc(1+i5)*fc(1+i6)) else if(i2.ge.0.and.i3.ge.0.and.i4.ge.0.and.i5.ge.0.and.i6.ge.0) 1 sum=sum 1 +1.0d0/(fc(1+i1)*fc(1+i2)*fc(1+i3)*fc(1+i4)*fc(1+i5)*fc(1+i6)) endif 1 continue ithrj=ar*sum endif 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= sqrt(2.0d0*j3+1.0d0)*rthrj(j1,j2,j3,m1,m2,-m3) else rcg=-sqrt(2.0d0*j3+1.0d0)*rthrj(j1,j2,j3,m1,m2,-m3) endif 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 Function Fac(n) implicit none integer*4 ii,n real*8 Fac,bb,aii bb=1.0d0 aii=1.0d0 do 1 ii=2,n aii=aii+1.0d0 1 bb=bb*aii fac=bb return end subroutine initfac(factorial,n) implicit none integer*4 i,n real*8 fac,factorial(*) do 1 i=1,n 1 factorial(i)=fac(i-1) 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(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 subroutine report(s) character*(*) s write(6,*)s stop end subroutine setpar(c,kcm,khz,lzmat,Jmin,Jmax,LG,lwr,AA,BB,CC,dw, 1AE,BE,CE,Te,WEXC,WMAX,WMIN,DD,NPOINT,u,lcm,tf,ta,uni, 1rtol,JEMIN,JEMAX,uN,debye) implicit none real*8 c,kcm,khz,AA,BB,CC,AE,BE,CE,Te,WEXC,WMAX,WMIN,u(3), 1DD,tf,ta,dw,rtol,uN,debye integer*4 Jmin,Jmax,NPOINT,lwr,JEMIN,JEMAX,i logical lzmat,LG,lcm,lex character*4 uni character*80 FN 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: khz=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: AA=0.0d0 BB=0.0d0 CC=0.0d0 AE=0.0d0 BE=0.0d0 CE=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 lcm=.false. 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,*)AA if(FN(1:2).eq.'BB')read(55,*)BB if(FN(1:2).eq.'CC')read(55,*)CC if(FN(1:2).eq.'AE')read(55,*)AE if(FN(1:2).eq.'BE')read(55,*)BE if(FN(1:2).eq.'CE')read(55,*)CE if(FN(1:2).eq.'DD')read(55,*)DD if(FN(1:4).eq.'WMAX')read(55,*)WMAX 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.'WMIN')read(55,*)WMIN if(FN(1:4).eq.'NPOI')read(55,*)NPOINT if(FN(1:4).eq.'GAUS')read(55,*)LG if(FN(1:4).eq.'TEMP')read(55,*)Te if(FN(1:3).eq.'LWR')read(55,*)lwr if(FN(1:3).eq.'LCM')read(55,*)lcm if(FN(1:3).eq.'LZM')read(55,*)lzmat if(FN(1:4).eq.'WEXC')read(55,*)WEXC if(FN(1:4).eq.'UXYZ')read(55,*)(u(i),i=1,3) if(FN(1:4).eq.'RTOL')read(55,*)rtol goto 551 555 close(55) endif if(lcm)then uni='CM-1' tf=kcm ta=1.0d9/c else uni=' GHz' tf=khz ta=1.0d0 endif dw=(WMAX-WMIN)/dble(NPOINT) return end subroutine makeg(g,gr,gi) c transforms cartesian g-tensor to spherical G-tensor implicit none real*8 g(3,3),gr(3,3),gi(3,3),s2 integer*4 f s2=2.0d0 c real and imaginary: do 1 f=1,3 gr(f,1)= g(f,1) /s2 gi(f,1)= -g(f,2)/s2 gr(f,2)= g(f,1) /s2 gi(f,2)= g(f,2)/s2 gr(f,3)= g(f,3) 1 gi(f,3)= 0.0d0 return end function skm(K) c skm=(-1)**K implicit none real*8 skm integer*4 K if(mod(K,2).eq.0)then skm=1.0d0 else skm=-1.0d0 endif 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,uni,tf,Te,qp,WEXC) implicit none real*8 e(*),tf,Te,qp,WEXC,p integer*4 i,n character*4 uni do 1 i=1,n p=exp(-(e(i)+WEXC)/tf/Te)/qp*dble(n) 1 write(6,600)i,e(i),uni,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 function coselJ(J,JE) c cosinus directional element = el + i eli c = = jj jk jm c J=part only implicit none integer*4 J,JE real*8 coselJ if(JE.eq.J+1)then coselJ=1.0d0/(dble(4*(J+1))*dsqrt(dble((2*J+1)*(2*J+3)))) else if(JE.eq.J-1)then coselJ=1.0d0/(dble(4*J)*dsqrt(dble(4*J**2-1))) else if(JE.eq.J)then coselJ=1.0d0/dble(4*J*(J+1)) else coselJ=0.0d0 endif endif endif return end subroutine cosel(J,K,M,G,f,JE,KE,ME,elr,eli) c cosinus directional element = el + i eli c = = jj jk jm implicit none integer*4 J,K,M,G,F,JE,KE,ME,s real*8 elr,eli,jj,jky,jmy,jkr,jki,jmr,jmi elr=0.0d0 eli=0.0d0 if(JE.eq.J+1)then jj=1.0d0/(dble(4*(J+1))*dsqrt(dble((2*J+1)*(2*J+3)))) else if(JE.eq.J-1)then jj=1.0d0/(dble(4*J)*dsqrt(dble(4*J**2-1))) else if(JE.eq.J)then jj=1.0d0/dble(4*J*(J+1)) else return endif endif endif if(f.eq.3)then c : if(K.eq.KE)then if(JE.eq.J+1)then jkr=2.0d0*dsqrt(dble((J+1)**2-K**2)) jki=0.0d0 else if(JE.eq.J-1)then jkr=-2.0d0*dsqrt(dble(J**2-K**2)) jki=0.0d0 else if(JE.eq.J)then jkr=dble(K+K) jki=0.0d0 else return endif endif endif else return endif else c and : if(KE.eq.K+1)then s=1 else if(KE.eq.K-1)then s=-1 else return endif endif if(JE.eq.J+1)then jky=-dble(s)*dsqrt(dble((J+s*K+1)*(J+s*K+2))) else if(JE.eq.J-1)then jky=-dble(s)*dsqrt(dble((J-s*K)*(J-s*K-1))) else if(JE.eq.J)then jky=dsqrt(dble(J*(J+1)-K*(K+s))) else return endif endif endif if(f.eq.2)then jkr=jky jki=0.0d0 else jkr=0.0d0 jki=jky*dble(s) endif endif if(G.eq.3)then c : if(M.eq.ME)then if(JE.eq.J+1)then jmr=2.0d0*dsqrt(dble((J+1)**2-M**2)) jmi=0.0d0 else if(JE.eq.J-1)then jmr=-2.0d0*dsqrt(dble(J**2-M**2)) jmi=0.0d0 else if(JE.eq.J)then jmr=dble(M+M) jmi=0.0d0 else return endif endif endif else return endif else c and : if(ME.eq.M+1)then s=1 else if(ME.eq.M-1)then s=-1 else return endif endif if(JE.eq.J+1)then jmy=-dble(s)*dsqrt(dble((J+s*M+1)*(J+s*M+2))) else if(JE.eq.J-1)then jmy=-dble(s)*dsqrt(dble((J-s*M)*(J-s*M-1))) else if(JE.eq.J)then jmy=dsqrt(dble(J*(J+1)-M*(M+s))) else return endif endif endif if(G.eq.2)then jmr=jmy jmi=0.0d0 else jmr=0.0d0 jmi=-jmy*dble(s) endif endif elr=jj*(jkr*jmr-jki*jmi) eli=jj*(jki*jmr+jkr*jmi) 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 subroutine magelx(J,K,M,G,JE,KE,ME,elr,eli) c magnetic element = el + i eli c = = jj jk jm implicit none integer*4 J,K,M,JE,KE,ME,s real*8 elr,eli,jj,jkr,jki,jmr,jmi,g(3,3),f,f1,f2 elr=0.0d0 eli=0.0d0 c : if(JE.eq.J+1)then jj=1.0d0/(dble(4*JE)*dsqrt(dble(4*JE**2-1))) else if(JE.eq.J-1)then jj=1.0d0/(dble(4*J)*dsqrt(dble(4* J**2-1))) else if(JE.eq.J)then jj=1.0d0/dble(4*J*(J+1)) else return endif endif endif c : if(ME.eq.M+1)then s=1 else if(ME.eq.M-1)then s=-1 else return endif endif if(JE.eq.J+1)then jmr=0.0d0 jmi=dsqrt(dble((J+s*M+1)*(J+s*M+2))) else if(JE.eq.J-1)then jmr=0.0d0 jmi=dsqrt(dble((JE-s*ME+1)*(JE-s*ME+2))) else if(JE.eq.J)then jmr=0.0d0 jmi=-dble(s)*dsqrt(dble(J*(J+1)-M*ME)) else return endif endif endif c : s=KE-K if(J.eq.JE)then if(abs(s).le.2)then if(abs(s).eq.2)then f=dsqrt(dble((J*(J+1)-K*(K+s))*(J*(J+1)-(K-s)*(K+s+s)))) jkr=-f*(g(1,1)-g(2,2))*0.5d0 jki= f*(g(1,2)+g(2,1))*0.5d0*dble(s) endif if(abs(s).eq.1)then f=dsqrt(dble((J*(J+1)-K*(K+s))))*dble(K+K+s) jkr=-f*(g(3,2)+g(2,3))*0.5d0 jki= f*(g(3,1)+g(1,3))*0.5d0*dble(s) endif if(abs(s).eq.0)then jkr=dble(J*(J+1)-K**2)*(g(1,1)+g(2,2))+2.0d0*dble(K**2)*g(3,3) jki=0.0d0 endif else return endif else if(J+1.eq.JE)then if(abs(s).le.2)then if(abs(s).eq.2)then f=dsqrt(dble((J*(J+1)-K*(K+s))*((J+s*K+2)*(K+s*K+3)))) 1 +dsqrt(dble(((J+1)*(J+2)-(K+s)*(K+s+s))*((J+s*K+1) 1 *(K+s*K+2)))) jkr= 0.25d0*f*(g(1,1)-g(2,2)) jki=-0.25d0*f*(g(1,2)-g(2,1)) endif if(abs(s).eq.1)then f1=dsqrt(dble((J*(J+1)-K*(K+s))*((J+1)**2-(K+s)**2))) 1 +dsqrt(dble(((J+1)*(J+2)-K*(K+s))*((J+1)**2-K**2))) f2=dble(2*K+s)*dsqrt(dble((J+s*K+1)*(J+s*K+2))) jkr= 0.5d0*f1*g(3,2)-0.5d0*f2*g(2,3)*dble(s) jki= 0.5d0*f1*g(3,1)-0.5d0*f2*g(1,3) endif if(abs(s).eq.0)then f1=-dsqrt(dble((J*(J+1)-K*(K-1))*(J+K)*(J+K+1))) 1 +dsqrt(dble((J*(J+1)-K*(K+1))*(J-K)*(J-K+1))) 1 +dsqrt(dble(((J+1)*(J+2)-K*(K-1))*(J-K+1)*(J-K+2))) 1 -dsqrt(dble(((J+1)*(J+2)-K*(K+1))*(J+K+1)*(J+K+2))) f2=dble(2*K)*dsqrt(dble((J+1)**2-K**2)) jkr=0.25d0*f1*(g(1,1)+g(2,2))-f2*g(3,3) jki=0.25d0*f1*(g(2,1)-g(1,2)) endif else return endif else return endif endif elr=jj*(jkr*jmr-jki*jmi) eli=jj*(jki*jmr+jkr*jmi) return end subroutine coselK(J,K,f,JE,KE,jkr,jki) c cosinus directional element = el + i eli c = = jj jk jm implicit none integer*4 J,K,F,JE,KE,s real*8 jky,jkr,jki jkr=0.0d0 jki=0.0d0 if(f.eq.3)then c : if(K.eq.KE)then if(JE.eq.J+1)then jkr=2.0d0*dsqrt(dble((J+1)**2-K**2)) jki=0.0d0 else if(JE.eq.J-1)then jkr=-2.0d0*dsqrt(dble(J**2-K**2)) jki=0.0d0 else if(JE.eq.J)then jkr=dble(K+K) jki=0.0d0 else return endif endif endif else return endif else c and : if(KE.eq.K+1)then s=1 else if(KE.eq.K-1)then s=-1 else return endif endif if(JE.eq.J+1)then jky=-dble(s)*dsqrt(dble((J+s*K+1)*(J+s*K+2))) else if(JE.eq.J-1)then jky=-dble(s)*dsqrt(dble((J-s*K)*(J-s*K-1))) else if(JE.eq.J)then jky=dsqrt(dble(J*(J+1)-K*(K+s))) else return endif endif endif if(f.eq.2)then jkr=jky jki=0.0d0 else jkr=0.0d0 jki=jky*dble(s) endif endif return end subroutine coselM(J,M,G,JE,ME,jmr,jmi) c cosinus directional element = el + i eli c = = jj jk jm implicit none integer*4 J,M,G,JE,ME,s real*8 jmy,jmr,jmi jmr=0.0d0 jmi=0.0d0 if(G.eq.3)then c : if(M.eq.ME)then if(JE.eq.J+1)then jmr=2.0d0*dsqrt(dble((J+1)**2-M**2)) jmi=0.0d0 else if(JE.eq.J-1)then jmr=-2.0d0*dsqrt(dble(J**2-M**2)) jmi=0.0d0 else if(JE.eq.J)then jmr=dble(M+M) jmi=0.0d0 else return endif endif endif else return endif else c and : if(ME.eq.M+1)then s=1 else if(ME.eq.M-1)then s=-1 else return endif endif if(JE.eq.J+1)then jmy=-dble(s)*dsqrt(dble((J+s*M+1)*(J+s*M+2))) else if(JE.eq.J-1)then jmy=-dble(s)*dsqrt(dble((J-s*M)*(J-s*M-1))) else if(JE.eq.J)then jmy=dsqrt(dble(J*(J+1)-M*(M+s))) else return endif endif endif if(G.eq.2)then jmr=jmy jmi=0.0d0 else jmr=0.0d0 jmi=-jmy*dble(s) endif endif return end