PROGRAM DGAUGE C LAST CHANGES 2-6-95 PETR BOUR IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (N0=1100) REAL*8 PP(N0,3,3),II(N0,3,3),JJ(N0,3,3),SV(N0,3*N0,3), 1 PV(N0,3,3),I0(N0,3,3),R(3,N0), 2 FRQ(3*N0),Z(N0),ELD(3),MGD(3),MDN(3),MDD(3),MDT(3), 3 ME(3),MM(3),AI(3,3),TI(3),PQ(N0*3,3),AQ(N0*3,3) INTEGER A,B PI=3.141592653589 FD =SQRT(388.9219) FR = 2586.285D-8 FC=131.14D-8 BOHR=0.52917715 Z0=0.0 WLIM=0.000001 WRITE(*,*)' COMBINATION OF POLARIZABILITY TENSORS' WRITE(*,*)' AND THE S-VECTORS' WRITE(*,6666) 6666 FORMAT(/, 1' INPUT FILES: F .INP ... VIBRATIONAL PARAMETERS',/, 1' FILE.TEN ... TENSORS',/, 3' OUTPUT FILE: DOG.TAB',/, 4' DIP.MOM',/) OPEN(18,FILE='DOG.TAB') OPEN(19,FILE='DIP.MOM') OPEN(15,FILE='FILE.TEN',STATUS='OLD') READ (15,*)NOAT OPEN(17,FILE='F.INP',STATUS='OLD') READ(17,*)NQ,NAC IF(NAC.NE.3*NOAT)THEN WRITE(*,*)' Number of atoms if FILE.TEN and F.INP is different !' STOP ENDIF DO 35 K=1,NOAT READ(17,*)Z(K),(R(I,K),I=1,3) DO 34 I=1,3 34 R(I,K)=R(I,K)/BOHR 35 CONTINUE WRITE(*,*)' Geometry read in' IMDMAX=0 READ(17,*) DO 45 LL=1,NOAT DO 45 I=1,NQ READ (17,*) IAT,IMD,SV(IAT,I,1),SV(IAT,I,2),SV(IAT,I,3) IF (IMD.GT.IMDMAX)IMDMAX=IMD 45 CONTINUE READ(17,*) READ(17,1717)(FRQ(I),I=1,NQ) 1717 FORMAT(6F11.3) WRITE(*,*)' S-matrix read in' C C READ DIPOLE DERIVATIVES = TENSOR P(A,B) C DO 10 L=1,NOAT DO 10 J=1,3 10 READ (15,*) (PP(L,J,I),I=1,3) C C READ TENSORS I(A,B) AND J(A,B) FOR THE 0 0 0 ORIGIN: C (=ELECTRONIC & NUCLEAR TERMS) C DO 220 L=1,NOAT DO 220 J=1,3 220 READ (15,*) (I0(L,J,I),I=1,3) DO 230 L=1,NOAT DO 230 J=1,3 230 READ (15,*) (JJ(L,J,I),I=1,3) C C READ VELOCITY REPRESENTATION OF DIPOLE DERIVATIVES C DO 100 L=1,NOAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) CLOSE(15) WRITE(*,*)' Dipole derivatives read in' DO 140 L=1,NOAT DO 140 A=1,3 DO 140 B=1,3 140 II (L,A,B) = I0(L,A,B) CALL VCDD0(II,PV,PP,R,Z,NOAT,N0) WRITE(*,*)' Tensors transformed' WRITE (18,501) CALL INERTIA(R,AI,TI,NOAT) WRITE (19,5019)(TI(I),I=1,3) 5019 FORMAT( 1' POLARIZATIONS ACCORDING TO THE INERTIA TENSOR AXES',/, 2' ( ',3F15.4, ')',/, 3' MODE ELECTRIC ', 4' MAGNETIC [%] ',/,80(1H-),/) 501 FORMAT (/,' FREQ. ', 1' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2/,'-----------------------------------------------------------') C DO 411 I=1,3*NOAT DO 411 A=1,3 AQ(I,A)=0.0d0 411 PQ(I,A)=0.0d0 DO 41 I=1,NQ W=FRQ(I) ABSW=ABS(W) WR=SQRT(ABSW) DO 56 B=1,3 MGD(B)=Z0 ELD(B)=Z0 MDN(B)=Z0 MDD(B)=Z0 MDT(B)=Z0 C DO 56 L=1,NOAT DO 56 A=1,3 C ELECTRONIC MOMENT TOTAL: ELD(B) = ELD(B) + PP(L,A,B)*SV(L,I,A) C MAGNETIC DIPOLE NUCLEAR SAME FOR DO/GO MDN(B) = MDN(B) + JJ(L,A,B) *SV(L,I,A) C MAGNETIC DIPOLE ELECTRONIC DO MDD(B) = MDD(B) + II(L,A,B)*SV(L,I,A) C MAGNETIC DIPOLE TOTAL DO MDT(B) = MDT(B) + (JJ(L,A,B)+ II(L,A,B))*SV(L,I,A) C MAGNETIC DIPOLE TOTAL GO 56 MGD(B) = MGD(B) + (I0(L,A,B)+JJ(L,A,B)) *SV(L,I,A) c c APT,AAT in normal modes: do 561 B=1,3 AQ(I,B)=MDT(B) 561 PQ(I,B)=ELD(B) C DO 51 B=1,3 IF(ABSW.GT.WLIM)ELD(B)=ELD(B)/WR*FD MDT(B)=MDT(B)*FC*WR MDN(B)=MDN(B)*FC*WR MDD(B)=MDD(B)*FC*WR 51 MGD(B)=MGD(B)*FC*WR C SUMD = ELD(1)*ELD(1)+ELD(2)*ELD(2)+ELD(3)*ELD(3) RDO = MDT(1)*ELD(1)+MDT(2)*ELD(2)+MDT(3)*ELD(3) RGO = MGD(1)*ELD(1)+MGD(2)*ELD(2)+MGD(3)*ELD(3) IMOD=IMDMAX-I+1 C DO 39 III=1,3 ME(III)=Z0 MM(III)=Z0 DO 39 J=1,3 MM(III)=MM(III)+AI(J,III)*MDT(J) 39 ME(III)=ME(III)+AI(J,III)*ELD(J) C CALL NORM(ME) CALL NORM(MM) C DO 4 III=1,3 MM(III)=MM(III)*MM(III)*100.0 4 ME(III)=ME(III)*ME(III)*100.0 C WRITE(19,1920)IMOD,W,(ME(IR),IR=1,3),(MM(IR),IR=1,3) 1920 FORMAT(I4,F7.1,2(3F10.0,' ')) WRITE(17,522)(ELD(IR),IR=1,3),(MDT(IR),IR=1,3) 522 FORMAT(6G13.5) IF(ABSW.LT.WLIM)SUMD=Z0 41 WRITE (18,523) IMOD,W,SUMD,RDO,RGO,SQRT(2.0*PI)*W*SUMD 523 FORMAT(I4,5G15.6) C WRITE(18,1819) WRITE(19,1819) 1819 FORMAT('-------------------------------------------------') CLOSE(17) CLOSE(18) CLOSE(14) CLOSE(19) call WRITEQTEN(N0,PQ,AQ,NOAT) WRITE(*,*)' >> PROGRAM TERMINATED <<' STOP END SUBROUTINE VCDD0(VCD,VEL,DIP,C,Z,NAT,N00) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT REAL*8(A-H,O-Z) INTEGER ALPHA,BETA,GAMMA,DELTA,LAMBDA PARAMETER (N0=1100) DIMENSION VCD(N00,3,3),VEL(N00,3,3),DIP(N00,3,3), 1C(3,N00),Z(N00) DO 90 BETA =1,3 DO 90 GAMMA=1,3 DO 90 DELTA=1,3 SKEW=0.25D0*EPS(BETA,GAMMA,DELTA) IF (DABS(SKEW).GT.1.0D-10)THEN DO 80 ALPHA=1,3 DO 80 LAMBDA=1,NAT 80 VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA) 1 -SKEW*C(GAMMA,LAMBDA) 2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA)) ENDIF 90 CONTINUE RETURN END FUNCTION EPS(I,J,K) REAL*8 EPS EPS=0.0D0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0D0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0D0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0D0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0D0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0D0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0D0 RETURN END SUBROUTINE INERTIA(R,A,TI,N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (N0=1100) DIMENSION R(3,N0),A(3,3),T(3,3),TI(3),CM(3) DO 2 I=1,3 CM(I)=0.0 DO 3 J=1,N 3 CM(I)=CM(I)+R(I,J) 2 CM(I)=CM(I)/N DO 1 I =1,3 DO 1 J =1,3 T(I,J)=0.0 DO 1 IS=1,N 1 T(I,J)=T(I,J)+(R(I,IS)-CM(I))*(R(J,IS)-CM(J)) CALL TRED12(3,T,A,TI,2,IERR) RETURN C DEL(I,J)=SUM(OVER JS) OF A(I,JS)*A(J,JS) C DIAGONAL TD(I,J) = A(IS,I)*T(IS,JS)*A(JS,j) END SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (N0=1100) DIMENSION A(N,N),Z(N,N),D(N),E(3*N0) COMMON/DIAG/E COMMON/INDEX/IND(3*N0) 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) RETURN 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END SUBROUTINE TQL12(N,Z,D,IERR,ICODE) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (N0=1100) DIMENSION Z(N,N),D(N),E(3*N0) REAL*8 MACHEP COMMON/DIAG/E COMMON/INDEX/IND(3*N0) 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 GO TO 1001 1000 IERR = L 1001 RETURN END SUBROUTINE NORM(A) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3) AN=DSQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3)) IF(ABS(AN).LT.1.0E-30)RETURN A(1)=A(1)/AN A(2)=A(2)/AN A(3)=A(3)/AN RETURN END SUBROUTINE WRITEQTEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(N0*3,3),A(N0*3,3) OPEN(15,FILE='FILE.Q.TEN') WRITE(15,1500) NAT,NAT-6,0 1500 FORMAT(3I5,' normal mode AAT AAP derivatives') DO 10 L=1,NAT*3 10 WRITE(15,1501) (P(L,I),I=1,3),L 1501 FORMAT(3F14.8,I5) DO 221 L=1,NAT*3 221 WRITE(15,1501) (A(L,I),I=1,3),L DO 230 L=1,NAT*3 230 WRITE(15,1501) (Z0,I=1,3),L DO 100 L=1,NAT*3 100 WRITE(15,1501) (P(L,I),I=1,3),L CLOSE(15) RETURN END