PROGRAM CIDR IMPLICIT none integer*4 k,nat real*8 w,A(27),AL(9),G(9),CM,ecm,P1,doc,EXCA,YDY,YDX, 1roa1,ram1,CID1,qq(3,3),am(3,3),alk(9),gk(9),ak(27), 1roav,ramv,dov real*8,allocatable:: r(:) integer*4,allocatable:: q(:) call readpol(AL,'POL.TTT',w) call wr('ALpha',AL) call readpol(G,'GP.TTT',w) call wr('G',G) call readpol(A,'A.TTT',w) call wra('A',A) c CM=219474.63d0 ecm=w*CM EXCA=1.0d7/ecm*10.0d0 write(6,601)EXCA 601 format(/,' Excitation ',F10.0,' A',/) call dr(w,AL,G,A,YDX,YDY,ram1,P1,roa1,CID1,doc) write(6,602)' ICP/SCP 180:',YDX,YDY,ram1,P1,roa1,CID1,doc 602 format(A13,/,/,'Rayleigh X Y TOT: ',3g12.4,' P:',f10.4,/, 1 ' ROA:',g12.4, 1' CID:',g12.4,' DOC:',g12.4) c rotational contributions: allocate(q(1),r(1)) call readx(0,q,r,nat) deallocate(q,r) allocate(q(nat),r(3*nat)) call readx(1,q,r,nat) c do moments of inertia: call qqd(qq,q,r,nat,am) c transform polarizabilites to them call tr(am,AL,A,G) ramv=0.0d0 roav=0.0d0 dov=0.0d0 do 1 k=1,3 write(6,600)k 600 format(' Molecular ',i1,'-axis:') c rotational tensor derivatives for each rotation: call rtd(k,alk,ak,gk,AL,A,G) call dr(w,alk,gk,ak,YDX,YDY,ram1,P1,roa1,CID1,doc) ramv=ramv+ram1/3.0d0 roav=roav+roa1/3.0d0 dov=dov+doc/3.0d0 1 write(6,602)' ICP/SCP 180:',YDX,YDY,ram1,P1,roa1,CID1,doc write(6,603)ramv,roav,roav/ramv,dov 603 format(' Average:',/, 1' Ram: ',g12.4,' Roa:',g12.4,' CID:',g12.4,' DOC:',f10.4) END subroutine tr(m,AL,A,G) IMPLICIT none integer*4 i,j,k,ii,jj,kk,jm,ijk real*8 m(3,3),AL(9),A(27),G(9),alt(9),gt(9),at(27) alt=0.0d0 gt=0.0d0 at=0.0d0 do 1 i=1,3 do 1 j=1,3 do 1 ii=1,3 do 1 jj=1,3 jm=3*(jj-1) alt(i+3*(j-1))=alt(i+3*(j-1))+m(ii,i)*m(jj,j)*AL(ii+jm) gt( i+3*(j-1))=gt( i+3*(j-1))+m(ii,i)*m(jj,j)*G( ii+jm) do 1 k=1,3 ijk=i+3*(j-1)+9*(k-1) do 1 kk=1,3 1 a(ijk)=a(ijk)+m(ii,i)*m(jj,j)*m(kk,k)*A( ii+jm+9*(kk-1)) AL=alt A=at G=gt return end subroutine setij(i,j,k) IMPLICIT none integer*4 i,j,k i=k+1 if(i.gt.3)i=1 j=i+1 if(j.gt.3)j=1 return end subroutine qqd(qq,q,r,nat,a) IMPLICIT none integer*4 nat,q(*),i,j,ia,ix,me,IERR real*8 qq(3,3),r(*),cm(3),m,ri,rj,e(3),a(3,3) parameter (me=89) real*4 mass(me) data mass/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/ m=0.0d0 do 1 ia=1,nat 1 m=m+dble(mass(q(ia))) c mass center cm=0.0d0 do 2 ix=1,3 do 2 ia=1,nat ri=r(ix+3*(ia-1)) 2 cm(ix)=cm(ix)+dble(mass(q(ia)))*ri/m c moment of inertia qq=0.0d0 do 3 i=1,3 do 3 j=1,3 do 3 ia=1,nat ri=r(i+3*(ia-1))-cm(i) rj=r(j+3*(ia-1))-cm(j) 3 qq(i,j)=qq(i,j)+dble(mass(q(ia)))*ri*rj c principal axes IERR=0 CALL TRED12(3,qq,A,e,2,IERR) write(6,600)m,cm,qq,e,A 600 format(' Mw = ',f12.3,' g/mol',/, 1' Center of mass:',3F10.2,' A',/, 1' Intertia tensor:',/, 13(3f12.2,/),/, 1' Eigenvalues:',3f12.3,/, 1' Eigenvectors:',/, 13(3f12.2,/),/) return end c =============================== subroutine readx(ic,q,r,nat) IMPLICIT none integer*4 ic,ia,ix,nat,q(*) real*8 r(*) open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat if(ic.eq.1)then do 1 ia=1,nat 1 read(9,*)q(ia),(r(ix+3*(ia-1)),ix=1,3) endif close(9) return end c =============================== SUBROUTINE readpol(A,ty,fr) IMPLICIT none real*8 A(*),fr integer*4 IX,IY,IZ character*(*) ty OPEN(90,FILE=ty,status='old') if(ty.eq.'A.TTT')then read(90,701)fr 701 format(42x,f11.6) read(90,*) do 1 IZ=1,3 read(90,*)((A(IY+(IX-1)*3+(IZ-1)*9),IY=1,IX),IX=1,3) do 1 IX=1,3 do 1 IY=1,IX-1 1 A(IX+(IY-1)*3+(IZ-1)*9)=A(IY+(IX-1)*3+(IZ-1)*9) else if(ty.eq.'POL.TTT')then read(90,702)fr 702 format(21x,f11.6) read(90,*) read(90,*)((A(IY+(IX-1)*3),IY=1,IX),IX=1,3) do 2 IX=1,3 do 2 IY=1,IX-1 2 A(IX+(IY-1)*3)=A(IY+(IX-1)*3) else if(ty.eq.'GP.TTT')then read(90,704)fr 704 format(15x,f11.6) read(90,*) read(90,*)((A(IY+(IX-1)*3),IY=1,3),IX=1,3) c inner index (IY) magnetic, last later c (does not matter for invariants) else write(90,705) 705 format(' Unknown polarizability') stop endif endif endif close(90) return end subroutine wr(s,a) implicit none real*8 a(9) character*(*) s write(6,*)s write(6,600)a 600 format(3(10x,3f10.3,/)) return end subroutine wra(s,a) implicit none real*8 a(*) integer*4 i,j,k character*(*) s write(6,*)s do 1 i=1,3 write(6,600)i 600 format(' Dipole ',i1,':') do 1 j=1,3 1 write(6,602)(a(9*(i-1)+3*(j-1)+k),k=1,3) 602 format(10x,3f10.3) return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT integer*4(I-N) DIMENSION A(N,N),Z(N,N),D(N) real*8,allocatable::E(:) allocate(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) .LE. 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 dr(w,AL,G,A,YDX,YDY,ram1,P1,roa1,CID1,doc) IMPLICIT none integer*4 I,J real*8 AL(*),G(*),A(*),YDX,YDY,ram1,P1,roa1,CID1,doc,SAL0,SAL1, 1SAG0,SAG1,SPAL,SA1,down,SpA3,beta2,AMU,BOHR,w, 1gpisvejc,betaa,betag,CM,ecm,EXCA AMU=1822.0d0 BOHR=0.529177d0 CM=219474.63d0 ecm=w*CM EXCA=1.0d7/ecm*10.0d0 SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 SPAL=0.0d0 DO 3 I=1,3 SA1=SA1+(AL(1+3*(I-1))*(A(9*(2-1)+3*(3-1)+I)-A(9*(3-1)+3*(2-1)+I)) 1 +AL(2+3*(I-1))*(A(9*(3-1)+3*(1-1)+I)-A(9*(1-1)+3*(3-1)+I)) 1 +AL(3+3*(I-1))*(A(9*(1-1)+3*(2-1)+I)-A(9*(2-1)+3*(1-1)+I)) 1 )/3.0d0 SPAL=SPAL+AL(I+3*(I-1)) DO 3 J=1,3 SAL0=SAL0+AL(I+3*(I-1))*AL(J+3*(J-1)) SAL1=SAL1+AL(I+3*(J-1))*AL(I+3*(J-1)) SAG0=SAG0+AL(I+3*(I-1))* G(J+3*(J-1)) 3 SAG1=SAG1+AL(I+3*(J-1))* G(I+3*(J-1)) c c polarized backscattering components (YDY+YDX=S180): YDY=3.0d0*SAL1-SAL0 YDX=4.0d0*SAL1+2.0d0*SAL0 if(YDX.gt.1.0d-9)then P1=YDY/YDX else P1=0.0d0 endif c calculate degree of circularity down= 7.0d0*SAL1+ 1.0d0*SAL0 doc=5.0d0*( 1.0d0*SAL1- 1.0d0*SAL0) if(down.ne.0.0d0)doc=doc/down SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4)) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) YDY=6.0d0*YDY*(AMU*BOHR**4) YDX=6.0d0*YDX*(AMU*BOHR**4) gpisvejc=(AMU*BOHR**5)*2.0d0*4.0d0*atan(1.0d0)/EXCA c betaa=SA1*3.0d0/2.0d0 *gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0 *gpisvejc c c ICP/SCP(180): c 8/c*(12*beta(G')^2 + 4*beta(A)^2): roa1=96.0d0*(betag+betaa/3.0d0) c ram1 = 6(7 al_ab al_ab + al_aa al_bb) ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) CID1=0.0d0 if(ram1.gt.1.0d-9)CID1=roa1/ram1 return end subroutine rtd(k,alk,ak,gk,AL,A,G) implicit none integer*4 i,j,k,jj,jm,km,kk real*8 alk(9),ak(27),gk(9),AL(*),A(*),G(*) alk=0.0d0 ak=0.0d0 gk=0.0d0 call setij(i,j,k) do 2 jj=1,3 jm=3*(jj-1) alk(i+jm)=alk(i+jm)+AL(j+jm) alk(j+jm)=alk(j+jm)-AL(i+jm) alk(jj+3*(i-1))=alk(jj+3*(i-1))+AL(jj+3*(j-1)) alk(jj+3*(j-1))=alk(jj+3*(j-1))-AL(jj+3*(i-1)) gk( i+jm)=gk( i+jm)+G( j+jm) gk( j+jm)=gk( j+jm)-G( i+jm) gk( jj+3*(i-1))=gk( jj+3*(i-1))+G( jj+3*(j-1)) gk( jj+3*(j-1))=gk( jj+3*(j-1))-G( jj+3*(i-1)) do 2 kk=1,3 km=9*(kk-1) ak(i+jm+km)=ak(i+jm+km)+A(j+jm+km) ak(j+jm+km)=ak(j+jm+km)-A(i+jm+km) ak(jj+3*(i-1)+km)=ak(jj+3*(i-1)+km)+A(jj+3*(j-1)+km) ak(jj+3*(j-1)+km)=ak(jj+3*(j-1)+km)-A(jj+3*(i-1)+km) ak(kk+jm+9*(i-1))=ak(kk+jm+9*(i-1))+A(kk+jm+9*(j-1)) 2 ak(kk+jm+9*(j-1))=ak(kk+jm+9*(j-1))-A(kk+jm+9*(i-1)) return end