program cubpolar implicit none integer*4 n,nat,nc(3),iargc,IERR,frame(4) real*8 p(3,3),I(3,3),UP(3,3),QP(3),UI(3,3),QI(3),cx(12),iso, 1length,cm(3) real*8,allocatable::r(:) integer*4,allocatable::q(:) character*80 fn c write(*,1010) 1010 FORMAT(' Polarizability visualization as Gaussian cube',/, 1 ' Input: gaussian output',/, 1 ' CUBPOLAR.INP',/, 1 ' Output: pol.cub polarizability ellipsoid',/) if(iargc().ne.1)call report('usage" cubpolar ') call getarg(1,fn) call readopt(nc,cx,iso,length,frame) call rg(fn,nat,q,r,p,0) n=3*nat allocate(r(n),q(nat)) call rg(fn,nat,q,r,p,1) call sc(nat,r,nc,cx,length) c moment of charge inertia: call si(nat,q,r,I,cm) CALL TRED12(3,I,UI,QI,2,IERR) call wr(I,QI,' Moment of inertia:') IF(IERR.NE.0)call report('i: diagonalization not ok') CALL TRED12(3,P,UP,QP,2,IERR) call wr(P,QP,' Polarizability:') IF(IERR.NE.0)call report('P: diagonalization not ok') call wrf('I.cub',nat,r,q,nc,cx,UI,QI,iso,cm,length,fn) call wrf('P.cub',nat,r,q,nc,cx,UP,QP,iso,cm,length,fn) call w2('PI.cub',nat,r,q,nc,cx,UP,UI,QP,QI,iso,cm,length,fn) call fr(frame,r,p,' Polarizability:') call fr(frame,r,I,' Moment of inertia:') end function sp(a,b) implicit none real*8 a(3),b(3),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine vp(a,b,c) implicit none real*8 a(3),b(3),c(3) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) return end subroutine norm(a) implicit none real*8 a(3),an,sp an=sp(a,a) if(an.gt.1.0d-20)an=1.0d0/dsqrt(an) a(1)=a(1)*an a(2)=a(2)*an a(3)=a(3)*an return end subroutine wr(t,d,s) implicit none real*8 t(3,3),d(3) character*(*) s write(6,*) write(6,*)s write(6,601)t 601 format(/,' Laboratory frame:',/,/,3(3E12.4,/),/) write(6,602)d 602 format(' Eigenvalues:',/,/,3E12.4,/) return end subroutine sc(nat,r,nc,cx,length) implicit none integer*4 nat,nc(3),ia,ix,jx real*8 cx(12),r(*),xmin(3),xmax(3),xa,length do 3 ix=1,3 xmin(ix)=r(ix+3*(1-1)) xmax(ix)=r(ix+3*(1-1)) do 3 ia=2,nat xa=r(ix+3*(ia-1)) if(xa.lt.xmin(ix))xmin(ix)=xa 3 if(xa.gt.xmax(ix))xmax(ix)=xa if(nc(1).eq.0)then cx=0.0d0 nc=50 do 1 ix=1,3 nc(ix)=50 cx(ix)=xmin(ix)-2.0d0 cx(3+3*(ix-1)+ix)=(xmax(ix)-xmin(ix)+4.0d0)/dble(nc(ix)) 1 write(6,600)ix,nc(ix),xmin(ix),xmax(ix),cx(ix), 1 (cx(3+3*(ix-1)+jx),jx=1,3) 600 format(i2,i4,' xmin:',f9.1,' xmax:',f8.1,' x0:', 1 f9.1,' cv:',3f8.1) endif if(length.lt.0.01d0)length=dsqrt((xmax(1)-xmin(1))**2 1+ (xmax(2)-xmin(2))**2 1+ (xmax(3)-xmin(3))**2)/4.0d0 return end subroutine wrf(cn,nat,r,q,nc,c,u,abc,iso,cm,length,fn) implicit none integer*4 nat,q(*),nc(3),ix,i,iy,iz real*8 r(*),c(12),x(3),xp,yp,zp,u(3,3),abc(3),iso,cm(3),fa, 1length,down real*8,allocatable::ii(:) character*(*) cn character*80 fn fa=iso*(length/abc(1))**2 allocate(ii(nc(3))) open(3,file=cn) write(3,3000)fn,nat,(c(ix),ix=1,3),nc(1),(c(i),i=4,6), 1nc(2),(c(i),i=7,9),nc(3),(c(i),i=10,12) 3000 format(a80,/,' SCF Total Density',/, 13(i5,3f12.6,/),i5,3f12.6) do 2 i=1,nat 2 write(3,3001)q(i),dble(q(i)),(r(ix+3*(i-1)),ix=1,3) 3001 format(i5,4f12.6) do 1 ix=1,nc(1) do 1 iy=1,nc(2) ii=0.0d0 do 3 iz=1,nc(3) x(1)=c(1)+dble(ix-1)*c(4)+dble(iy-1)*c(7)+dble(iz-1)*c(10) x(2)=c(2)+dble(ix-1)*c(5)+dble(iy-1)*c(8)+dble(iz-1)*c(11) x(3)=c(3)+dble(ix-1)*c(6)+dble(iy-1)*c(9)+dble(iz-1)*c(12) xp=+u(1,1)*(x(1)-cm(1))+u(2,1)*(x(2)-cm(2))+u(3,1)*(x(3)-cm(3)) yp=+u(1,2)*(x(1)-cm(1))+u(2,2)*(x(2)-cm(2))+u(3,2)*(x(3)-cm(3)) zp=+u(1,3)*(x(1)-cm(1))+u(2,3)*(x(2)-cm(2))+u(3,3)*(x(3)-cm(3)) down=(xp/abc(1))**2+(yp/abc(2))**2+(zp/abc(3))**2 3 if(down.gt.1.0d-20)ii(iz)=fa/down 1 write(3,4000)ii 4000 format(6E13.5) close(3) return end subroutine w2(cn,nat,r,q,nc,c,u1,u2,abc1,abc2,iso,cm,length,fn) implicit none integer*4 nat,q(*),nc(3),ix,i,iy,iz real*8 r(*),c(12),x(3),xp,yp,zp,u1(3,3),abc1(3),iso,cm(3),fa1, 1length,down,abc2(3),u2(3,3),fa2 real*8,allocatable::ii(:,:) character*(*) cn character*80 fn fa1=dsqrt(iso)*(length/abc1(1))**2 fa2=dsqrt(iso)*(length/abc2(1))**2 allocate(ii(nc(3),2)) open(3,file=cn) write(3,3000)fn,-nat,(c(ix),ix=1,3),2,nc(1),(c(i),i=4,6), 1nc(2),(c(i),i=7,9),nc(3),(c(i),i=10,12) 3000 format(a80,/,' MO coefficients',/, 1i5,3f12.6,i5,3(/,i5,3f12.6)) do 2 i=1,nat 2 write(3,3001)q(i),dble(q(i)),(r(ix+3*(i-1)),ix=1,3) 3001 format(i5,4f12.6) write(3,3002)2,1,2 3002 format(3i5) do 1 ix=1,nc(1) do 1 iy=1,nc(2) ii=0.0d0 do 3 iz=1,nc(3) x(1)=c(1)+dble(ix-1)*c(4)+dble(iy-1)*c(7)+dble(iz-1)*c(10) x(2)=c(2)+dble(ix-1)*c(5)+dble(iy-1)*c(8)+dble(iz-1)*c(11) x(3)=c(3)+dble(ix-1)*c(6)+dble(iy-1)*c(9)+dble(iz-1)*c(12) xp=u1(1,1)*(x(1)-cm(1))+u1(2,1)*(x(2)-cm(2))+u1(3,1)*(x(3)-cm(3)) yp=u1(1,2)*(x(1)-cm(1))+u1(2,2)*(x(2)-cm(2))+u1(3,2)*(x(3)-cm(3)) zp=u1(1,3)*(x(1)-cm(1))+u1(2,3)*(x(2)-cm(2))+u1(3,3)*(x(3)-cm(3)) down=(xp/abc1(1))**2+(yp/abc1(2))**2+(zp/abc1(3))**2 if(down.gt.1.0d-20)ii(iz,1)= fa1/down xp=u2(1,1)*(x(1)-cm(1))+u2(2,1)*(x(2)-cm(2))+u2(3,1)*(x(3)-cm(3)) yp=u2(1,2)*(x(1)-cm(1))+u2(2,2)*(x(2)-cm(2))+u2(3,2)*(x(3)-cm(3)) zp=u2(1,3)*(x(1)-cm(1))+u2(2,3)*(x(2)-cm(2))+u2(3,3)*(x(3)-cm(3)) down=(xp/abc2(1))**2+(yp/abc2(2))**2+(zp/abc2(3))**2 3 if(down.gt.1.0d-20)ii(iz,2)=-fa2/down 1 write(3,4000)((ii(iz,i),i=1,2),iz=1,nc(3)) 4000 format(6E13.5) close(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(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: 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 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: 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 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C 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(*),E(*) 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 C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: 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) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: 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)) C :::::::::: FORM VECTOR :::::::::: 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 C 200 CONTINUE C 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 C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END subroutine report(s) character*(*) s write(6,*)s stop end subroutine readopt(nc,cx,iso,length,frame) implicit none integer*4 ix,nc(3),frame(4) real*8 cx(12),bohr,iso,length character*5 key bohr=0.529177d0 iso=0.0004d0 nc=0 cx=0.0d0 length=0.0d0 frame=0 open(9,file='CUBPOLAR.OPT') 1 read(9,90,end=99,err=99)key 90 format(a5) if(key(1:3).eq.'FRA')read(9,*)frame if(key(1:3).eq.'ISO')read(9,*)iso if(key(1:3).eq.'LEN')read(9,*)length if(key(1:4).eq.'CAGE')then read(9,*)(cx(ix),ix=1,3) read(9,*)nc(1),(cx(ix),ix=4,6) read(9,*)nc(2),(cx(ix),ix=7,9) read(9,*)nc(3),(cx(ix),ix=10,12) do 2 ix=1,12 2 cx(ix)=cx(ix)/bohr endif goto 1 99 close(9) length=length/bohr return end subroutine si(nat,q,r,I,cm) implicit none integer*4 nat,q(*),ix,jx,ia,MENDELEV real*8 r(*),I(3,3),cm(3),qt,ma real*4 amas parameter (MENDELEV=89) dimension 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/ do 1 ix=1,3 cm(ix)=0.0d0 qt=0.0d0 do 11 ia=1,nat ma=dble(amas(q(ia))) qt=qt+ma 11 cm(ix)=cm(ix)+ma*r(ix+3*(ia-1)) 1 cm(ix)=cm(ix)/qt do 2 ix=1,3 do 2 jx=1,3 I(ix,jx)=0.0d0 do 2 ia=1,nat ma=dble(amas(q(ia))) 2 I(ix,jx)=I(ix,jx)+ma*(r(ix+3*(ia-1))-cm(ix)) 1 *(r(jx+3*(ia-1))-cm(jx)) return end subroutine rg(fn,nat,q,r,p,ic) implicit none integer*4 nat,i,ic,q(*),ix,jx real*8 r(*),p(3,3),bohr character*(*) fn character*80 s80 bohr=0.529177d0 open(9,file=fn,status='old') 1 read(9,80,end=88,err=88)s80 80 format(a80) if(s80(27:44).eq.'Input orientation:')then do 2 i=1,4 2 read(9,*) if(ic.eq.0)then i=0 3 read(9,80)s80 if(s80(2:2).ne.'-')then i=i+1 goto 3 endif nat=i write(6,*)nat,' atoms' goto 88 else do 4 i=1,nat 4 read(9,*)q(i),q(i),r(1+3*(i-1)),(r(ix+3*(i-1)),ix=1,3) do 41 i=1,nat do 41 ix=1,3 41 r(ix+3*(i-1))=r(ix+3*(i-1))/bohr write(6,*)'geometry read' endif endif if(s80(3:23).eq.'Exact polarizability:'.and.ic.ne.0)then read(s80(24:71),*)((p(ix,jx),jx=1,ix),ix=1,3) do 5 ix=1,3 do 5 jx=1,ix 5 p(jx,ix)=p(ix,jx) write(6,*)'polarizability read' goto 88 endif goto 1 88 close(9) return end subroutine fr(frame,r,p,s) c tensors in coordinate system defined by 4 atoms implicit none integer*4 frame(4),i,j,k,l,ix,jx real*8 r(*),p(3,3),pp(3,3),u(3,3),a(3),b(3),c(3),t(3) character*(*) s if(frame(1).eq.0)return i=frame(1) j=frame(2) k=frame(3) l=frame(4) do 1 ix=1,3 a(ix)=r(ix+3*(j-1))-r(ix+3*(i-1)) 1 t(ix)=r(ix+3*(l-1))-r(ix+3*(k-1)) call norm(a) call vp(t,a,b) call norm(b) call vp(b,a,c) do 2 ix=1,3 u(1,ix)=a(ix) u(2,ix)=b(ix) 2 u(3,ix)=c(ix) do 3 ix=1,3 do 3 jx=1,3 pp(ix,jx)=0.0d0 do 3 i=1,3 do 3 j=1,3 3 pp(ix,jx)=pp(ix,jx)+u(ix,i)*p(i,j)*u(jx,j) write(6,'(a)')s write(6,601)pp 601 format(/,' Molecular frame:',/,/,3(3E12.4,/),/) return end