PROGRAM FTRY92 IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) DIMENSION BB(MX,MX),ZM(MX),ZL(MX,MX),QI(MX),MODE(MX),AMODE(MX), 1IFF((MX*(MX+1))/2),ZIM(MX/3,5),V(MX,MX),BD(MX,MX),PE(MX), 2TITLE(20),XXX((MX*(MX+1))/2),AS(3,MX),ZINTY(MX),CHAR(MX/3),IPE(MX) CHARACTER*4 CODE,FTRY,ATOM,RMSA,INTY CHARACTER*1 OK,OK1,OK2 COMMON /ITERA/AMULT(MX),NCORRESP(MX),X((MX*(MX+1))/2),G(MX,MX), 1Q(MX),ZNU(MX,5),Z(MX,MX),P(MX),A(MX,MX),AA(MX,MX), 2Y(MX,MX),FF(MX,MX),SCFAC(MX),DEL COMMON /ITERI/NFC,NOB,NOAT,NOSC COMMON /INDEX/IND(MX) COMMON /DIAG/E(MX) DATA MAXNOB/435/,IRMS/0/,IERR/0/ FTRY='FTRY' ATOM='ATOM' RMSA='RMSA' INTY='INTY' DEL=0.05D0 WRITE(*,*)' FORCE FIELD REFINEMENT PROGRAM' WRITE(*,*)' PC version' WRITE(*,*)' Petr Bour, Prague 9-92' OPEN(16,FILE='FTRY.OUT') WRITE(16,*)' FTRY OUTPUT FILE' WRITE(16,*)' ----------------' DO 22 I=1,435 22 SCFAC(I)=1.0D0 ITRY=0 OPEN(15,FILE='FTRY.INP') 10 READ(15,1070,END=690)CODE,IPB,IPF,IPG,IPL,IPP,IPA,IPS 1,MODMIN,MODMAX OK1='Y' IF(CODE.EQ.'AUTO')THEN WRITE(*,*)' AUTOMATIC RUN ' CODE='INTY' OK1='N' ENDIF IF(CODE.EQ.'RMSA')GO TO 499 c OPEN(8,FILE='INTY.FC',STATUS='OLD') c READ(8,*) c READ(8,*)N c NFC=N c READ(8,555)(X(I),I=1,N) c555 FORMAT(6F12.6) c CLOSE(8) c c Use rather the binary stream for FC: OPEN(8,FILE='INTY.SCR',STATUS='OLD',FORM='UNFORMATTED') READ(8)N,(X(I),I=1,N) CLOSE(8) NFC=N open(10,file='0.txt') do i=1,N write(10,*)i,x(i) enddo close(10) c c Read-in binary B-mat c OPEN(9,FILE='B.MAT',STATUS='OLD',FORM='UNFORMATTED') c READ(9)N,M,((BB(I,J),J=1,M),I=1,N) c CLOSE(9) c NOB=N c NA=M OPEN(9,FILE='B.MAT') read(9,*)N read(9,*)M NOB=N NA=M NOAT=NA/3 do 14 i=1,NOB do 14 j=1,NA 14 BB(i,j)=0.0d0 read(9,*) do 13 I=1,NOB read(9,*) read(9,*)nn do 13 ii=1,nn 13 read(9,*)ia,ix,BB(I,ix+3*(ia-1)) CLOSE(9) write(6,*)'B.MAT read',NOAT,NOB OPEN(8,FILE='INTY.FC') read(8,*) N=0 N1=1 199 N3=N1+4 IF(N3.GT.NOB)N3=NOB read(8,*) DO 1301 LN=N1,NOB c read(8,*)LNDUM,(X(N+J-N1+1),J=N1,MIN(LN,N3)) read(8,*)LNDUM,(X((LN*(LN-1))/2+J),J=N1,MIN(LN,N3)) 1301 N=N+MIN(LN,N3)-N1+1 N1=N1+5 IF(N3.LT.NOB)GOTO 199 close(8) write(6,*)N,' force constants in INTY.FC' NFC=N open(10,file='1.txt') do i=1,N write(10,*)i,x(i) enddo close(10) c NTOT=N READ(15,1020) NOSC READ(15,1040) (SCFAC(I),I=1,NOSC) READ(15,1020) (NCORRESP(I),I=1,NOSC) WRITE(*,*)NOSC,' scale factors' NOAT=NA/3 DO 7 I=1,NOAT 7 CHAR(I)=0.0D0 CHATOT=0.0D0 IF(CODE.EQ.'INTY') READ(15,1040) (CHAR(I),I=1,NOAT) DO 77 I=1,NOAT 77 CHATOT=CHATOT+CHAR(I) WRITE(*,777)CHATOT 777 FORMAT(' TOTAL CHARGE',F10.6) READ(15,1020)NM DO 20 J=1,NM 20 READ(15,1040)(ZIM(I,J),I=1,NOAT) ISW=0 JMIN=1 DO 50 I=1,NOB DO 40 J=1,NOB 40 FF(I,J)=0.D0 JMAX=JMIN+I-1 DO 41 J=JMIN,JMAX 41 IFF(J)=J 50 JMIN=JMAX+1 60 READ(15,1090)I,J,T IF(I.LE.0)GO TO 85 IF(IPF.NE.0)GO TO 65 IF(ISW.EQ.0)WRITE(16,1100)I,J,T IF(ISW.EQ.1)WRITE(16,1101)I,J,T IF(ISW.EQ.2)WRITE(16,1102)I,J,T ISW=ISW+1 IF(ISW.GT.2)ISW=0 65 IF(J.LE.I.AND.I.LE.NOB)GO TO 80 70 WRITE(16,1110) STOP 80 FF(I,J)=T FF(J,I)=T K=I*(I-1)/2+J IF(IFF(K).NE.0)GO TO 70 IFF(K)=29999 GO TO 60 85 JMAX=NOB*(NOB+1)/2 DO 150 J=1,NM 150 READ(15,1040)(ZNU(I,J),I=1,NOB) CLOSE(15) WRITE(16,*) NOB,NOAT WRITE(16,*)' NUMBER OF MOLECULES: ',NM DO 21 J=1,NM DO 21 I=1,NOAT WRITE(16,1271) ZIM(I,J) 21 ZIM(I,J)=1.D0/ZIM(I,J) DO 90 J=1,JMAX 90 IF(IFF(J).EQ.29999)IFF(J)=0 IF(IPF.NE.0)GO TO 111 KMIN=1 KMAX=MIN0(NOB,30) IMIN=0 95 CONTINUE DO 100 I=KMIN,NOB IMIN=IMIN+I-1 JMIN=IMIN+KMIN 100 JMAX=IMIN+MIN0(I,KMAX) IMIN=(KMAX-1)*KMAX/2 KMIN=KMAX+1 KMAX=KMAX+MIN0(NOB-KMAX,30) IF(KMIN.LE.NOB)GO TO 95 WRITE(16,1130) 111 KKK=1 ITRY=ITRY+1 WRITE(16,*)' TRIAL ',ITRY N=NFC WRITE(16,*) ' THE SCALE FACTORS' WRITE(16,*)' COORD SF# SF FROM TO' DO 131 I=1,NOB AMULT(I)=SCFAC(KKK) IF(KKK.EQ.1)WRITE(16,1616) I,KKK,SCFAC(KKK),1,NCORRESP(KKK) IF(KKK.GT.1)WRITE(16,1616) I,KKK,SCFAC(KKK),NCORRESP(KKK-1), 1NCORRESP(KKK) 1616 FORMAT(2I4,G15.6,I5,' - ',I5) 131 IF(I.EQ.NCORRESP(KKK)) KKK=KKK+1 WRITE(16,1020) NOSC WRITE(16,1040) (SCFAC(I),I=1,NOSC) WRITE(16,1020) (NCORRESP(I),I=1,NOSC) K=0 DO 130 I=1,NOB DO 130 J=1,I K=K+1 L=IFF(K) IF(L.EQ.0)GO TO 130 LL=IABS(L) IF(LL.LE.N)GO TO 120 WRITE(16,1140)I,J,LL,N,IFF(K) STOP 120 T=X(LL)*DSQRT(AMULT(I)*AMULT(J)) XXX(LL)=T IF(L.LT.0)T=-T FF(I,J)=T FF(J,I)=T Z(I,J)=T Z(J,I)=T 130 CONTINUE OPEN(35,FILE='SCALED.FC') WRITE(35,*)' SCALED FORCE CONSTANTS' NQ=INT( (-1.0+SQRT(REAL(1+8*N)))*0.5 ) WRITE(35,3337)N,NQ,NQ 3337 FORMAT(/,I6,' INTRINSIC FORCE CONSTANTS',I5,' x (',I5,' + 1) / 2') M=6 MI=-5 207 MI=MI+6 MJ=M-5 IF(M.GT.NQ)M=NQ WRITE(35,28)(I,I=MJ,M) 28 FORMAT(8X,6(I7,5X)) DO 206 I=MI,NQ MW=M IF(I.LT.MW)MW=I 206 WRITE(35,27)I,(FF(J,I),J=MJ,MW) 27 FORMAT(I5,3X,6F12.6) IF (M.LT.NQ)THEN M=M+6 GOTO 207 ENDIF CLOSE(35) IRMS=1 DO 1 I=1,180 P(I)=0.0D0 1 IND(I)=I IF(IPF.NE.0)GO TO 140 WRITE(16,1200) CALL MATOUT(FF,MAXNOB,NOB,NOB) 140 CONTINUE CALL TRED12(NOB,FF,Z,P,2,IERR) IF(IERR.EQ.0)GO TO 200 190 WRITE(*,1220)IERR STOP 200 CONTINUE DO 490 L=1,NM DO 210 I=1,NOB T=P(I) IF(T.LT.1.D-8)T=0.D0 T=DSQRT(T) IF(DABS(T).GT.0.0D-10)P(I)=1.0D0/T DO 210 J=1,NOB 210 A(I,J)=T*Z(J,I) OPEN(3,FILE='SCRATCH.3',STATUS='UNKNOWN',FORM='UNFORMATTED') REWIND 3 WRITE(3) ((A(II,JJ),II=1,NOB),JJ=1,NOB) NF=0 FWAVE=0.D0 DO 220 I=1,NOAT T=ZIM(I,L) J=3*I ZM(J)=T ZM(J-1)=T 220 ZM(J-2)=T DO 240 I=1,NOB DO 240 J=I,NOB TP=0.D0 DO 230 K=1,NA 230 TP=TP+BB(I,K)*BB(J,K)*ZM(K) G(I,J)=TP 240 G(J,I)=TP DO 270 I=1,NOB DO 270 J=1,NOB T=0.D0 DO 260 K=1,NOB 260 T=T+G(I,K)*A(J,K) 270 AA(I,J)=T DO 290 I=1,NOB DO 290 J=1,I T=0.D0 DO 280 K=1,NOB 280 T=T+A(I,K)*AA(K,J) 290 Y(I,J)=T WRITE(16,1210) IF(CODE.NE.'INTY') WRITE(16,1280) L DO 2 I=1,NOB 2 IND(I)=NOB-I+1 CALL TRED12(NOB,Y,Y,Q,-2,IERR) IF(IERR.NE.0)GO TO 190 RMSW=0.0D0 DO 310 I=1,NOB T=Q(I) IF(DABS(T).LE.1.D-8)T=0.D0 QI(I)=0.D0 SIGN=1.0D0 IF (T.GT.1.0D-8) SIGN=T/DABS(T) T=DABS(T) IF(DABS(T).GT.1.0D-10)QI(I)=1.0D0/T T=DSQRT(T)*1302.828D0 Q(I)=T TP=ZNU(I,L) NF=NF+1 QQ=T-TP IF (TP.GT.0.0001)RMSW=RMSW+QQ*QQ FWAVE=FWAVE+QQ*QQ T=SIGN*T IF(CODE.NE.'INTY') WRITE(16,1230) I, T,TP,QQ AS(1,I)=T AS(2,I)=TP AS(3,I)=QQ 310 CONTINUE OPEN(11,FILE='SCRATCH.1',STATUS='UNKNOWN',FORM='UNFORMATTED') REWIND(11) WRITE(11) (Q(I),I=1,NOB) FWAVE=DSQRT(FWAVE/DFLOAT(NF)) RMSW=DSQRT(RMSW/DFLOAT(NF)) IF (CODE.EQ.'FTRY')WRITE(16,1370)FWAVE WRITE(16,*)' DEVIATION FROM NON-ZERO FREQUENCIES:' WRITE(16,1370)RMSW SO1=0.0D0 SO2=0.0D0 DO 5151 J=1,NOB SO1=SO1+AS(1,J)*AS(2,J) 5151 SO2=SO2+AS(1,J)*AS(1,J) AK=SO1/SO2 WRITE(*,5555) 5555 FORMAT(1X,'************ NEW FREQUENCIES ******************') ILINE=0 DO 311 J=NOB,1,-1 ILINE=ILINE+1 IF (ILINE.EQ.20)THEN WRITE(*,*)' INPUT ANY CHARACTER TO CONTINUE ...' IF(OK1.NE.'N')READ(*,'(A1)')OK ILINE=0 ENDIF 311 WRITE(*,6666)J,NOB-J+1,(AS(I,NOB-J+1),I=1,3) 6666 FORMAT(1X, 1 I5,' (',I3,') ...... ',F8.2,' -',F8.2,'/',F8.2) WRITE(*,5551) 5551 FORMAT(1X,'************************************************') WRITE(*,6667)FWAVE,RMSW,AK 6667 FORMAT(1X,'>>>>>> ',F8.2,' << >> ',F8.2,' <<<<<< ',F8.6,/, 11X,' NEW SCALEFACTORS (Y/N) ?') OK='N' IF(OK1.NE.'N')READ(*,'(A)')OK IF (OK.EQ.'Y'.OR.OK.EQ.'y') THEN DO 601 I=1,NOSC WRITE(*,6668)I,SCFAC(I) 6668 FORMAT(1X,'NEW VALUE FOR THE ',I3,'. SCALE FACTOR (OLD ', 1 f5.3,' ) :') 601 READ(*,*)SCFAC(I) GOTO 111 ENDIF WRITE(*,*)' AUTOMATIC REFINEMENT (Y/N) ?' OK='N' IF(OK1.NE.'N')READ(*,'(A)')OK IF (OK.EQ.'Y'.OR.OK.EQ.'y')THEN 333 WRITE(*,*)' MAXIMUM NUMBER OF ITERATIONS:' READ(*,*)NIT IF (NIT.LT.0)THEN NIT=-NIT WRITE(*,*)' STEP:' READ(*,*)DEL ENDIF WRITE(*,*)' TYPE OF THE FITTING PROCEDURE (1,2): ' READ(*,*)ITYPE IF(ITYPE.NE.2)CALL IMPROVE(NIT,S0) IF(ITYPE.EQ.2)CALL IMPROVE2(NIT,S0) WRITE(*,6001)S0 6001 FORMAT(' Average error: ',F12.6) WRITE(*,*)' AGAIN (Y/N) ?' READ(*,'(A)')OK IF (OK.EQ.'y'.OR.OK.EQ.'Y')GOTO 333 GOTO 111 ENDIF C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ DO 340 J=1,NOB T=P(J) DO 340 I=1,NOB 340 Z(I,J)=Z(I,J)*T DO 360 I=1,NOB DO 360 J=1,NOB T=0.D0 DO 350 K=1,NOB 350 T=T+Z(I,K)*Y(K,J) A(I,J)=T 360 ZL(I,J)=T CALL INVERT(A,NOB,TP,IERR,I,J) IF(IERR.EQ.0)GO TO 370 WRITE(16,1060) GO TO 410 370 DO 400 I=1,NOB TP=0.D0 DO 390 K=1,NOB SUM=0.D0 DO 380 M=1,NOB 380 SUM=SUM+G(K,M)*A(I,M) 390 TP=TP+A(I,K)*SUM IF(TP.LT.1.D-10)TP=0.D0 TP=DSQRT(TP) DO 400 J=1,NOB 400 ZL(J,I)=ZL(J,I)*TP IF(IPL.NE.0)GO TO 410 IF(CODE.EQ.'ATOM') WRITE(16,1170) IF(CODE.EQ.'ATOM') CALL MATOUT(ZL,MAXNOB,NOB,NOB) IF(CODE.EQ.'FTRY') WRITE(16,1170) IF(CODE.EQ.'FTRY') CALL MATOUT(ZL,MAXNOB,NOB,NOB) 410 DO 430 I=1,NOB T=QI(I) DO 430 J=1,NOB SUM=0.D0 DO 420 K=1,NOB 420 SUM=SUM+ZL(K,I)*FF(K,J) 430 A(I,J)=SUM*T IF(IPP.NE.0)GO TO 460 DO 450 J=1,NOB SUM=0.D0 DO 440 I=1,NOB TP=ZL(I,J) TP=FF(I,I)*TP*TP SUM=SUM+TP 440 AA(I,J)=TP IF(DABS(SUM).GT.1.D-5)SUM=1.D0/SUM DO 450 I=1,NOB 450 AA(I,J)=AA(I,J)*SUM IF(CODE.EQ.'ATOM')WRITE(16,1250) 1250 FORMAT(////'0POTENTIAL ENERGY DISTRIBUTION MATRIX (INTERNAL ', $ 'COORDS BY NORMAL COORDS):') IF(CODE.EQ.'ATOM')CALL MATOUT(AA,MAXNOB,NOB,NOB) OPEN(19,FILE='POT.EN',STATUS='UNKNOWN') WRITE(19,*)' POTENTIAL ENERGY DISTRIBUTION [in %]' WRITE(19,*)' MODE coordinate(PED) (for PED > 5%)' WRITE(19,1920) 1920 FORMAT(80(1H-)) DO 454 J=1,NOB AMODE(J)=0.0D0 454 MODE(J)=0 c c Loop over normal modes: DO 451 J=1,NOB NPE=0 c c Loop over coordinates DO 453 II=1,NOB AMAX=0.0D0 IMAX=0 c c Second loop over coordinates - find maximum and zero after DO 452 I=1,NOB IF(AA(I,J).GT.AMAX)THEN AMAX=AA(I,J) IMAX=I ENDIF 452 CONTINUE AA(IMAX,J)=0.0D0 AMAX=AMAX*100.0D0 c IF (AMAX.GE.5.0D0)THEN NPE=NPE+1 PE(NPE)=AMAX IPE(NPE)=IMAX ENDIF c c record the absolute maximum:(mode J, coord IMAX) IF(AMODE(IMAX).LT.AMAX)THEN MODE(IMAX)=J AMODE(IMAX)=AMAX ENDIF 453 CONTINUE NPE=MIN(NPE,5) IF(NPE.EQ.1)WRITE(19,1922)NOB-J+1,AS(1,J),IPE(1),PE(1) IF(NPE.EQ.2)WRITE(19,1922)NOB-J+1,AS(1,J),IPE(1),PE(1), 1IPE(2),PE(2) IF(NPE.EQ.3)WRITE(19,1922)NOB-J+1,AS(1,J),IPE(1),PE(1), 1IPE(2),PE(2),IPE(3),PE(3) IF(NPE.EQ.4)WRITE(19,1922)NOB-J+1,AS(1,J),IPE(1),PE(1), 1IPE(2),PE(2),IPE(3),PE(3),IPE(4),PE(4) IF(NPE.EQ.5)WRITE(19,1922)NOB-J+1,AS(1,J),IPE(1),PE(1), 1IPE(2),PE(2),IPE(3),PE(3),IPE(4),PE(4),IPE(5),PE(5) 1922 FORMAT(I5,F9.1, 5(I5,'(',F4.1,')') ) 451 CONTINUE WRITE(19,1920) DO 455 I=1,NOB 455 WRITE(19,1919)I,NOB+1-MODE(I),AMODE(I) 1919 FORMAT(' Coordinate',I4,' is most present in mode',I4,' (by', 1F5.1,'%).') WRITE(19,1920) CLOSE(19) 460 IF(IPA.EQ.2)GO TO 490 DO 480 I=1,NA T=ZM(I) DO 480 J=1,NOB SUM=0.D0 DO 470 K=1,NOB 470 SUM=SUM+BB(K,I)*A(J,K) 480 AA(I,J)=SUM*T WRITE(*,*)' S -MATRIX COMPUTED' CALL SMATOUT(AA,AS,NOB,NA,MODMIN,MODMAX,OK1) IF(CODE.EQ.'ATOM') WRITE(16,1240) IF(CODE.EQ.'ATOM')CALL MATOUT(AA,MAXNOB,NA,NOB) DO 485 I=1,NOB ZMUX = 0.0D0 ZMUY = 0.0D0 ZMUZ = 0.0D0 DO 484 J=1,NOAT JJX = 3*J - 2 JJY = 3*J - 1 JJZ = 3*J ZMUX = ZMUX + CHAR(J)*AA(JJX,I) ZMUY = ZMUY + CHAR(J)*AA(JJY,I) ZMUZ = ZMUZ + CHAR(J)*AA(JJZ,I) 484 CONTINUE ZINTY(I) = ZMUX*ZMUX + ZMUY*ZMUY + ZMUZ*ZMUZ 485 CONTINUE OPEN(2,FILE='SCRATCH.2',STATUS='UNKNOWN',FORM='UNFORMATTED') WRITE(2) (ZINTY(I),I=1,NOB) CALL FPC(Q,CHAR,AA,NOAT,NOB) IF(IPA.EQ.0.OR.CODE.EQ.'FTRY'.OR.CODE.EQ.'INTY')GOTO 490 OPEN(7,FILE='AA.PLT') WRITE(7,1400) NA,NOB DO 489 J=1,NOB DO 488 I=1,NA IF(DABS(AA(I,J)).LE.1.0D-10) AA(I,J) = 0.0D0 488 CONTINUE WRITE(7,1410) (AA(I,J),I=1,NA) 489 CONTINUE CLOSE(7) 490 CONTINUE IF(CODE.EQ.'ATOM') GO TO 499 IF(CODE.EQ.'FTRY')STOP IF(CODE.NE.'INTY') GO TO 499 REWIND 2 BGZTY = 0.0D0 DO 492 I=1,NM READ(2) (ZINTY(J),J=1,NOB) DO 492 J=1,NOB 492 BGZTY = DMAX1(BGZTY,ZINTY(J)) FAC=1.0D0 IF (DABS(BGZTY).GT.1.0D-10) FAC = 10.0D0/BGZTY REWIND 11 REWIND 2 REWIND 3 READ(3) ((A(I,J),I=1,NOB),J=1,NOB) NF = 0 FWAVE = 0.0D0 WRITE(16,1420) DO 494 L=1,NM WRITE(16,1280) L READ(2) (ZINTY(J),J=1,NOB) READ(11)(Q(I),I=1,NOB) DO 494 J=1,NOB ZINTY(J) = ZINTY(J)*FAC IF(DABS(ZNU(J,L)).LT.1.0D-10) GO TO 493 QQ = Q(J) - ZNU(J,L) NF = NF + 1 FWAVE = FWAVE + QQ*QQ WRITE(16,1440) Q(J),ZNU(J,L),QQ,ZINTY(J) GO TO 494 493 WRITE(16,1430) Q(J),ZINTY(J) 494 CONTINUE IF(NF.GT.0) FWAVE = DSQRT(FWAVE/DFLOAT(NF)) IF(NF.GT.0) WRITE(16,1370) FWAVE 499 IF(IPS.LE.0.OR.CODE.EQ.'INTY') STOP IF(IPS.EQ.1)READ(15,1000)TITLE,M,L,((BD(I,J),J=1,NA),I=1,M) IF (CODE.EQ.'ATOM')READ(15,1040)TK IF(IERR.NE.0)GO TO 10 IF(IRMS)520,680,500 500 DO 510 I=1,NOB QQ=Q(I) IF(DABS(QQ).LT.1.0D-10)GO TO 510 T=DEXP(0.719418D0*QQ/TK) TP=1.D0/T Q(I)=(16.858039D0/QQ)*((T+TP)/(T-TP)) 510 CONTINUE WRITE(16,1330) 520 GO TO (580,530,560),IPS 530 WRITE(16,1160) M=NOB DO 550 I=1,M T=0.D0 DO 540 K=1,NOB TP=ZL(I,K) 540 T=T+TP*TP*Q(K) 550 P(I)=DSQRT(T) GO TO 670 560 WRITE(16,1300) M=NA DO 570 I=1,M T=ZM(I) DO 570 J=1,NOB 570 V(I,J)=T*BB(J,I) GO TO 620 580 WRITE(16,1340)TITLE CALL MATOUT(BD,MAXNOB,M,NA) DO 590 J=1,NA T=ZM(J) DO 590 I=1,M 590 BD(I,J)=BD(I,J)*T DO 610 I=1,M DO 610 J=1,NOB T=0.D0 DO 600 K=1,NA 600 T=T+BD(I,K)*BB(J,K) 610 V(I,J)=T 620 DO 640 I=1,M DO 640 J=1,NOB T=0.D0 DO 630 K=1,NOB 630 T=T+V(I,K)*A(J,K) 640 AA(I,J)=T DO 660 I=1,M T=0.D0 DO 650 K=1,NOB TP=AA(I,K) 650 T=T+TP*TP*Q(K) 660 P(I)=DSQRT(T) 670 WRITE(16,1320)TK,(I,P(I),I=1,M) IRMS=-1 GO TO 10 680 WRITE(16,1310) GO TO 10 690 WRITE(16,1260) CLOSE(7) CLOSE(16) CLOSE(17) CLOSE(11) CLOSE(2) CLOSE(3) STOP 1000 FORMAT(10A8/10A8/2I4/(9F14.10)) 1020 FORMAT(20I4) 1040 FORMAT(6G12.6) 1060 FORMAT('0INVERSION OF L MATRIX FAILED: CAUSE IS LINEAR', $ ' DEPENDENCE OF THE EIGENVECTORS') 1070 FORMAT(A4,9I4) 1090 FORMAT(2I4,G12.6) 1100 FORMAT('0F(',I3,',',I3,') FIXED AT',G15.7) 1101 FORMAT('+',44X,'F(',I3,',',I3,') FIXED AT',G15.7) 1102 FORMAT('+',88X,'F(',I3,',',I3,') FIXED AT',G15.7) 1110 FORMAT('0CONSTANT VALUE FOR F(I,J) IS: NOT IN LOWER TRIANGLE, ', $ 'BOTH A CONSTANT AND A PARAMETER, OR DEFINED TWICE'/) 1130 FORMAT('0IFF(I,J)>0: INSERT PARAMETER NUMBER SHOWN' 1/9X,'=0: F(I,J)', 2 'J) IS A CONSTANT'/9X,'<0: INSERT - OF PARAMETER NUMBER SHOWN') 1140 FORMAT('0* ELEMENT (',I5,',',I5,') OF F MATRIX REQUIRES NON-', $ 'EXISTENT FORCE CONSTANT NO',I5//) 1160 FORMAT(/'0B MATRIX USED AS BD MATRIX, FOR RMS AMPLITUDES OF ', $ 'THE INTERNAL COORDINATES') 1170 FORMAT(////'0G-NORMALIZED L MATRIX (INTERNAL COORDS BY NORMAL ', $ 'COORDS):') 1200 FORMAT(//'0F MATRIX:') 1210 FORMAT(/////'0 CALCULATED OBSERVED DIFFERENCE'/) 1220 FORMAT('0*** DIAGONALIZATION ERROR *** IERR =',I4) 1230 FORMAT(2X,I4,3F15.2) 1240 FORMAT(////'0ATOM DISPLACEMENT MATRIX (NA BY NORMAL COORDS):') 1260 FORMAT('1*** NORMAL TERMINATION'//) 1271 FORMAT(1X,6F12.6) 1280 FORMAT('0MOLECULE',I4/) 1300 FORMAT(/'0UNIT MATRIX USED AS BD MATRIX, FOR RMS CARTESIAN ', $ 'DISPLACEMENTS') 1310 FORMAT('0ILLEGAL USE OF RMSA OPTION (NO PREVIOUS ATOM RUN)') 1320 FORMAT(/'0RMS AMPLITUDES (ANGSTROMS) AT',F8.2,' DEG KELVIN:'/ $ ('0',6(I3,G15.6,4X))) 1330 FORMAT(///) 1340 FORMAT(//2(1X,10A8/)/' BD MATRIX (ATOM PAIRS BY NA):') 1370 FORMAT(///,' RMS ERROR (WAVENUMBERS):',G13.4) 1400 FORMAT(' ATOMIC DISPLACEMENTS MATRIX: NA ROWS BY NOB COLUMNS ', 1 /,' OUTPUT IN COLUMN BY COLUMN FORMAT ',/,2I4) 1410 FORMAT(3D23.16) 1420 FORMAT('1',/' CALCULATED OBSERVED DIFFERENCE', $ ' IR INTENSITY',/) 1430 FORMAT(2X,F15.2,35X,F15.6) 1440 FORMAT(2X,3F15.2,5X,F15.6) END SUBROUTINE MATOUT(A,N,IR,IC) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION A(N,N) K=-11 10 K=K+12 L=MIN0(K+11,IC) WRITE(16,1000)(J,J=K,L) DO 20 I=1,IR 20 WRITE(16,1010) I,(A(I,J),J=K,L) IF(L.LT.IC) GO TO 10 RETURN 1000 FORMAT(//3X,12I10) 1010 FORMAT('0',I4,2X,12F10.6) END SUBROUTINE INVERT(A,MA,DETA,IER,IROW,ICOL) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) DIMENSION A(MX,MX),IR(MX),IC(MX) IER=0 DO 10 I=1,MA IR(I)=0 10 IC(I)=0 DETA=1.D0 DO 90 M=1,MA TEMP=0.D0 DO 30 K=1,MA IF(IR(K).NE.0)GO TO 30 DO 20 L=1,MA IF(IC(L).NE.0)GO TO 20 T=DABS(A(K,L)) IF(T.LT.TEMP)GO TO 20 I=K J=L TEMP=T 20 CONTINUE 30 CONTINUE IF(M.EQ.1)THRES=TEMP*1.D-14 PIV=A(I,J) DETA=DETA*PIV IR(I)=J IC(J)=I IF(TEMP.GT.THRES)GO TO 40 IF(IER.NE.0)GO TO 90 IER=-1 IROW=J ICOL=I GO TO 90 40 PIV=1.D0/PIV DO 50 K=1,MA 50 A(I,K)=A(I,K)*PIV A(I,J)=PIV DO 70 K=1,MA IF(K.EQ.I)GO TO 70 PIV1=A(K,J) DO 60 L=1,MA 60 A(K,L)=A(K,L)-PIV1*A(I,L) A(K,J)=PIV1 70 CONTINUE DO 80 K=1,MA 80 A(K,J)=-PIV*A(K,J) A(I,J)=PIV 90 CONTINUE DO 120 I=1,MA K=IC(I) M=IR(I) IF(K.EQ.I)GO TO 120 DETA=-DETA DO 100 L=1,MA TEMP=A(K,L) A(K,L)=A(I,L) 100 A(I,L)=TEMP DO 110 L=1,MA TEMP=A(L,M) A(L,M)=A(L,I) 110 A(L,I)=TEMP IC(M)=K IR(K)=M 120 CONTINUE RETURN END SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) DIMENSION A(MX,MX),Z(MX,MX),D(MX),E(MX) COMMON/DIAG/E COMMON/INDEX/IND(MX) 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) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 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) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) DIMENSION Z(MX,MX),D(MX),E(MX) REAL*8 MACHEP COMMON/DIAG/E COMMON/INDEX/IND(MX) 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 IY=IND(I) IND(I)=IND(K) IND(K)=IY 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 SMATOUT(S,AS,NINT,NAT3,IRANGE,JRANGE,OK1) C dX = S dQ S(nat3xnint) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) DIMENSION S(MX,MX),AS(3,MX) CHARACTER*1 OK1 logical lex COMMON /INDEX/IND(MX) IRANGE=0 WRITE(*,*)' RANGE OF MODES FOR THE INP FILE: ' IF(OK1.NE.'N')READ(*,*) IRANGE,JRANGE NI=NINT IF(IRANGE.NE.0)NI=JRANGE-IRANGE+1 OPEN(34,FILE='F.INP') WRITE(34,10)NI,NAT3,NAT3/3 10 FORMAT(3I4) inquire(file='FILE.XYZ',exist=lex) if(lex)then OPEN(35,FILE='FILE.XYZ') READ(35,*) DO 3 I=1,NAT3/3 READ(35,*)IAT,X,Y,Z 3 WRITE(34,11)IAT,X,Y,Z 11 FORMAT(I4,3F12.6) else OPEN(35,FILE='FILE.X',status='old') READ(35,*) READ(35,*) DO 31 I=1,NAT3/3 READ(35,*)IAT,X,Y,Z 31 WRITE(34,11)IAT,X,Y,Z endif CLOSE(35) WRITE(34,14) 14 FORMAT(1X,'Atom Mode X-disp. Y-disp. Z-disp.') DO 1 I=1,NAT3/3 IU=3*(I-1) DO 1 J=IRANGE,JRANGE 1 WRITE(34,15)I,J,S(IU+1,J ),S(IU+2,J ),S(IU+3,J ) 15 FORMAT(2I4,3F10.6) WRITE(34,11)NI WRITE(34,13)(AS(1,I),I=IRANGE,JRANGE) 13 format(6F11.3) CLOSE(34) RETURN END SUBROUTINE IMPROVE2(NIT,S0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) COMMON /ITERI/NFC,NOB,NOAT,NOSC COMMON /ITERA/AMULT(MX),NCORRESP(MX),X((MX*(MX+1))/2),G(MX,MX), 1Q(MX),ZNU(MX,5),Z(MX,MX),P(MX),A(MX,MX),AA(MX,MX), 2Y(MX,MX),FF(MX,MX),SCFAC(MX),DEL COMMON/DIAG/B(MX),AKSI(MX),GR(MX) STEP=0.00001 ANORM0=0.25 ALIMIT=0.0005 SLIMIT=0.01 IIT=0 S0=SUM(0) 1 IIT=IIT+1 IF (IIT.GT.NIT)RETURN IF (MOD(IIT,10).EQ.0)THEN WRITE(16,*)' ITERATION ',IIT WRITE(16,*)' THE SCALE FACTORS:' WRITE(16,1616)(SCFAC(I),I=1,NOSC) 1616 FORMAT(6G12.6) ENDIF DO 2 I=1,NOSC X0=SCFAC(I) SCFAC(I)=X0-STEP SM=SUM(0) SCFAC(I)=X0+STEP SP=SUM(0) SCFAC(I)=X0 B(I)=(SP+SM-2.0*S0)/STEP/STEP GR(I)=(SP-SM)/2.0/STEP 2 B(I)=1.0/B(I) ANORM=0.0 DO 3 I=1,NOSC AKSI(I)=0.0 IF(B(I).GT.0.0)THEN AKSI(I)=B(I)*GR(I) ELSE AKSI(I)=GR(I)/ABS(GR(I))*DEL ENDIF 3 ANORM=ANORM+AKSI(I)*AKSI(I) ANORM=SQRT(ANORM) IF(ANORM.GT.ANORM0)THEN DO 5 I=1,NOSC 5 AKSI(I)=AKSI(I)*ANORM0/ANORM ANORM=ANORM0 ENDIF DO 4 I=1,NOSC 4 SCFAC(I)=SCFAC(I)-AKSI(I) 8 SUMNEW=SUM(0) WRITE(16,1617)ANORM,SUMNEW,S0 1617 FORMAT(' NORM: ',F12.6,' SUM: ',F12.6,' (OLD: ',F12.6,' )') IF(SUMNEW.GT.S0)THEN WRITE(16,*)' ITERATIVE BACK-STEP' DO 6 I=1,NOSC SCFAC(I)=SCFAC(I)+AKSI(I)/2.0 6 AKSI(I)=AKSI(I)/2.0 ANORM=ANORM/2.0 IF(ANORM.LT.STEP)THEN WRITE(16,*)' IMPOSSIBLE TO CONTINUE THE ITERATION' RETURN ENDIF GOTO 8 ENDIF IF(ANORM.LT.ALIMIT.OR.S0-SUMNEW.LT.SLIMIT)THEN WRITE(16,*)' ITERATION FINISHED' RETURN ENDIF S0=SUMNEW GOTO 1 END SUBROUTINE IMPROVE(NIT,S0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) COMMON /ITERI/NFC,NOB,NOAT,NOSC COMMON /ITERA/AMULT(MX),NCORRESP(MX),X((MX*(MX+1))/2),G(MX,MX), 1Q(MX),ZNU(MX,5),Z(MX,MX),P(MX),A(MX,MX),AA(MX,MX), 2Y(MX,MX),FF(MX,MX),SCFAC(MX),DEL IIT=0 1 IIT=IIT+1 IF (IIT.GT.NIT)RETURN IF (MOD(IIT,10).EQ.0)THEN WRITE(16,*)' ITERATION ',IIT WRITE(16,*)' THE SCALE FACTORS:' WRITE(16,1616)(SCFAC(I),I=1,NOSC) 1616 FORMAT(6G12.6) ENDIF ICHANGE=0 DO 2 I=1,NOSC X0=SCFAC(I) XP=X0+DEL XM=X0-DEL S0=SUM(0) SCFAC(I)=XM SM=SUM(0) SCFAC(I)=XP SP=SUM(0) SCFAC(I)=X0 IF (SP.LT.S0)THEN S0=SP SCFAC(I)=XP ICHANGE=1 ENDIF IF (SM.LT.S0)THEN S0=SM SCFAC(I)=XM ICHANGE=1 ENDIF 2 CONTINUE IF (ICHANGE.EQ.0)DEL=DEL/2.0D0 WRITE(*,6002)DEL 6002 FORMAT(' Scale factors varied +/- ',F12.8) IF (DEL.GT.0.0001D0)GOTO 1 RETURN END FUNCTION SUM(IDUMM) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) COMMON /ITERA/AMULT(MX),NCORRESP(MX),X((MX*(MX+1))/2),G(MX,MX), 1Q(MX),ZNU(MX,5),Z(MX,MX),P(MX),A(MX,MX),AA(MX,MX), 2Y(MX,MX),FF(MX,MX),SCFAC(MX),DEL COMMON /ITERI/NFC,NOB,NOAT,NOSC COMMON /INDEX/IND(MX) EQUIVALENCE(N,NFC) IERR=IDUMM KKK=1 DO 131 I=1,NOB AMULT(I)=SCFAC(KKK) 131 IF(I.EQ.NCORRESP(KKK)) KKK=KKK+1 K=0 DO 130 I=1,NOB DO 130 J=1,I K=K+1 130 FF(I,J)=X(K)*DSQRT(AMULT(I)*AMULT(J)) DO 1 I=1,180 1 IND(I)=I CALL TRED12(NOB,FF,Z,P,2,IERR) IF(IERR.EQ.0)GO TO 200 WRITE(*,*)' DIAGONALIZATION ERROR' STOP 200 CONTINUE DO 210 I=1,NOB T=P(I) IF(T.LT.1.D-8)T=0.D0 T=DSQRT(T) IF(DABS(T).GT.1.0D-10)P(I)=1.D0/T DO 210 J=1,NOB 210 A(I,J)=T*Z(J,I) DO 270 I=1,NOB DO 270 J=1,NOB T=0.D0 DO 260 K=1,NOB 260 T=T+G(I,K)*A(J,K) 270 AA(I,J)=T DO 290 I=1,NOB DO 290 J=1,I T=0.D0 DO 280 K=1,NOB 280 T=T+A(I,K)*AA(K,J) 290 Y(I,J)=T DO 2 I=1,NOB 2 IND(I)=NOB-I+1 CALL TRED12(NOB,Y,Y,Q,-1,IERR) IF(IERR.NE.0)STOP FWAVE=0 DO 310 I=1,NOB T=DABS(Q(I)) T=DSQRT(T)*1302.828D0 TP=ZNU(I,1) QQ=T-TP IF (TP.GT.0.0001)THEN FWAVE=FWAVE+QQ*QQ ENDIF 310 CONTINUE NF=NOB FWAVE=DSQRT(FWAVE/DFLOAT(NF)) SUM=FWAVE RETURN END SUBROUTINE FPC (V,CHAR,S,NOAT,NOB) C PETR BOUR: LAST UPDATE 4-28-92 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=436) logical lex DIMENSION V(MX), S(MX,MX), X(3,MX/3), 1 CHAR(MX/3), DS(3), SCRT(3) DCOR = 388891.4D0 RCOR = 12221.7D0 inquire(file='FILE.XYZ',exist=lex) if(lex)then OPEN(35,FILE='FILE.XYZ',STATUS='OLD') READ(35,*) DO 51 J=1,NOAT 51 READ(35,*)IAT,(X(I,J),I=1,3) else OPEN(35,FILE='FILE.X',STATUS='OLD') READ(35,*) READ(35,*) DO 511 J=1,NOAT 511 READ(35,*)IAT,(X(I,J),I=1,3) endif CLOSE(35) NA=3*NOAT NVIBS=NA-6 OPEN(13,FILE='FPC.TAB') WRITE(13,50) 50 FORMAT(3X,'FREQ',8X,'DIP. STRN.*E39',6X,'ROT. STRN.*E45'/ 1 3X,'CM-1',2X,2(6X,' (FR CM)**2 '),/,'---------------') DO 42 I=1,NVIBS DO 784 IKY=1,3 DS(IKY)=0.0D0 SCRT(IKY)=0.0D0 DO 784 JR=1,NOAT DS(IKY)=DS(IKY)+CHAR(JR)*S(3*(JR-1)+IKY,I) DO 784 IB=1,3 DO 784 IG=1,3 784 SCRT(IKY)=SCRT(IKY)+EPS(IKY,IB,IG)*X(IB,JR)*S(3*(JR-1)+IG,I) 1 *CHAR(JR) DT =DS(1)*DS (1)+DS(2)*DS (2)+DS(3)*DS (3) ROT=DS(1)*SCRT(1)+DS(2)*SCRT(2)+DS(3)*SCRT(3) ROTC =+ROT*RCOR DIPOL=0.0D0 IF (V(I).GT.0) DIPOL = DT*DCOR/V(I) 42 WRITE(13,85) I,V(I),DIPOL,ROTC,DT,ROT 85 FORMAT(I4,F8.2,' ',F13.6,' ',f13.6,' ',2G14.5) WRITE(13,86) 86 FORMAT('---------------------------------------------') CLOSE(13) 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