PROGRAM FTRY IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) real*8,allocatable::BB(:,:),SCFAC(:),CHAR(:),ZIM(:),ZNU(:), 1FO(:,:),FF(:,:),AMULT(:),ZM(:),ZL(:,:),QI(:),AMODE(:),V(:,:), 2BD(:,:),PE(:),AS(:,:),ZINTY(:),G(:,:),Q(:),Z(:,:),P(:), 3A(:,:),AA(:,:),Y(:,:),AKSI(:),GR(:),R(:,:),B(:),YY(:,:), 4SS(:,:),ZINT0(:),AREM(:,:),Q0(:) integer*4,allocatable::NCORRESP(:),MODE(:),IPE(:) CHARACTER*4 CODE CHARACTER*1 OK,OK1 CHARACTER*120 s120 character*2,allocatable::a1(:),a2(:),a3(:) character*1,allocatable::bt(:),ty(:) WRITE(*,*) WRITE(*,*)' FORCE FIELD REFINEMENT PROGRAM' WRITE(*,*)' Fortran 95 version' WRITE(*,*)' Petr Bour, Prague 1992-2015' WRITE(*,*) ia1=iargc() OPEN(16,FILE='FTRY.OUT') WRITE(16,*)' FTRY OUTPUT FILE' WRITE(16,*)' ----------------' DEL=0.02D0 IRMS=0 IERR=0 ITRY=0 c NOB=1 NA=1 allocate(BB(NOB,NA)) call rb(0,NOB,NA,BB) deallocate(BB) allocate(BB(NOB,NA)) call rb(1,NOB,NA,BB) NOAT=NA/3 write(6,*)'B.MAT read,',NOAT,' atoms,',NOB,' coordinates' allocate(CHAR(NOAT),ZIM(NOAT),ZNU(NOB),Q0(NOB), 1FO(NOB,NOB),FF(NOB,NOB), 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),ZINT0(NOB),AREM(NOB,NOB)) write(6,*)'Matrices allocated' BD=0.0d0 YY=0.0d0 OPEN(15,FILE='FTRY.INP') 10 READ(15,1070,END=690)CODE,(MODMIN,i=1,7),MODMIN,MODMAX 1070 FORMAT(A4,9i4) READ(15,*) NOSC allocate(SCFAC(NOSC),NCORRESP(NOSC),GR(NOSC),AKSI(NOSC), 1B(NOSC)) READ(15,1040) (SCFAC(I),I=1,NOSC) 1040 FORMAT(6F12.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 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) ZIM(I)=1.D0/ZIM(I) ZM(3*I )=ZIM(I) ZM(3*I-1)=ZIM(I) 21 ZM(3*I-2)=ZIM(I) WRITE(*,877)am 877 FORMAT(' TOTAL MASS',F15.2) c call rdfc(FO,NOB) c ITRY=0 111 ITRY=ITRY+1 WRITE(16,*)' TRIAL ',ITRY WRITE(16,*) ' THE SCALE FACTORS' WRITE(16,*)' COORD SF# SF FROM TO' call asm(AMULT,SCFAC,NOB,NOSC,NCORRESP) WRITE(16,1020) NOSC WRITE(16,1040) (SCFAC(I),I=1,NOSC) WRITE(16,1020) (NCORRESP(I),I=1,NOSC) call scaf(NOB,FF,FO,AMULT) call wrff('SCALED.FC',FF,NOB) 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) AREM=A NF=0 FWAVE=0.D0 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) YY(I,J)=T 290 YY(J,I)=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 Q0=Q 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)FWAVE,RMSW,AK 5551 FORMAT(1X,'************************************************',/, 1 1X,'>>>>>> ',F8.2,' << >> ',F8.2,' <<<<<< ',F8.6,/) OK='N' if(OK1.ne.'N')then if(ia1.lt.1)then WRITE(*,6667) 6667 FORMAT(1X,' NEW SCALEFACTORS (Y/N) ?') READ(*,'(A)')OK else call getarg(1,OK) endif endif 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 OK='N' if(OK1.ne.'N')then if(ia1.lt.2)then WRITE(*,*)' AUTOMATIC REFINEMENT (Y/N) ?' READ(*,'(A)')OK else call getarg(2,OK) endif endif 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, 1 AMULT,NCORRESP,FO,G,Q,Z,P,A,AA,Y,FF,SCFAC,DEL,YY) IF(ITYPE.EQ.2)CALL IMPROVE2(NA,NIT,S0,ZNU,NOB,NOSC, 1 AMULT,NCORRESP,FO,G,Q,Z,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 DO 430 J=1,NOB A(I,J)=0.D0 DO 420 K=1,NOB 420 A(I,J)=A(I,J)+ZL(K,I)*FF(K,J) 430 A(I,J)=A(I,J)*QI(I) 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 1 i=1,NOB write(17,1717)i,AS(1,i) 1717 format(i6,f10.2,$) do 2 J=1,NOB 2 write(17,17171)AA(J,i) 17171 format(f10.6,$) 1 write(17,*) close(17) open(8,file='FILE.UMA') read(8,*) read(8,*) read(8,*)nat,nc allocate(ty(nc),A1(nc),A2(nc),A3(nc),bt(nc)) ty='?' bt='?' A1=' X' A2=' X' A3=' X' do 7 J=1,nat 7 read(8,*) do 8 J=1,nc read(8,800)s120 800 format(a120) read(s120,*)itype if(itype.eq.1)then ty(J)='v' A1(J)=s120(65:66) A2(J)=s120(74:75) bt(J)=s120(72:72) else if(itype.eq.2)then ty(J)='d' A1(J)=s120(63:64) A2(J)=s120(72:73) A3(J)=s120(81:82) else if(itype.eq.4)then ty(J)='t' A1(J)=s120(66:67) A2(J)=s120(77:78) read(8,*) read(8,*) else goto 444 endif endif endif 8 continue 444 close(8) OPEN(20,FILE='POT.TXT') OPEN(19,FILE='POT.EN') WRITE(19,1921) WRITE(20,1921) 1921 format(' POTENTIAL ENERGY DISTRIBUTION [in %]',/, 1 ' MODE coordinate(PED) (for PED > 5%)',/, 180(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=AA(1,J) IMAX=1 c c Second loop over coordinates - find maximum and zero after DO 452 I=2,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(PE(1).gt.99.9d0)PE(1)=99.9d0 WRITE(19,1923)NOB-J+1,AS(1,J) 1923 FORMAT(I5,F9.1,$) WRITE(20,1924)J,AS(1,J) 1924 FORMAT(I5,F9.0,' ',$) do 9 I=1,NPE K=IPE(I) c is this type already recorded?: do 11 II=1,I-1 KK=IPE(II) 11 if(ty(K).eq.ty(KK).and.A1(K).eq.A1(KK).and.A2(K).eq.A2(KK).and. 1A3(K).eq.A3(KK).and.bt(K).eq.bt(KK))goto 9 write(20,1926)ty(K) 1926 FORMAT(a1,$) if(ty(K).eq.'v')then write(20,1927) 1927 FORMAT('(',$) if(A1(K).eq.' H')then call wr1(A2(K)) call wr1(A1(K)) else call wr1(A1(K)) call wr1(' '//bt(K)) call wr1(A2(K)) endif write(20,1928) 1928 FORMAT(')',$) endif if(ty(K).eq.'d')then write(20,1927) call wr1(A1(K)) call wr1(A2(K)) call wr1(A3(K)) write(20,1928) endif if(ty(K).eq.'t')then write(20,1927) call wr1(A1(K)) call wr1(A2(K)) write(20,1928) endif if(I.lt.NPE)write(20,1929) 1929 FORMAT(', ',$) WRITE(19,1925)K,PE(i) 1925 FORMAT(I5,'(',F4.1,')',$) 9 continue WRITE(19,*) WRITE(20,*) 451 CONTINUE WRITE(19,1920) 1920 format(80(1H-)) DO 455 I=1,NOB 455 WRITE(19,1919)I,NOB+1-MODE(I),AS(1,MODE(I)),AMODE(I) 1919 FORMAT(' Coordinate',I4,' is most present in mode',I4, 1f10.2, 'cm-1 (by',F5.1,'%).') WRITE(19,1920) CLOSE(19) CLOSE(20) allocate(SS(NA,NOB)) DO 480 I=1,NA DO 480 J=1,NOB SS(I,J)=0.D0 DO 480 K=1,NOB 480 SS(I,J)=SS(I,J)+BB(K,I)*A(J,K) DO 470 I=1,NA DO 470 J=1,NOB 470 SS(I,J)=SS(I,J)*ZM(I) WRITE(*,*)' S -MATRIX COMPUTED' CALL SMATOUT(SS,AS,NOB,NA,MODMIN,MODMAX,OK1,ia1) IF(CODE.EQ.'ATOM') WRITE(16,1240) IF(CODE.EQ.'ATOM')CALL MATOUT(SS,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)*SS(JJX,I) ZMUY = ZMUY + CHAR(J)*SS(JJY,I) 484 ZMUZ = ZMUZ + CHAR(J)*SS(JJZ,I) 485 ZINTY(I) = ZMUX*ZMUX + ZMUY*ZMUY + ZMUZ*ZMUZ ZINT0 = ZINTY CALL FPC(Q,CHAR,SS,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(SS(I,J)).LE.1.0D-10) SS(I,J) = 0.0D0 488 CONTINUE WRITE(7,1410) (SS(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 ZINTY=ZINT0 BGZTY = 0.0D0 DO 492 I=1,1 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 A=AREM NF = 0 FWAVE = 0.0D0 WRITE(16,1420) ZINTY=ZINT0 DO 494 L=1,1 WRITE(16,1280) L Q=Q0 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') 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) integer*4 ,allocatable::IC(:),IR(:) 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,NOB,NAT3,IRANGE,JRANGE,OK1,ia1) C dX = S dQ S(nat3xnint) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION S(NAT3,NOB),AS(3,NOB) CHARACTER*1 OK1 CHARACTER*10 s10 IRANGE=1 JRANGE=NOB IF(OK1.NE.'N')then if(ia1.lt.4)then WRITE(*,*)' RANGE OF MODES FOR THE INP FILE: ' READ(*,*) IRANGE,JRANGE else call getarg(3,s10) read(s10,*)IRANGE call getarg(4,s10) read(s10,*)JRANGE endif endif if(IRANGE.lt.1.or.IRANGE.gt.NOB) 1call report('Invalid lower limit') if(JRANGE.lt.1.or.JRANGE.gt.NOB) 1call report('Invalid upper limit') IF(IRANGE.NE.0)NOB=JRANGE-IRANGE+1 OPEN(34,FILE='F.INP') WRITE(34,10)NOB,NAT3,NAT3/3 10 FORMAT(3I6) 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(I6,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+k,J),k=1,3) 15 FORMAT(2I6,3F12.6) WRITE(34,11)NOB 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,FO,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(*),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),FO(NOB,NOB) STEP=0.00001d0 ANORM0=0.25d0 ALIMIT=0.0005d0 SLIMIT=0.01d0 IIT=0 S0=SSM(ZNU,AMULT,NCORRESP,FO,G,NOB,NOSC,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,FO,G,NOB,NOSC,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) SCFAC(I)=X0+STEP SP=SSM(ZNU,AMULT,NCORRESP,FO,G,NOB,NOSC,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,FO,G,NOB,NOSC,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) WRITE(*,6001)(SCFAC(i),i=1,NOSC) 6001 FORMAT(10f8.4) WRITE(6,1617)ANORM,SUMNEW,S0 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,FO,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(*),G(NOB,NOB),Q(*), 1Z(NOB,NOB),P(*),A(NOB,NOB),AA(NA,NOB),Y(NOB,NOB),FF(NOB,NOB), 1SCFAC(*),YY(NOB,NOB),FO(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,FO,G,NOB,NOSC,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) SCFAC(I)=XM SM=SSM(ZNU,AMULT,NCORRESP,FO,G,NOB,NOSC,Q, 1Z,P,A,AA,Y,FF,SCFAC,YY) SCFAC(I)=XP SP=SSM(ZNU,AMULT,NCORRESP,FO,G,NOB,NOSC,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(*,6001)(SCFAC(i),i=1,NOSC) 6001 FORMAT(10f8.4) WRITE(*,6002)S0,DEL 6002 FORMAT('ERR,DEL: ',f12.4,F9.5) IF (DEL.GT.0.0001D0)GOTO 1 RETURN END FUNCTION SSM(ZNU,AMULT,NCORRESP,FO,G,NOB,NOSC,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(*),G(NOB,NOB),Q(*), 1Z(NOB,NOB),P(*),A(NOB,NOB),AA(NOB,NOB),Y(NOB,NOB),FF(NOB,NOB), 1SCFAC(*),YY(NOB,NOB),FO(NOB,NOB) IERR=0 call asm(AMULT,SCFAC,NOB,NOSC,NCORRESP) call scaf(NOB,FF,FO,AMULT) 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) DIMENSION V(*), S(3*NOAT,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 integer*4 I,J,K 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 subroutine wr1(A) character*2 A if(A(1:1).eq.' ')then write(20,201)A(2:2) 201 FORMAT(a1,$) else write(20,202)A 202 FORMAT(a2,$) endif return end subroutine rb(ic,NOB0,NA,B) implicit none integer*4 ic,NA,i,NOB0,nn,ii,ia,ix real*8 b(NOB0,NA) open(9,file='B.MAT') read(9,*)NOB0 read(9,*)na if(ic.eq.0)goto 99 b=0.0d0 read(9,*) do 13 i=1,NOB0 read(9,*) read(9,*)nn do 13 ii=1,nn 13 read(9,*)ia,ix,b(i,ix+3*(ia-1)) 99 close(9) return end subroutine rdfc(FO,NOB) implicit none integer*4 NOB,N,N1,N3,LN,J,I real*8 FO(NOB,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 read(8,*)FO(LN,N1),(FO(LN,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) DO 31 I=1,NOB DO 31 J=I+1,NOB 31 FO(I,J)=FO(J,I) write(6,*)N,' force constants from INTY.FC' return end subroutine wrff(s,FF,NOB) implicit none integer*4 NOB,M,MI,MJ,I,MW,J real*8 FF(NOB,NOB) character*(*) s OPEN(35,FILE=s) WRITE(35,3337)(NOB*(NOB+1))/2,NOB,NOB 3337 FORMAT(' SCALED FORCE CONSTANTS',/,/, 1I6,' INTRINSIC FORCE CONSTANTS',I5,' x (',I5,' + 1) / 2') M=6 MI=-5 207 MI=MI+6 MJ=M-5 IF(M.GT.NOB)M=NOB WRITE(35,28)(I,I=MJ,M) 28 FORMAT(8X,6(I7,5X)) DO 206 I=MI,NOB MW=M IF(I.LT.MW)MW=I 206 WRITE(35,27)I,(FF(J,I),J=MJ,MW) 27 FORMAT(I5,6e13.5) IF (M.LT.NOB)THEN M=M+6 GOTO 207 ENDIF CLOSE(35) return end subroutine asm(AMULT,SCFAC,NOB,NOSC,NCORRESP) implicit none integer*4 NOB,NOSC,K,I,NCORRESP(*),ISTART real*8 AMULT(NOB),SCFAC(*) AMULT=1.0d0 do 1 K=1,NOSC if(K.eq.1)then ISTART=1 else ISTART=NCORRESP(K-1)+1 endif do 1 i=ISTART,NCORRESP(K) 1 AMULT(i)=SCFAC(K) return end subroutine report(i) character*(*) i write(6,*)i stop end subroutine scaf(NOB,FF,FO,AMULT) implicit none integer*4 NOB,I,J real*8 FF(NOB,NOB),FO(NOB,NOB),AMULT(*) DO 130 I=1,NOB DO 130 J=1,I FF(I,J)=FO(I,J)*DSQRT(AMULT(I)*AMULT(J)) 130 FF(J,I)=FF(I,J) return end