PROGRAM FTRY C Harmonic force field diagonalization, c works in cartesian coordinates, for one molecule only IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100,MX=MX3/3) CHARACTER*1 OK CHARACTER*4 JOB LOGICAL LBMAT DIMENSION ZIM(MX),AS(3,MX3),CHAR(MX),FINT(MX3,MX3) COMMON/ITERA/AMULT(MX3),SCFAC(MX),ZM(MX3),FCAR(MX3,MX3), 1FNEW(MX3,MX3),SMAT(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX),NA,NM(MX) DATA IERR/0/ C wtol=0.01d0 c1=1302.828d0 nlinepar=iargc() ili=0 OPEN(16,FILE='FTRY.OUT') WRITE(*,60000) WRITE(16,60000) 60000 FORMAT(' FORCE FIELD REFINEMENT PROGRAM',/, 1 ' PC version - in cartesian coordinates',/, 2 ' Scaling of masses!',/,/, 3 ' Petr Bour, Prague 2-95',/,/) C OPEN(17,FILE='FILE.UMA') READ(17,*) READ(17,*) READ(17,*)NOAT NA=3*NOAT CLOSE(17) C C Read-in cartesian FF OPEN(20,FILE='FILE.FC',STATUS='OLD') CALL READFF(MX3,NA,FCAR) CLOSE(20) C C Read-in scalefactors: OPEN(15,FILE='FTRY.INP') READ(15,1500)JOB 1500 FORMAT(A4) READ(15,*)NOSC DO 222 I=1,NOSC 222 READ(15,*) SCFAC(I),NM(I),(NC(I,J),J=1,NM(I)) 1040 FORMAT(6G12.6) WRITE(*,*)' Scale factors read in...' C C Read-in atomic charges for FPC: c READ(15,1040) (CHAR(I),I=1,NOAT) READ(15,*) (CHAR(I),I=1,NOAT) CHATOT=0.0d0 DO 77 I=1,NOAT 77 CHATOT=CHATOT+CHAR(I) WRITE(*,777)CHATOT,NOAT 777 FORMAT(' TOTAL CHARGE',F20.6,/I6,' atoms') C C Read-in atomic masses: READ(15,*)NMOL IF(NMOL.NE.1)THEN WRITE(16,*)' Only one molecule allowed here !' CLOSE(16) CLOSE(15) STOP ENDIF READ(15,1040)(ZIM(I),I=1,NOAT) WRITE(16,1040)(ZIM(I),I=1,NOAT) BM=0.0d0 DO 102 IA=1,NOAT 102 BM=BM+ZIM(IA) WRITE(*,778)BM 778 FORMAT(' TOTAL MASS ',F20.6) DO 220 I=1,NOAT DO 220 J=1,3 220 ZM(3*I-3+J)=1/SQRT(ZIM(I)) C C Read-in experimental frequencies: READ(15,*)NNU IF(NNU.EQ.0)NNU=NA-6 READ(15,1040)(ZNU(I),I=1,NNU) DO 211 I=NNU+1,NA 211 ZNU(I)=0.0d0 CLOSE(15) WRITE(*,*)' FTRY.INP READ IN' C C WRITE(*,*)'Project out zero vibrations (y/n)?' OK='Y' if(nlinepar.gt.ili)then ili=ili+1 call getarg(ili,OK) write(6,*)OK else IF(JOB.NE.'AUTO')READ(*,'(A)')OK endif IF(OK.NE.'n'.and.OK.NE.'N')CALL PROJECT(FCAR,MX3,NA) C C Generate the scale constants at the start of iteration loop ITRY=0 111 ITRY=ITRY+1 WRITE(16,*)' TRIAL ',ITRY WRITE(16,*)NOSC, ' THE SCALE FACTORS' WRITE(16,*)' SF ATOMS' DO 200 I=1,NOSC 200 WRITE(16,10400) SCFAC(I),(NC(I,J),J=1,NM(I)) 10400 FORMAT(F9.6,100I3) C CALL FAST(NA,FCAR,FNEW,ZM) IF(NOSC.GT.0)THEN DO 201 I=1,NA 201 AMULT(I)=1.0d0 DO 202 I=1,NOSC SF=1.0d0/SQRT(SCFAC(I)) DO 202 J=1,NM(I) DO 202 K=1,3 202 AMULT(3*(NC(I,J)-1)+K)=SF CALL FAST(NA,FNEW,FINT,AMULT) DO 203 I=1,NA DO 203 J=1,NA 203 FNEW(I,J)=FINT(I,J) ENDIF WRITE(*,*)' New cartesian force field done ' WRITE(16,*)' TRIAL ',ITRY IMOD=2 CALL TRED12(NA,FNEW,SMAT,Q,IMOD,IERR) IF(IERR.NE.0)STOP WRITE(*,*)' Diagonalization OK' c RMSW=0.0D0 DO 310 I=1,NA C Make imaginary frequencies negative: T=Q(I) SIGN=1.0d0 IF(T.LT.0.0)SIGN=-1.0d0 T=SIGN*SQRT(ABS(T))*c1 Q(I)=T TP=ZNU(I) DEL=T-TP if(ABS(TP).le.wtol)DEL=0.0d0 RMSW=RMSW+DEL**2 WRITE(16,1230) I, T,TP,DEL 1230 FORMAT(I6,3F15.2) AS(1,I)=T AS(2,I)=TP AS(3,I)=T-TP 310 CONTINUE RMSW=DSQRT(RMSW/DFLOAT(NA)) WRITE(16,1370)RMSW 1370 format(' RMPS DEVIATION FROM NON-ZERO FREQUENCIES:',f10.4) SO1=0.0D0 SO2=0.0D0 DO 5151 J=1,NA SO1=SO1+AS(1,J)*AS(2,J) 5151 SO2=SO2+AS(1,J)*AS(1,J) AK=SO1/SO2 C WRITE(*,5555) 5555 FORMAT(1X,'************ NEW FREQUENCIES ******************') ILINE=0 DO 311 J=NA,1,-1 ILINE=ILINE+1 IF (ILINE.EQ.20)THEN IF(JOB.NE.'AUTO')THEN WRITE(*,*)' INPUT ANY CHARACTER TO CONTINUE ...' c READ(*,'(A1)')OK ENDIF ILINE=0 ENDIF a3=0.0 if(AS(2,NA-J+1).gt.wtol)a3=AS(3,NA-J+1) 311 WRITE(*,6666)J,NA-J+1,(AS(I,NA-J+1),I=1,2),a3 6666 FORMAT(1X, 1 I5,' (',I4,') ...... ',F8.2,' -',F8.2,'/',F8.2) WRITE(*,5551) 5551 FORMAT(1X,'************************************************') c c find if any experimental frequencies: ie=0 do 1 iq=1,NA 1 if(AS(2,iq).gt.wtol)ie=ie+1 if(ie.gt.0)then WRITE(6,6667)ie,RMSW,AK 6667 FORMAT(i6,'experimental frequencies found',/, 1 ' >>>>><< >> ',F8.2,' <<<<<< ',F8.6,/, 1 ' NEW SCALEFACTORS (Y/N/E) ?') OK='N' if(nlinepar.gt.ili)then ili=ili+1 call getarg(ili,OK) write(6,*)OK else IF(JOB.NE.'AUTO') READ(*,'(A)')OK 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 IF(OK.EQ.'E'.or.ok.eq.'e')then write(6,*)'Direct FF scaling to experiment, interpolate (Y/N)?' if(nlinepar.gt.ili)then ili=ili+1 call getarg(ili,OK) write(6,*)OK else read(5,'(A)')OK endif if(OK.eq.'y'.or.OK.eq.'Y')then do 2 iq=1,NA if(AS(2,iq).lt.0.01d0)then i1=0 do 3 is=iq-1,1,-1 if(AS(2,is).gt.wtol)then i1=is goto 33 endif 3 continue 33 i2=0 do 4 ie=iq+1,NA if(AS(2,ie).gt.wtol)then i2=ie goto 44 endif 4 continue 44 if(i1.gt.0.and.i2.gt.0)AS(2,iq)=AS(2,is) 1 +(AS(1,iq)-AS(1,is))/(AS(1,ie)-AS(1,is))*(AS(2,ie)-AS(2,is)) if(i1.gt.0.and.i2.eq.0)AS(2,iq)=AS(2,is)+AS(1,iq)-AS(1,is) if(i1.eq.0.and.i2.gt.0)AS(2,iq)=AS(2,ie)+AS(1,iq)-AS(1,ie) if(AS(1,iq).lt.wtol)AS(2,iq)=0.0d0 write(6,6001)iq,AS(2,iq) 6001 format(i6,' set to ',f8.2) endif 2 continue endif do 6 i=1,NA do 6 j=1,NA FCAR(i,j)=0.0d0 do 6 jp=1,NA 6 FCAR(i,j)=FCAR(i,j) 1 +SMAT(i,jp)*(AS(2,jp)/c1)**2*SMAT(j,jp)/ZM(i)/ZM(j) OPEN(20,FILE='NEW.FC') CALL WRITEFF(MX3,NA,FCAR) CLOSE(20) write(6,*)'NEW.FC written' stop endif WRITE(*,*)' AUTOMATIC REFINEMENT (Y/N) ?' OK='N' if(nlinepar.gt.ili)then ili=ili+1 call getarg(ili,OK) write(6,*)OK else IF(JOB.NE.'AUTO') READ(*,'(A)')OK endif IF (OK.EQ.'y'.OR.OK.EQ.'Y')THEN DEL=0.05d0 333 WRITE(*,*)' MAXIMUM NUMBER OF ITERATIONS:' READ(*,*)NIT IF (NIT.LT.0)THEN NIT=-NIT WRITE(*,*)' STEP:' READ(*,*)DEL ENDIF WRITE(*,*)' ITERATION IN PROCESS ...' CALL IMPROVE(NIT,S0,NOSC,DEL) WRITE(*,*)' AGAIN (Y/N) ?' READ(*,'(A)')OK IF (OK.EQ.'Y'.OR.OK.EQ.'y')GOTO 333 GOTO 111 ENDIF endif C C Calculate the true mass-weighted S-matrix: DO 210 IA=1,NA T=ZM(IA) DO 210 IM=1,NA 210 SMAT(IA,IM)=SMAT(IA,IM)*T CALL SMATOUT(SMAT,AS,NA,JOB,ili,nlinepar) WRITE(*,*)' S-matrix written out ...' CALL FPC(Q,CHAR,SMAT,NOAT,3*NOAT) WRITE(*,*)' FPC.TAB written ...' CALL TERMO(Q,NA) C C Read-in B-matrix: INQUIRE(FILE='B.MAT',EXIST=LBMAT) C IF(LBMAT)THEN OPEN(15,FILE='B.MAT',FORM='UNFORMATTED') C SMAT used for B-matrix: C FNEW used for B^-1: READ(15) NOBDUM,NA,((SMAT(I,J),J=1,NA),I=1,NA), 1((FNEW(I,J),J=1,NA),I=1,NA) CLOSE(15) WRITE(*,*)' B and B^-1 read in ...' CALL CIFF(NA,FCAR,FNEW,FINT,3,AMULT) ENDIF C IF(NOSC.EQ.0)STOP C SMAT used for B-matrix: C FC IF(LBMAT)CALL CIFF(NA,FINT,SMAT,FNEW,4,ZM) WRITE(*,*) WRITE(*,*)' << Program finished OK >>' STOP END C SUBROUTINE SMATOUT(S,AS,NAT3,JOB,ili,nlinepar) C dX = S dQ S(nat3xnint) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100) CHARACTER*4 JOB CHARACTER*4 num DIMENSION S(MX3,NAT3),AS(3,NAT3) WRITE(*,*)' RANGE OF MODES FOR THE INP FILE (Use the bracket#): ' IRANGE=1 JRANGE=NAT3 if(nlinepar.gt.ili+1)then ili=ili+1 call getarg(ili,num) write(6,*)num read(num,*)IRANGE ili=ili+1 call getarg(ili,num) write(6,*)num read(num,*)JRANGE else IF(JOB.NE.'AUTO')READ(*,*) IRANGE,JRANGE endif NI=NAT3 NI=JRANGE-IRANGE+1 c c look at the phase and make the biggest element always positive: DO 101 J=IRANGE,JRANGE amx=S(1,J) DO 102 I=2,NAT3 102 if(abs(S(I,J)).gt.abs(amx))amx=S(I,J) if(amx.lt.0)then do 103 I=1,NAT3 103 S(I,J)=-S(I,J) endif 101 continue c OPEN(34,FILE='F.INP') WRITE(34,10)NI,NAT3,NAT3/3 10 FORMAT(3I5) OPEN(35,FILE='FILE.X',STATUS='OLD') 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(' Atom Mode X-disp. Y-disp. Z-disp.') DO 1 I=1,NAT3/3 DO 1 J=IRANGE,JRANGE IU=3*(I-1) K=NAT3-J+1 1 WRITE(34,15)I,K,S(IU+1,J),S(IU+2,J),S(IU+3,J) 15 FORMAT(2I5,3F11.6) WRITE(34,11)NI WRITE(34,13)(AS(1,I),I=IRANGE,JRANGE) 13 format(6F11.3) CLOSE(34) RETURN END SUBROUTINE IMPROVE(NIT,S0,NOSC,DEL) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100) COMMON/ITERA/AMULT(MX3),SCFAC(MX3/3), 1ZM(MX3),FCAR(MX3,MX3),FNEW(MX3,MX3), 2SMAT(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX3/3,MX3/3), NA, 3NM(MX3/3) 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) S0=SUMk(NOSC) SCFAC(I)=X0-DEL SM=SUMk(NOSC) SCFAC(I)=X0+DEL SP=SUMk(NOSC) SCFAC(I)=X0 WRITE(16,1600)X0,X0-DEL,X0+DEL 1600 format(3g15.5) WRITE(16,1600)S0,SM,SP IF (SP.LT.S0)THEN S0=SP SCFAC(I)=X0+DEL ICHANGE=1 ENDIF IF (SM.LT.S0)THEN S0=SM SCFAC(I)=X0-DEL ICHANGE=1 ENDIF 2 CONTINUE IF (ICHANGE.EQ.0)DEL=DEL/2.0D0 WRITE(*,6000)IIT,DEL,S0 6000 FORMAT(I4,' del = ',F15.9,', error = ',F15.4) IF (DEL.GT.0.0001)GOTO 1 RETURN END FUNCTION SUMk(NOSC) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX3=3*1100) COMMON/ITERA/AMULT(MX3),SCFAC(MX3/3), 1ZM(MX3),FCAR(MX3,MX3),FNEW(MX3,MX3), 2SMAT(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX3/3,MX3/3), NA, 3NM(MX3/3) DIMENSION FT(MX3,MX3) wtol=0.01d0 c1=1302.828d0 IERR=0 DO 132 I=1,NA 132 AMULT(I)=1.0d0 DO 202 I=1,NOSC SF=1.0d0/SQRT(SCFAC(I)) DO 202 J=1,NM(I) DO 202 K=1,3 202 AMULT(3*(NC(I,J)-1)+K)=SF CALL FAST(NA,FNEW,FT,AMULT) IMOD=1 CALL TRED12(NA,FT,SMAT,Q,IMOD,IERR) IF(IERR.NE.0)STOP SUM=0.0d0 DO 310 I=1,NA T=Q(I) SIGN=1.0d0 IF(T.LT.0.0d0)SIGN=-1.0d0 T=SIGN*SQRT(ABS(T))*c1 TP=ZNU(I) DEL=T-TP if(ABS(TP).le.wtol)DEL=0.0d0 310 SUM=SUM+(T-TP)**2 SUMk=SUM RETURN END SUBROUTINE FPC (V,CHAR,S,NOAT,NA3) C PETR BOUR: LAST UPDATE 4-28-92 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100,MX=MX3/3) DIMENSION V(MX3), S(MX3,NA3), X(3,MX), 1 CHAR(MX), DS(3), SCRT(3) DCOR = 388891.4D0 RCOR = 12221.7D0 OPEN(35,FILE='FILE.X',STATUS='OLD') READ(35,*) READ(35,*)NOAT DO 51 J=1,NOAT 51 READ(35,*)IATDUM,(X(I,J),I=1,3) CLOSE(35) NA=3*NOAT NVIBS=NA C INCLUDE TRANSLATIONS AND ROTATIONS, TOO 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 (ABS(V(I)).GT.0) DIPOL = DT*DCOR/ABS(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) PARAMETER (MX3=3*1100) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) COMMON/DIAG/E 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 C SUBROUTINE TQL12(N,Z,D,IERR,ICODE) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100) DIMENSION Z(MX3,N),D(MX3),E(MX3) REAL*8 MACHEP COMMON/DIAG/E 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 C SUBROUTINE CIFF(NA,A,B,C,IC,U) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100) DIMENSION A(MX3,NA),B(MX3,NA),C(MX3,NA),TEM(MX3,MX3), 1U(NA),UM(MX3) C C IC: 1 ... calculate intrinsic FF C = U.Bt . A .B . U C FINT = SF.B-1t.FCAR.B-1.SF C 2 ... calculate mass-weighted cartesian FF C FCAR=sqrt(M)-1.Bt .FINT.B.sqrt(M)-1 C 3 ... calculate Intrinsic FF and write it C 4 ... calculate cartesian FF and write it in right units C FCAR= Bt .FINT.B C DO 1 I=1,NA 1 UM(I)=U(I) IF(IC.EQ.4)THEN DO 2 I=1,NA 2 UM(I)=1.0d0 ENDIF DO 212 J=1,NA UJ=UM(J) DO 212 I=1,NA TEM(I,J)=0.0d0 DO 212 II=1,NA 212 TEM(I,J)=TEM(I,J)+A(I,II)*B(II,J)*UJ DO 214 I=1,NA UI=UM(I) DO 214 J=1,NA C(I,J)=0.0d0 DO 214 II=1,NA 214 C(I,J)=C(I,J)+B(II,I)*TEM(II,J)*UI IF(IC.EQ.3)THEN OPEN(7,FILE='INTY.FC') WRITE(7,3337)(NA*(NA+1))/2,NA,NA 3337 FORMAT(/,I6,' INTRINSIC FORCE FIELD',I5,' x (',I5,' + 1) / 2') M=6 MI=-5 2071 MI=MI+6 MJ=M-5 WRITE(7,28)(I,I=MJ,MIN(NA,M)) 28 FORMAT(8X,6(I7,5X)) DO 226 I=MI,NA 226 WRITE(7,7000)I,(C(J,I),J=MJ,MIN(I,M)) 7000 FORMAT(I5,3X,6F12.6) IF (M.LT.NA)THEN M=M+6 GOTO 2071 ENDIF CLOSE(7) WRITE(*,*)' INTY.FC with intrinsic coordinates written ...' ENDIF IF(IC.EQ.4)THEN OPEN(20,FILE='NEW.FC') CALL WRITEFF(MX3,NA,C) CLOSE(20) WRITE(*,*)' NEW.FC written out ...' ENDIF RETURN END C SUBROUTINE READFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) C CONST=4.359828d0/0.5291772d0/0.5291772d0 DO 3 I=1,N DO 3 J=I,N 3 FCAR(J,I)=FCAR(J,I)*CONST DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' RETURN END SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) CONST=4.359828d0/0.5291772d0**2 DO 6 I=1,N DO 6 J=1,N 6 FCAR(I,J)=FCAR(I,J)/CONST N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) RETURN END SUBROUTINE FAST(NA,FCAR,FNEW,ZM) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX3=3*1100) DIMENSION FCAR(MX3,NA),FNEW(MX3,NA),ZM(NA) DO 214 I=1,NA UI=ZM(I) DO 214 J=I,NA FNEW(I,J)=FCAR(I,J)*ZM(J)*UI 214 FNEW(J,I)=FNEW(I,J) RETURN END C SUBROUTINE TERMO(W,NA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION W(NA) C Thermodynamics according to Mopac notation Cl=1.4388d0 ELIM=0.1d0 C 1cm-1=1.4388 K C C Zero temperature energy e0=sum over hwi/2 E0=0.0d0 DO 1 I=1,NA IF(I.GT.NA-6)GOTO 1 IF(W(I).LE.ELIM)GOTO 1 E0=E0+0.5d0*W(I) 1 CONTINUE E0=E0*2.8601d0 C cal/mol C C Room temperature T=298.15d0 Cl=1.4388d0/T Qv=0.0d0 E0T=0.0d0 S=0.0d0 SP=0.0d0 C=0.0d0 DO 2 I=1,NA IF(I.GT.NA-6)GOTO 2 IF(W(I).LE.ELIM)GOTO 2 EWJ=0.0d0 IF(W(I)*Cl.LT.200.0d0)EWJ=EXP(-W(I)*Cl) Qv=Qv+1.0d0/(1.0d0-EWJ) E0T=E0T+W(I)*EWJ/(1.0d0-EWJ) S=S+W(I)*EWJ/(1.0d0-EWJ) SP=SP+LOG(1.0d0-EWJ) C=C+W(I)**2*EWJ/(1.0d0-EWJ)**2 2 CONTINUE E0T=E0T*2.8601d0 S=S*2.8601d0/T-SP*1.987216d0 C cal/mol/K C=C*Cl**2*1.987216d0 C cal/mol/K WRITE(16,16000)E0,E0+E0T,S,C,Qv,E0+E0T-T*S 16000 FORMAT(' Thermochemistry parameters',/,80(1H-),/, 1' Vibrational contributions:',/, 1' Linear molecules not implemented!',/, 2' Internal energy 0K : U0 = ',F20.3,' cal/mol',/, 2' Internal energy 298K : U = ',F20.3,' cal/mol',/, 3' Entropy 298K : S = ',F20.3,' cal/mol/K',/, 4' Heat capacity 298K : C = ',F20.3,' cal/mol/K',/, 5' Partition function : Qv = ',F20.3,/, 6' Gibbs free energy : G = ',F20.3,' cal/mol/K',/, 780(1H-),/,/) RETURN END C SUBROUTINE PROJECT(F,N0,N) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100,MX=MX3/3) DIMENSION F(N0,N) C C Using the common here to save memory for A and TEM COMMON/ITERA/AMULT(MX3),SCFAC(MX), 1ZM(MX3),FCAR(MX3,MX3),A(MX3,MX3), 2TEM(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX),NA,NM(MX) WRITE(*,*)'Projecting transl/rotations from the force field...' CALL DOMA(A,N0,N) DO 1 I=1,N DO 1 J=1,N TEM(I,J)=0.0d0 DO 1 II=1,N 1 TEM(I,J)=TEM(I,J)+A(I,II)*F(II,J) DO 2 I=1,N DO 2 J=1,N F(I,J)=0.0d0 DO 2 II=1,N 2 F(I,J)=F(I,J)+TEM(I,II)*A(J,II) RETURN END SUBROUTINE DOMA(A,NDIM,N) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*1100,MX=MX3/3) DIMENSION C(3,MX),A(NDIM,N) C C Using the common here to save memory for TEM COMMON/ITERA/AMULT(MX3),SCFAC(MX), 1ZM(MX3),FCAR(MX3,MX3),AD(MX3,MX3), 2TEM(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX),NA,NM(MX) AMACH=0.00000000000001d0 NAT=N/3 OPEN(4,FILE='FILE.X',STATUS='OLD') READ(4,*) READ(4,*)NAT DO 11 I=1,NAT 11 READ(4,*)ID,(C(II,I),II=1,3) CLOSE(4) DO 12 I=1,N DO 12 J=1,6 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 A(3*I-2,4)=-C(2,I) A(3*I-1,4)= C(1,I) A(3*I-1,5)=-C(3,I) A(3*I ,5)= C(2,I) A(3*I-2,6)= C(3,I) 1 A(3*I ,6)=-C(1,I) N6=6 CALL ORT(A,NDIM,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1000.0d0*AMACH)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE WRITE(*,*)N6,' coordinates projected' DO 2 I=1,N DO 2 J=1,N TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END C ============================================================ SUBROUTINE ORT(A,NDIM,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(NDIM,N6) AMACH=0.00000000000001d0 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/SQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP RETURN END