PROGRAM DC IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NMAX=160) REAL*8 M DIMENSION A(3),B(3),C(3), UA(NMAX,3), 1CM(NMAX,3),U(3,NMAX,NMAX),M(3,NMAX,NMAX), 2TAB(3,NMAX,NMAX), DTAB(NMAX,NMAX), V(NMAX,NMAX), T(3), 3NS(NMAX),XT(3,NMAX), 4EVAL(NMAX), COEF(NMAX,NMAX), 5OMEGA(NMAX,NMAX),E(NMAX),UM(NMAX,3), 6DM(NMAX,3),AMM(NMAX,3),AMMI(NMAX,3),RRV(3,NMAX) CHARACTER*1 OK PI = 3.141592654D0 WRITE(6,66666)1993 66666 FORMAT(80(1H-),/,15X,'DYNAMIC COUPLING FOR VCD',/,80(1H-),/,/, 1 15X,20H Petr Bour's program,/, 25X,I10, ' version, UOCHB, Prague',/,/, 35X,'INPUT : DC.SC OUTPUT: DC.TAB',/, 55X,' FILE.XYZ DIP.MOM',/, 65X,' EIGEN.VEC',/) OPEN(22,FILE='DC.SC',STATUS='OLD') READ(22,*)NG DO 1 I=1,NG READ(22,*) JD READ(22,*)XT(1,I),XT(2,I),XT(3,I) READ(22,*)NS(I) READ(22,*) READ(22,*)JD,OMEGA(I,1) DO 2 J=2,NS(I) READ(22,*)JD,OMEGA(I,J) 2 OMEGA(I,J)=OMEGA(I,J)-OMEGA(I,1) READ(22,*) DO 1 J=2,NS(I) 1 READ(22,*)JD,U(1,I,J),U(2,I,J),U(3,I,J),M(1,I,J),M(2,I,J),M(3,I,J) CLOSE(22) IOSC=0 DO 4 I=1,NG DO 4 J=2,NS(I) IOSC=IOSC+1 E(IOSC)=OMEGA(I,J)*8065.54093D0 DO 4 K=1,3 CM(IOSC,K)=XT(K,I)*10.0D0 C SC INPUT FILE USES NANOMETERS C WORKING UNITS : DEBYES AND ANGSTROMS UA(IOSC,K)=U(K,I,J) UM(IOSC,K)=M(K,I,J) 4 CONTINUE N=IOSC WRITE(6,*)' ',N,' OSCILLATORS ',NG,' GROUPS' WRITE(6,*) WRITE(6,*)' Relative permitivity:' READ(5,*)EPS0 DO 150 I = 1,N DO 150 J = 1,N DO 150 K = 1,3 150 TAB(K,I,J) = CM(J,K) - CM(I,K) DO 180 I = 1,N DO 180 J = 1,N DO 200 K = 1,3 200 C(K) = TAB(K,I,J) 180 DTAB(I,J) =DSQRT (DOT(C,C)) DO 210 I = 1,N DO 220 J = I,N EPSR=1.0D0 IF ((I .EQ. J).OR.(DTAB(I,J).LT.0.00001))THEN V(I,J) = 0.0D0 ELSE DO 230 K =1,3 B(K) = UA(I,K) C(K) = UA(J,K) 230 T(K) = TAB(K,I,J) DOTIJ = DOT(B,C) DOTIT = DOT(B,T) DOTJT = DOT(C,T) DTABSQ = DTAB(I,J)*DTAB(I,J) V1 = DOTIJ / DTAB(I,J) / DTABSq V2= DOTIT*DOTJT/ DTAB(I,J) / DTABSQ / DTABSQ V(I,J) = V1 - V2*3.0D0 EPSR=EPS0 IF (EPS0.LT.0.0D0) EPSR=-EPS0*DTAB(I,J) ENDIF V(I,J)=V(I,J)*5034.1096D0/EPSR IF (I.EQ.J)V(I,J)=V(I,I)+E(I) C NOW THE UNITS OF THE ENERGY SHOULD BE CM^-1 V(J,I)=V(I,J) 220 CONTINUE 210 CONTINUE 212 WRITE(6,*)' I,J FOR THE H(I,J) {[0,0] TO CONTINUE, [-100,0]', 1',[-200,0] OPTIONS}:' READ(5,*)I,J IF(I.EQ.-100)THEN WRITE(6,*)' 1-2,3-4,5-6 ETC. INTERACTIONS CANCELED' DO 67 I=1,N-1,2 V(I,I+1)=0.0 67 V(I+1,I)=0.0 GOTO 212 ENDIF IF(I.EQ.-200)THEN WRITE(6,*)' ODD-EVEN AND EVEN-ODD INTERACTIONS CANCELED' DO 68 I=1,N DO 68 J=1,I IF(MOD(I-J,2).NE.0)THEN V(I,J)=0.0 V(J,I)=0.0 ENDIF 68 CONTINUE GOTO 212 ENDIF IF(I*J.EQ.0)GOTO 211 WRITE(6,*)V(IABS(I),IABS(J)) IF (I.LT.0) THEN I=IABS(I) J=IABS(J) WRITE(6,*)' NEW V(I,J):' READ(5,*)V(I,J) V(J,I)=V(I,J) ENDIF GOTO 212 211 WRITE(6,*)' DIAGONALIZING THE HAMILTONIAN' IERR=0 ICODE=2 CALL TRED(NMAX,N,V,COEF,EVAL,ICODE,IERR) IF (IERR.NE.0)THEN WRITE(6,*)' DIAGONALIZATION ERROR ***************' STOP ENDIF WRITE(6,*)' DIAGONALIZATION SUCCESSFUL' WRITE(6,*)' CALCULATING SPECTRAL PARAMETERS ..' DO 370 L = 1,N CON=PI*EVAL(L) C CON=PI DO 371 I=1,3 DM (L,I)=0.0 AMM (L,I)=0.0 371 AMMI(L,I)=0.0 DO 370 I = 1,N CLI=COEF(I,L) DO 372 J=1,3 B(J)=UA(I,J) 372 A(J)=CM(I,J) CALL CROSS(A,B,C) DO 370 J=1,3 DM (L,J)=DM (L,J)+CLI*UA(I,J) AMM (L,J)=AMM (L,J)+CLI*CON*C(J) C AMM (L,J)=AMM (L,J)+CLI*CON*C(J)*E(I) 370 AMMI(L,J)=AMMI(L,J)+CLI*UM(I,J) WRITE(6,*)' NORMALIZE THE RESULTS (Y/N/S) ? ' READ(5,'(A)')OK NN=1 IF (OK.EQ.'Y'.OR.OK.EQ.'Y')NN=NG IF (OK.EQ.'S'.OR.OK.EQ.'S')THEN WRITE(6,*)' NN: ' READ(5,*)NN ENDIF Z=1.0/REAL(NN) WRITE(6,*)' DC.TAB ... FILE WITH THE SPECTRUM' OPEN(23,FILE='DC.TAB') WRITE(23,2323) 2323 FORMAT( 1 ' MODE D RTOT RINT', 2/,' CM-1 DEBYE2 10-44 CgS ', 3/,80(1H-)) DTOT=0.0 RTOT=0.0 DO 8 I=1,N EI=EVAL(I) D=0.0 R=0.0 RI=0.0 DO 9 II=1,3 D =D +DM (I,II)*DM(I,II) R =R +AMM (I,II)*DM(I,II) 9 RI=RI+AMMI(I,II)*DM(I,II) DTOT=DTOT+D*Z RTOT=RTOT+(R+RI)*Z 8 WRITE(23,2324)I,EI,D*Z,Z*(R+RI),Z*RI 2324 FORMAT(I4,F10.3,4G16.6) WRITE(23,2325) 2325 FORMAT(80(1H-)) WRITE (23,450) DTOT,RTOT,NN 450 FORMAT ('DSUM=',G16.6,' ','RSUM=',G16.6,/, 1'SPECTRUM OF ALL SYSTEM DEVIDED BY ',I4) CLOSE(23) CALL WREI(NMAX,N,COEF,DM,AMM,AMMI,EVAL,RRV) WRITE(6,*)' >>> PROGRAM TERMINATED <<<' STOP END SUBROUTINE CROSS (U,V,W) IMPLICIT REAL*8(A-H,O-Z) DIMENSION U(3),V(3),W(3) W(1) = U(2)*V(3)-V(2)*U(3) W(2) = U(3)*V(1)-V(3)*U(1) W(3) = U(1)*V(2)-V(1)*U(2) RETURN END SUBROUTINE WREI(NM,N,COEF,DM,AMM,AMMI,EVAL,R) IMPLICIT REAL*8(A-H,O-Z) DIMENSION COEF(NM,NM),DM(NM,3),EVAL(NM),AMM(NM,3),AMMI(NM,3), 1R(3,NM) REAL*8 ME(3),MM(3),AI(3,3),TI(3) WRITE(6,*)' EIGEN.VEC ... FILE WITH THE EIGENVECTORS' OPEN(26,FILE='EIGEN.VEC') WRITE(26,*)' THE NORMALIZED EIGENVECTORS (BASE MODE):' NQ=N M=8 MI=-7 207 MI=MI+8 MJ=M-7 IF(M.GT.NQ)M=NQ WRITE(26,28)(I,I=MJ,M) 28 FORMAT(8X,8(I7,2X)) DO 206 I=1,NQ 206 WRITE(26,27)I,(COEF(J,I),J=MJ,M) 27 FORMAT(I5,3X,8F9.5) IF (M.LT.NQ)THEN M=M+8 GOTO 207 ENDIF CLOSE(26) DO 281 I=1,N DO 281 J=1,3 281 AMM(I,J)=AMM(I,J)+AMMI(I,J) WRITE(6,*)' DIP.MOM ... FILE WITH MODE POLARIZATIONS' OPEN(3,FILE='FILE.XYZ',STATUS='OLD') READ(3,*) IA=0 9999 READ(3,*)ITYPE IF(ITYPE.NE.9999)THEN IA=IA+1 BACKSPACE(3) READ(3,*)ITYPE,R(1,IA),R(2,IA),R(3,IA) GOTO 9999 ENDIF CLOSE(3) CALL INERTIA(R,AI,TI,IA,NM) OPEN(26,FILE='DIP.MOM') WRITE(26,*)' PROJECTIONS OF TRANSITION DIPOLES:' WRITE(26,2627)(TI(I),I=1,3) 2627 FORMAT(3F15.5,/, 1'MODE DX DY DZ [%] / MX MY MZ [%]') DO 2280 I = 1,N DO 39 II=1,3 ME(II)=0.0 MM(II)=0.0 DO 39 J=1,3 MM(II)=MM(II)+AI(J,II)*AMM(I,J) 39 ME(II)=ME(II)+AI(J,II)*DM (I,J) ANORM=ME(1)**2+ME(2)**2+ME(3)**2 BNORM=MM(1)**2+MM(2)**2+MM(3)**2 DO 404 IX=1,3 IF(BNORM.GT.1.0D-38)MM(IX)=MM(IX)**2/BNORM*100.0 404 IF(ANORM.GT.1.0D-38)ME(IX)=ME(IX)**2/ANORM*100.0 2280 WRITE (26,2626)I,EVAL(I),(ME(IX),IX=1,3),(MM(IX),IX=1,3) 2626 FORMAT(I4,F8.1,3F6.1,' /',3F6.1) CLOSE(26) RETURN END FUNCTION DOT(U,V) IMPLICIT REAL*8(A-H,O-Z) DIMENSION U(3),V(3) DOT=0.0 DO 101 I=1,3 101 DOT = DOT+U(I)*V(I) RETURN END SUBROUTINE INERTIA(R,A,TI,N,NMAX) IMPLICIT REAL*8(A-H,O-Z) DIMENSION R(3,NMAX),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 TRED(3,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 TRED(NMAX,N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(NMAX,NMAX),Z(NMAX,NMAX),D(NMAX),E(160) COMMON/DIA/E COMMON/INDE/IND(160) 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 TQL (NMAX,N,Z,D,IERR,ICODE) RETURN 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END SUBROUTINE TQL(NMAX,N,Z,D,IERR,ICODE) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Z(NMAX,NMAX),D(NMAX),E(160) REAL*8 MACHEP COMMON/DIA/E COMMON/INDE/IND(160) 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