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 (MX=2000,MX3=3*MX) CHARACTER*1 OK CHARACTER*4 JOB LOGICAL LBMAT DIMENSION ZIM(MX),CHR(MX),FINT(MX3,MX3),FCAR(MX3,MX3) COMMON/ITERA/AMULT(MX3),SCFAC(MX),ZM(MX3), 1FNEW(MX3,MX3),SMAT(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX),NA,NM(MX), 2AS(3,MX3) C wtol=0.01d0 c1=1302.828d0 it0=mclock() 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 if(NA.gt.MX3)then write(6,*)'too many atoms' stop endif do i=1,3 do j=1,NA AS(i,j)=0.0d0 enddo enddo 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) (CHR(I),I=1,NOAT) READ(15,*) (CHR(I),I=1,NOAT) CHATOT=0.0d0 DO 77 I=1,NOAT 77 CHATOT=CHATOT+CHR(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' write(6,6009)mclock()-it0 6009 format(' Time: ',i10,' msec') 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) write(6,6009)mclock()-it0 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 ie=0 ADUM=SUMk(NOSC,1,ie,FCAR) write(6,6009)mclock()-it0 if(ie.gt.0)then WRITE(6,6667) 6667 FORMAT(' 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,FCAR) 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,i1,i2) WRITE(*,*)' S-matrix written out ...' CALL FPC(Q,CHR,SMAT,NOAT) WRITE(*,*)' FPC.TAB written ...' CALL TERMO(Q,NA,i1,i2) 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 >>' write(6,6009)mclock()-it0 STOP END C SUBROUTINE SMATOUT(S,AS,NAT3,JOB,ili,nlinepar,IRANGE,JRANGE) C dX = S dQ S(nat3xnint) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=2000,MX3=3*MX) CHARACTER*4 JOB CHARACTER*6 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(3I7) 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(I7,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(2I7,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,FCAR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=2000,MX3=3*MX) dimension FCAR(MX3,MX3) COMMON/ITERA/AMULT(MX3),SCFAC(MX), 1ZM(MX3),FNEW(MX3,MX3), 2SMAT(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX), NA, 3NM(MX), 2AS(3,MX3) 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 ie=0 DO 2 I=1,NOSC X0=SCFAC(I) S0=SUMk(NOSC,0,ie,FCAR) SCFAC(I)=X0-DEL SM=SUMk(NOSC,0,ie,FCAR) SCFAC(I)=X0+DEL SP=SUMk(NOSC,0,ie,FCAR) 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.0001d0)GOTO 1 RETURN END FUNCTION SUMk(NOSC,iwr,iee,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=2000,MX3=3*MX) dimension FCAR(MX3,MX3) COMMON/ITERA/AMULT(MX3),SCFAC(MX), 1ZM(MX3),FNEW(MX3,MX3), 2SMAT(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX), NA, 3NM(MX), 2AS(3,MX3) DIMENSION FINT(MX3,MX3) c wtol=0.01d0 c1=1302.828d0 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 i=1,NA do j=1,NA FNEW(i,j)=FINT(i,j) enddo enddo ENDIF if(IWR.gt.0)WRITE(6,*)' New cartesian force field done ' IMOD=2 CALL TRED12(NA,FNEW,SMAT,Q,IMOD,IERR) IF(IERR.NE.0)STOP if(iwr.gt.0)WRITE(6,*)' Diagonalization OK' c RMSW=0.0D0 ie=0 DO 310 I=1,NA T=Q(I) SIGN=1.0d0 C Make imaginary frequencies negative: IF(T.LT.0.0d0)SIGN=-1.0d0 T=SIGN*SQRT(ABS(T))*c1 Q(I)=T TP=ZNU(I) DEL=T-TP if(ABS(TP).le.wtol)then DEL=0.0d0 else ie=ie+1 endif RMSW=RMSW+DEL**2 WRITE(16,1230) I, T,TP,DEL 1230 FORMAT(I6,3F15.2) AS(1,I)=T AS(2,I)=TP 310 AS(3,I)=T-TP iee=ie if(ie.gt.0)then RMSW=DSQRT(RMSW/DFLOAT(ie)) else if(NA.gt.0)RMSW=DSQRT(RMSW/DFLOAT(NA)) endif 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 if(iwr.gt.0)WRITE(*,5555) 5555 FORMAT(1X,'************ NEW FREQUENCIES ******************') DO 311 J=NA,1,-1 a3=0.0d0 if(AS(2,NA-J+1).gt.wtol)a3=AS(3,NA-J+1) 311 if(iwr.gt.0)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) if(iwr.gt.0)WRITE(*,5551) 5551 FORMAT(1X,'************************************************') if(iwr.gt.0)WRITE(6,6667)ie,RMSW,AK 6667 FORMAT(i6,'experimental frequencies found',/, 1 ' >>>>><< >> ',F8.2,' <<<<<< ',F8.6) c SUMk=RMSW RETURN END SUBROUTINE FPC(V,CHR,S,NOAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=2000,MX3=3*MX) real*8 V(*),S(MX3,MX3),X(3,MX),CHR(*),DS(3),DR(3) DCOR = 388891.4D0 RCOR = 12221.7D0 OPEN(35,FILE='FILE.X',STATUS='OLD') READ(35,*) READ(35,*) DO 51 J=1,NOAT 51 READ(35,*)IATDUM,(X(I,J),I=1,3) CLOSE(35) 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,3*NOAT DO 784 IKY=1,3 DS(IKY)=0.0D0 DR(IKY)=0.0D0 DO 784 JR=1,NOAT DS(IKY)=DS(IKY)+CHR(JR)*S(3*(JR-1)+IKY,I) DO 784 IB=1,3 DO 784 IG=1,3 784 DR(IKY)=DR(IKY)+EPS(IKY,IB,IG)*X(IB,JR)*S(3*(JR-1)+IG,I) 1 *CHR(JR) DT =DS(1)*DS(1)+DS(2)*DS(2)+DS(3)*DS(3) ROT=DS(1)*DR(1)+DS(2)*DR(2)+DS(3)*DR(3) IF(V(I).GT.1.0d0)then WRITE(13,85) I,V(I),DT*DCOR/ABS(V(I)),ROT*RCOR 85 FORMAT(I4,F8.2,' ',F13.6,' ',f13.6) endif 42 continue 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 (MX=2000,MX3=3*MX) 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 (MX=2000,MX3=3*MX) 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 (MX=2000,MX3=3*MX) 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 (MX=2000,MX3=3*MX) 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,i1,i2) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION W(NA) character*2 key logical lex C Thermodynamics according to Mopac notation Cl=1.4388d0 w0=0.0d0 T=298.15d0 C 1cm-1=1.4388 K c inquire(file='THERMO.OPT',exist=lex) if(lex)then open(99,file='THERMO.OPT') 3 read(99,990,end=99,err=99)key 990 format(a2) if(key.eq.'w0')read(99,*)w0 if(key.eq.'TE')read(99,*)T goto 3 99 close(99) endif C Ecm=0.0d0 Cl=1.4388d0/T Qv=0.0d0 Ucm=0.0d0 Scm=0.0d0 C=0.0d0 c DO 2 I=i1,i2 o=dabs(max(W(I),0.0d0)) ot=o+1.0d-9 EWJ=0.0d0 EWJt=0.0d0 IF(o*Cl.LT.200.0d0)EWJ=EXP(-o*Cl) IF(o*Cl.LT.200.0d0)EWJt=EXP(-ot*Cl*(1.0d0+w0/ot)) Qv=Qv+1.0d0/(1.0d0-EWJ) Ecm=Ecm+o* 0.5d0 Ucm=Ucm+o*(0.5d0+EWJ/(1.0d0-EWJ)) Scm=Scm+ot*EWJt/(1.0d0-EWJt)/T-0.69503476d0*Log(1.0d0-EWJt) 2 C=C+o**2*EWJ/(1.0d0-EWJ)**2 c Gcm=Ucm-T*Scm Eau=Ecm/219470.0d0 Uau=Ucm/219470.0d0 Sau=Scm/219470.0d0 Gau=Gcm/219470.0d0 Ukcal=Uau*627.5d0 Skcal=Sau*627.5d0 Gkcal=Gau*627.5d0 Ekcal=Eau*627.5d0 C cal/mol/K C=C*Cl**2*1.987216d0 C cal/mol/K WRITE(16,16001)T,w0 16001 FORMAT(/,' Thermochemistry parameters',/, 1 ' Temperature :',f12.3,' K',/, 1 ' w0 :',f12.3,' cm-1',/, 224x,' kcal/mol cm-1 hartree',/,80(1H-)) WRITE(16,16000)'Internal energy 0K U0:',Ekcal,Ecm,Eau WRITE(16,16000)'Internal energy TK U :',Ukcal,Ucm,Uau WRITE(16,16000)'Entropy S :',Skcal,Scm,Sau WRITE(16,16000)'Gibbs G :',Gkcal,Gcm,Gau WRITE(16,16000)'Heat capacity C :',C/1000.0d0 16000 FORMAT(a22,F14.5,F12.1,F12.6) WRITE(16,16002) 16002 format(80(1H-)) RETURN END C SUBROUTINE PROJECT(F,N0) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=2000,MX3=3*MX) DIMENSION F(N0,N0) C C Using the common here to save memory for A and TEM COMMON/ITERA/AMULT(MX3),SCFAC(MX), 1ZM(MX3),A(MX3,MX3), 2TEM(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX),NA,NM(MX), 2AS(3,MX3) N=NA WRITE(*,*)'Projecting transl/rotations from the force field...' CALL DOMA 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 IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=2000,MX3=3*MX) DIMENSION C(3,MX) C C Using the common here to save memory for TEM COMMON/ITERA/AMULT(MX3),SCFAC(MX), 1ZM(MX3),A(MX3,MX3), 2TEM(MX3,MX3),Q(MX3),ZNU(MX3),NC(MX,MX),NA,NM(MX), 2AS(3,MX3) N=NA 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,MX3,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 function mclock() integer*4 mclock mclock=0 return end