PROGRAM FTRY95 IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) real*8,allocatable::BB(:,:),SCFAC(:),CHAR(:),ZIM(:),ZNU(:), 1X(:),FF(:,:),AMULT(:),ZM(:),ZL(:,:),QI(:),AMODE(:),V(:,:), 2BD(:,:),PE(:),AS(:,:),ZINTY(:),G(:,:),Q(:),Z(:,:),P(:), 3A(:,:),AA(:,:),Y(:,:),AKSI(:),GR(:),R(:,:),B(:),YY(:,:) integer*4,allocatable::NCORRESP(:),MODE(:),IFF(:),IPE(:) CHARACTER*4 CODE CHARACTER*1 OK,OK1 WRITE(*,*) WRITE(*,*)' FORCE FIELD REFINEMENT PROGRAM' WRITE(*,*)' Fortran 95 version' WRITE(*,*)' Petr Bour, Prague 1992-2015' WRITE(*,*) OPEN(16,FILE='FTRY.OUT') WRITE(16,*)' FTRY OUTPUT FILE' WRITE(16,*)' ----------------' DEL=0.05D0 IRMS=0 IERR=0 ITRY=0 c OPEN(9,FILE='B.MAT') read(9,*)NOB read(9,*)NA NOAT=NA/3 allocate(BB(NOB,NA),CHAR(NOAT),ZIM(NOAT),ZNU(NOB), 1X((NOB*(NOB+1))/2),FF(NOB,NOB),IFF((NOB*(NOB+1))/2), 1AMULT(NOB),Z(NOB,NOB),ZM(NA),ZL(NOB,NOB),QI(NOB),MODE(NOB), 1AMODE(NOB),V(NA,NOB),BD(NA,NA),PE(NOB),IPE(NOB),AS(3,NOB), 1ZINTY(NOB),G(NOB,NOB),Q(NOB),P(NOB),A(NOB,NOB),YY(NOB,NOB), 1AA(NOB,NOB),Y(NOB,NOB),R(3,NOAT)) do 141 i=1,NOB do 141 j=1,NOB 141 YY(i,j)=0.0d0 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) JMAX=NOB*(NOB+1)/2 write(6,*)'B.MAT read,',NOAT,' atoms,',NOB,' coordinates' OPEN(15,FILE='FTRY.INP') 10 READ(15,1070,END=690)CODE,MODMIN,MODMAX READ(15,1020) NOSC allocate(SCFAC(NOSC),NCORRESP(NOSC),GR(NOSC),AKSI(NOSC), 1B(NOSC)) READ(15,1040) (SCFAC(I),I=1,NOSC) 1040 FORMAT(6G12.6) READ(15,1020) (NCORRESP(I),I=1,NOSC) READ(15,1040) (CHAR(I),I=1,NOAT) READ(15,*) READ(15,1040)(ZIM(I),I=1,NOAT) READ(15,*) READ(15,1040)(ZNU(I),I=1,NOB) CLOSE(15) WRITE(6,*)'FTRY.INP read' WRITE(*,*)NOSC,' scale factors' if(NOSC.eq.0)then WRITE(*,*)' .... (SF = 1 for all)' deallocate(SCFAC,NCORRESP) allocate(SCFAC(1),NCORRESP(1)) SCFAC(1)=1.0d0 NCORRESP(1)=NOB endif write(6,1070)CODE,MODMIN,MODMAX c CODE IPB,IPF,IPG,IPL,IPP,IPA,IPS c CODE 0 1 0 1 0 1 0 1070 FORMAT(A4 , 4x, 4x, 4x, 4x, 4x, 4x, 4x,2I4) OK1='Y' IPS=0 IF(CODE.EQ.'AUTO')THEN CODE='INTY' OK1='N' ENDIF CHATOT=0.0D0 DO 77 I=1,NOAT 77 CHATOT=CHATOT+CHAR(I) WRITE(*,777)CHATOT 777 FORMAT(' TOTAL CHARGE',F10.2) am=0.0d0 DO 21 I=1,NOAT am=am+ZIM(I) 21 ZIM(I)=1.D0/ZIM(I) WRITE(*,877)am 877 FORMAT(' TOTAL MASS',F15.2) c 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 read(8,*)X((LN*(LN-1))/2+N1),(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' c 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 ITRY=0 111 ITRY=ITRY+1 WRITE(16,*)' TRIAL ',ITRY WRITE(16,*) ' THE SCALE FACTORS' WRITE(16,*)' COORD SF# SF FROM TO' KKK=1 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)) IF(L.LT.0)T=-T FF(I,J)=T FF(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 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,1 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 ZM(3*I )=ZIM(I) ZM(3*I-1)=ZIM(I) 220 ZM(3*I-2)=ZIM(I) 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 YY(I,J)=T WRITE(16,1210) IF(CODE.NE.'INTY') WRITE(16,1280) L CALL TRED12(NOB,YY,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) NF=NF+1 QQ=T-TP IF (TP.GT.0.0001d0)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 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(NA,NIT,S0,ZNU,NOB,NOSC, 1AMULT,NCORRESP,X,G,Q, 1Z,P,A,AA,Y,FF,SCFAC,DEL,YY) IF(ITYPE.EQ.2)CALL IMPROVE2(NA,NIT,S0,ZNU,NOB,NOSC, 1AMULT,NCORRESP,X,G,Q, 1Z,P,A,AA,Y,FF,SCFAC,AKSI,GR,DEL,B,YY) 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 c if(I.eq.NOB.and.J.eq.NOB)write(6,*)K,Z(I,K),Y(K,J) 350 T=T+Z(I,K)*Y(K,J) A(I,J)=T 360 ZL(I,J)=T c write(6,*)A(1,1),A(NOB,NOB) c stop 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 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 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,NOB,NOB,NOB) OPEN(17,FILE='POT_all.LST') do i=1,NOB write(17,1717)i,AS(1,i) 1717 format(i6,f10.2,$) do J=1,NOB write(17,17171)AA(J,i) 17171 format(f10.6,$) enddo write(17,*) enddo close(17) 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) 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,NOB,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,R) IF(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,1 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,1 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)).LT.1.0D-10) GO TO 493 QQ = Q(J) - ZNU(J) NF = NF + 1 FWAVE = FWAVE + QQ*QQ WRITE(16,1440) Q(J),ZNU(J),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(CODE.EQ.'INTY') STOP IF (CODE.EQ.'ATOM')READ(15,1040)TK IF(IERR.NE.0)GO TO 10 IF(IRMS.lt.0)then goto 520 else if(IRMS.eq.0)then goto 680 else goto 500 endif endif c 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) 1160 FORMAT(/'0B MATRIX USED AS BD MATRIX, FOR RMS AMPLITUDES OF ', $ 'THE INTERNAL COORDINATES') 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) 1300 FORMAT(/'0UNIT MATRIX USED AS BD MATRIX, FOR RMS CARTESIAN ', $ 'DISPLACEMENTS') 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) DO 590 J=1,NA T=ZM(J) DO 590 I=1,M 590 BD(I,J)=BD(I,J)*T CALL MATOUT(BD,NOB,M,NA) 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 1020 FORMAT(20I4) 1060 FORMAT('0INVERSION OF L MATRIX FAILED: CAUSE IS LINEAR', $ ' DEPENDENCE OF THE EIGENVECTORS') 1140 FORMAT('0* ELEMENT (',I5,',',I5,') OF F MATRIX REQUIRES NON-', $ 'EXISTENT FORCE CONSTANT NO',I5//) 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'//) 1280 FORMAT('0MOLECULE',I4/) 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) DIMENSION A(MA,MA) real*8,allocatable::IR(:),IC(:) allocate(IR(MA),IC(MA)) 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 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) DIMENSION S(NINT,NINT),AS(3,NINT) CHARACTER*1 OK1 logical lex IRANGE=1 JRANGE=NAT3-6 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) OPEN(35,FILE='FILE.X') READ(35,*) READ(35,*)NAT DO 3 I=1,NAT READ(35,*)IAT,X,Y,Z 3 WRITE(34,11)IAT,X,Y,Z 11 FORMAT(I4,3F12.6) 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(NA,NIT,S0,ZNU,NOB,NOSC, 1AMULT,NCORRESP,X,G,Q, 1Z,P,A,AA,Y,FF,SCFAC,AKSI,GR,DEL,B,YY) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION ZNU(*),AMULT(*),NCORRESP(*),X(*),G(NOB,NOB),Q(*), 1Z(NOB,NOB),P(*),A(NOB,NOB),AA(NA,NOB),Y(NOB,NOB),FF(NOB,NOB), 1SCFAC(*),AKSI(*),GR(*),B(*),YY(NOB,NOB) STEP=0.00001d0 ANORM0=0.25d0 ALIMIT=0.0005d0 SLIMIT=0.01d0 IIT=0 S0=SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) 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=SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) SCFAC(I)=X0+STEP SP=SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) SCFAC(I)=X0 B(I)=(SP+SM-2.0d0*S0)/STEP/STEP GR(I)=(SP-SM)/2.0d0/STEP 2 B(I)=1.0d0/B(I) ANORM=0.0d0 DO 3 I=1,NOSC AKSI(I)=0.0d0 IF(B(I).GT.0.0d0)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=SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) 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.0d0 6 AKSI(I)=AKSI(I)/2.0d0 ANORM=ANORM/2.0d0 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(NA,NIT,S0,ZNU,NOB,NOSC, 1AMULT,NCORRESP,X,G,Q, 1Z,P,A,AA,Y,FF,SCFAC,DEL,YY) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION ZNU(*),AMULT(*),NCORRESP(*),X(*),G(NOB,NOB),Q(*), 1Z(NOB,NOB),P(*),A(NOB,NOB),AA(NA,NOB),Y(NOB,NOB),FF(NOB,NOB), 1SCFAC(*),YY(NOB,NOB) 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=SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) SCFAC(I)=XM SM=SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) SCFAC(I)=XP SP=SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) 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 SSM(ZNU,AMULT,NCORRESP,X,G,NOB,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION ZNU(*),AMULT(*),NCORRESP(*),X(*),G(NOB,NOB),Q(*), 1Z(NOB,NOB),P(*),A(NOB,NOB),AA(NOB,NOB),Y(NOB,NOB),FF(NOB,NOB), 1SCFAC(*),YY(NOB,NOB) IERR=0 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)) 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 CALL TRED12(NOB,Y,YY,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) QQ=T-TP IF (TP.GT.0.0001d0)THEN FWAVE=FWAVE+QQ*QQ ENDIF 310 CONTINUE NF=NOB FWAVE=DSQRT(FWAVE/DFLOAT(NF)) SSM=FWAVE RETURN END SUBROUTINE FPC (V,CHAR,S,NOAT,NOB,X) C PETR BOUR: LAST UPDATE 4-28-92 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) logical lex DIMENSION V(*), S(NOB,NOB), X(3,NOAT), 1 CHAR(*), DS(3), SCRT(3) DCOR = 388891.4D0 RCOR = 12221.7D0 OPEN(35,FILE='FILE.X',STATUS='OLD') READ(35,*) READ(35,*) DO 511 J=1,NOAT 511 READ(35,*)X(1,J),(X(I,J),I=1,3) 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 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 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