PROGRAM NEW4BIG C Harmonic force field diagonalization, c works in cartesian coordinates, for one molecule only c for giant FFs IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*21000,MX=MX3/3,NV0=5000) CHARACTER*1 OK CHARACTER*4 JOB DIMENSION ZIM(MX),CHAR(MX) REAL*4 SMAT COMMON/ITERA/AMULT(MX3),SCFAC(NV0),ZM(MX3), 1SMAT(MX3,NV0),Q(MX3),ZNU(MX3),NC(NV0,MX),NA,NM(MX), 2AS(3,MX3),EMAT(NV0) C 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!',/,/, 2 ' For giant force fields',/,/, 3 ' Petr Bour, Prague 2-08',/,/) C OPEN(17,FILE='FILE.X') READ(17,*) READ(17,*)NOAT CLOSE(17) NA=3*NOAT if(NA.gt.MX3)call report('too many atoms') do 2 i=1,3 do 2 j=1,NA 2 AS(i,j)=0.0d0 C C Read-in scalefactors: OPEN(15,FILE='FTRY.INP') READ(15,1500)JOB 1500 FORMAT(A4) READ(15,*)NOSC if(NOSC.gt.NV0)call report('too many scale factors') 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)call report(' Only one molecule allowed!') 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) DO 220 I=1,NOAT DO 220 J=1,3 220 ZM(3*I-3+J)=1/SQRT(ZIM(I)) WRITE(*,778)BM 778 FORMAT(' TOTAL MASS ',F20.6) 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 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 ievs1=0 ievs2=0 ADUM=SUMk(NOSC,1,ie,ievs1,ievs2) if(ie.gt.0)then WRITE(6,6667) 6667 FORMAT(' NEW SCALEFACTORS (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 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 endif C Calculate the true mass-weighted S-matrix and writ eon disk: DO 210 IA=1,NA T=ZM(IA) DO 210 IM=1,ievs2-ievs1+1 210 SMAT(IA,IM)=SMAT(IA,IM)*real(T) CALL SMATOUT(SMAT,NA,ievs1,ievs2,EMAT) WRITE(*,*)' S-matrix written out ...' c c fixed partial charge model: CALL FPC(EMAT,CHAR,SMAT,NOAT,ievs1,ievs2) WRITE(*,*)' FPC.TAB written ...' CALL TERMO(EMAT,ievs2-ievs1+1) C close(16) call report('<< Program finished OK >>') END C SUBROUTINE SMATOUT(S,NAT3,ievs1,ievs2,EMAT) C dX = S dQ S(nat3xnint) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*21000) real*4 S(MX3,*),amx DIMENSION EMAT(*) NI=ievs2-ievs1+1 c c look at the phase and make the biggest element always positive: DO 101 J=1,NI 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(3I6) 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=1,NI 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)(EMAT(I),I=1,NI) 13 format(6F11.3) CLOSE(34) RETURN END FUNCTION SUMk(NOSC,iwr,iee,ievs1,ievs2) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NV0=5000,MX3=3*21000) dimension j33(MX3),a33(MX3) REAL*4 SMAT COMMON/ITERA/AMULT(MX3),SCFAC(NV0), 1ZM(MX3), 2SMAT(MX3,NV0),Q(MX3),ZNU(MX3),NC(NV0,MX3/3), NA, 3NM(MX3/3), 2AS(3,MX3),EMAT(NV0) logical lex,LUP,lvar character*5 s5 CONST=4.359828d0/0.5291772d0/0.5291772d0 c DO 201 I=1,NA 201 AMULT(I)=1.0d0 IF(NOSC.GT.0)THEN 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 ENDIF c WRITE(*,*)' Making scaled mass-weighed FF' OPEN(20,FILE='AFILE.FC',STATUS='OLD',form='unformatted') OPEN(21,FILE='A3.SCR',form='unformatted') c OPEN(22,FILE='A3.TXT') c write(22,*)NA do 20 i=1,NA fi=ZM(i)*AMULT(i)*CONST read(20)nn,(j33(j),j=1,nn),(a33(j),j=1,nn) write(21)nn,(j33(j),j=1,nn),(a33(j)*fi*ZM(j33(j)) 1*AMULT(j33(j)) ,j=1,nn) c write(22,*)nn c do 20 j=1,nn c write(22,*)j33(j),a33(j) 20 continue close(20) close(21) close(22) WRITE(*,*)' Cartesian FF rewritten' wtol=0.01d0 c1=1302.828d0 c M: subspace dimension: c NE number of ev to b efound c ELMAX .. max eigenvalue c IERR c MV number of ev than has been found c lvar variable subspace dimension c ievs1,ievs2 - interval of evs written into S-matrix c IVST starting eigenvalue c IPULAY use the pulay gradient vectors c ibl eigenvalues in the block at once c ibla start its always with unit vector M=10 ELMAX=9000.0d0 ipr=0 NE=NA GRADCONV=1.0d-12 ECONV=1.0d-12 LUP=.false. lvar=.false. ievs1=1 ievs2=NE OLIMIT=1.0d-6 MAXIT=10000 IVST=1 ibl=0 ibla=0 inquire(file='NEW4BIG.OPT',exist=lex) if(lex)then open(77,file='NEW4BIG.OPT') 111 read(77,770,end=112,err=112)s5 770 format(a5) if(s5(1:5).eq.'ELMAX')read(77,*)ELMAX if(s5(1:3).eq.'IPR')read(77,*)ipr if(s5(1:3).eq.'LUP')read(77,*)LUP if(s5(1:2).eq.'NE')read(77,*)NE if(s5(1:3).eq.'GRA')read(77,*)GRADCONV if(s5(1:3).eq.'ECO')read(77,*)ECONV if(s5(1:4).eq.'LVAR')read(77,*)lvar if(s5(1:4).eq.'IEVS')read(77,*)ievs1,ievs2 if(s5(1:4).eq.'MSUB')read(77,*)M if(s5(1:4).eq.'OLIM')read(77,*)OLIMIT if(s5(1:5).eq.'MAXIT')read(77,*)MAXIT if(s5(1:4).eq.'IVST')read(77,*)IVST if(s5(1:5).eq.'PULAY')read(77,*)IPULAY if(s5(1:5).eq.'BLOCK')read(77,*)ibl if(s5(1:5).eq.'BLANK')read(77,*)ibla goto 111 112 close(77) endif if(ievs2-ievs1+1.gt.NV0)call report('too many evs to record') if(ibl.eq.0)then call MITIN3(NA,M,NE,ELMAX,IERR,MV,LUP,lvar,ipr,c1, 1 GRADCONV,ievs1,ievs2,OLIMIT,ECONV,MAXIT,IVST,IPULAY,ibla) else call MITINB(NA,M,NE,ELMAX,IERR,MV,LUP,ipr,c1, 1 GRADCONV,ievs1,ievs2,OLIMIT,ECONV,MAXIT,IVST,IPULAY,ibl) endif IF(IERR.NE.0)STOP if(iwr.gt.0)WRITE(6,*)' Diagonalization OK' c RMSW=0.0D0 ie=0 DO 310 I=1,MV 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=MV,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 I7,' (',I7,') ...... ',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,CHAR,S,NOAT,ievs1,ievs2) C PETR BOUR: LAST UPDATE 4-28-92/2008 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*21000,MX=MX3/3) real*4 S(MX3,*) DIMENSION V(*),X(3,MX),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) NVIBS=ievs2-ievs1+1 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)*dble(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)* 1dble(S(3*(JR-1)+IG,I))*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(MX3,N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) 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 (MX3,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(MX3,N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(MX3,N),D(MX3),E(MX3) 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 C SUBROUTINE TERMO(W,NA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION W(*) 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 MITIN3(N,M,NE,ELMAX,IERR,MV,LUP,lvar,ipr,c1, 1GRADCONV,ievs1,ievs2,OLIMIT,ECONV,MAXIT,IVST,IPULAY,ibla) C C by Petr Bour, Prague 1994/2008 C C This subroutine finds few of lowest eigenvectors for a C very large matrix A, elements done on the fly C Diagonal terms are already divided by two in A c N matrix dimension C M is soft dimension of the subspace where direct diagonalization is done C NE is number of eigenvectors to be found C MV is the number of eigenvectors found C C For details look at C Alexander V. Mitin, J.Comp.Chem., Vol. 15, No.7, (1994), 747-751 C (Warning: The paper is full of missprints; the recommended procedure C is modified here.) C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*21000,MX=MX3/3,NV0=5000) PARAMETER (NSS=50,NM0=MX3) DIMENSION X(NM0*NSS),E(NSS),F(NSS,NSS),AIJV(NM0),DUM(1), 1AXV(NM0*NSS),Z(NSS,NSS),Y(NM0),YNEW(NM0),YNEW2(NM0), 1j33(NM0),a33(NM0),EWORK(MX3),DD(MX3) LOGICAL HW,lvar,lup COMMON/OPTS/HW,IERM REAL*4 SMAT COMMON/ITERA/AMULT(MX3),SCFAC(NV0),ZM(MX3), 1SMAT(MX3,NV0),D(MX3),ZNU(MX3),NC(NV0,MX),NA,NM(MX), 2AS(3,MX3),EMAT(NV0) c MV=0 HW=ipr.gt.1 ELAST=-1.0d6 GRADOLD=-1.0d9 IERR=0 ANORM=0.0d0 DIAG=0.0d0 OPEN(45,FILE='CFI.SCR',FORM='UNFORMATTED') AF=1.0d0/DSQRT(DBLE(N)) DO 291 I=1,N 291 YNEW2(I)=AF itime0=mclock() amax=0.0d0 open(49,file='A3.SCR',form='unformatted') DO 2229 I=1,N read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) call vz(AIJV,0.0d0,I) do 26 ii=1,nn 26 AIJV(j33(ii))=a33(ii) DIAG=DIAG+ABS(AIJV(I)) DD(I)=AIJV(I) DO 2239 J=1,I if(ABS(AIJV(J)).gt.amax)amax=AIJV(J) 2239 ANORM=ANORM+ABS(AIJV(J)) 2229 continue write(6,*)'maximum component of A : ',amax write(6,*)'passed dimension of subspace: ',M IF(N.LT.M.OR.NSS.LT.M)THEN M=MIN(N,NSS) WRITE(6,*)' Sub-space dimension decreased to = ',M ENDIF MSAVE=M ANORM=2.0d0*ANORM/(DBLE(N*(N-1))) DIAG=DIAG/DBLE(N)*2.0d0 ALCONV=GRADCONV*MAX(DIAG,ANORM)*DBLE(N) WRITE(6,3003)ANORM,DIAG, ALCONV,ECONV,OLIMIT,MAXIT,N,M,IVST,NE 3003 FORMAT(' Average magnitude of matrix el. :',F20.10,/, 1 ' Average diagonal element :',F20.10,/, 2 ' Convergence criterium grad :',F20.10,/, 2 ' Convergence criterium E :',F20.10,/, 3 ' Cut-off norm for subspace :',F20.10,/, 4 ' Maximum number of iterations :',I20,/, 5 ' Dimension of the matrix space :',I20,/, 6 ' Dimension of the iterative space:',I20,/, 8 ' Total eigenvalue :',I20,/, 8 ' Number of vectors to be found :',I20) if(LUP)then WRITE(6,4004)ELMAX 4004 format( ' Lower limit for eigenvalues :',F20.10,/,/) else WRITE(6,4001)ELMAX 4001 format( ' Upper limit for eigenvalues :',F20.10,/,/) endif itime0=mclock() OPEN(4,FILE='E.SCR',FORM='UNFORMATTED') NDONE=0 READ(4,END=444)NDONE DO 23 IEE=1,NDONE READ(4)D(IEE) 23 if(ipr.gt.0)write(6,*)D(IEE) close(4) write(6,*)NDONE, ' eigenvalues found on disk' 444 IERM=NDONE+1 if(IVST.GT.IERM)then write(6,6009) 6009 format('Requested starting eigenvalue must be reduced') else IERM=IVST endif WRITE(6,6008)IERM 6008 FORMAT(' Calculation started from eigenvalue number ',I4) C Starting dimension of the subspace: MSTART=M c c Loop over eigenvectors: DO 5555 IE=IERM,MIN(N,NE) c C Estimate starting eigenvalue and eigenvector as unit vector: EL=0.0d0 call VZ(YNEW,AF,N) C c for advanced, take previous guess IF(IE.GT.IERM.and.ibla.gt.0)then EL=D(IE-1) do 299 I=1,N 299 YNEW(I)=YNEW2(I) endif c if already done, take it if(IE.LE.NDONE)then REWIND 45 DO 91 IN=1,IE-1 91 READ(45) read(45)(YNEW(I),I=1,N) endif c if(HW)write(6,*)'guess:' if(HW)write(6,*)EL if(HW)WRITE(6,3004)(YNEW(J),J=1,N) c C Orthonormalize to previous eigenvectors Y1..YIE-1: REWIND 45 DO 71 IN=1,IE-1 READ(45)(Y(I),I=1,N) SN=spv(Y,YNEW,N,0,0) DO 71 I=1,N 71 YNEW(I)=YNEW(I)-SN*Y(I) call normv(YNEW,1,N) DO 51 I=1,N 51 X(I)=YNEW(I) c ITER=0 DELL=0.0d0 M=MSTART C C Start of the iterative cycle: C 9999 ITER=ITER+1 IF(ITER.GT.MAXIT)THEN WRITE(6,*)' Maximum number of iterations exceeded' IERR=8888 GOTO 8888 ENDIF IF(M+IE-1.GT.N)THEN WRITE(6,*)' Subspace-dimension has to be decreased' M=N-IE+1 IF(M.LT.1)THEN WRITE(6,*)' No subspace for the eigenvector left !' WRITE(6,*)' Try to increase the subspace-dimension.' GOTO 6666 ENDIF ENDIF C Form gradient vectors Xk,k=2..M, new basis: {X1,..,Xk} C Xk+1=AXk-ELXk (Xk,K>1 generated by X1) call vz(AXV,0.0d0,N*M) DO 3 K=2,M rewind 49 K2N=(K-2)*N K1N=(K-1)*N DO 20 I=1,N c read non-zero components read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) DO 20 ii=1,nn J=j33(ii) AXV(J+K2N)=AXV(J+K2N)+a33(ii)*X(I+K2N) 20 AXV(I+K2N)=AXV(I+K2N)+a33(ii)*X(J+K2N) c C Gradient vectors DO 2031 I=1,N 2031 X(I+K1N)=AXV(I+K2N)-EL*X(I+K2N) if(IPULAY.eq.1.and.K.eq.M)then DO 2131 I=1,N X(I+K1N)=AXV(I+K2N)-EL*X(I+K2N) 2131 if(DABS(DD(I)-EL).gt.1.0d-6)X(I+K1N)=AXV(I+K2N)/(DD(I)-EL) endif c C Orthogonalize vector XK to vectors XIN, IN=1..K-1 (including trial eigv): DO 7 IN=1,K-1 IN1=(IN-1)*N SN=spv(X,X,N,IN1,K1N) DO 7 I=1,N 7 X(I+K1N)=X(I+K1N)-SN*X(I+IN1) c C Orthogonalize vector XK to previous eigenvectors Xin, in=1..IE-1 REWIND 45 DO 72 IN=1,IE-1 READ(45)(Y(I),I=1,N) SN=spv(Y,X,N,0,K1N) DO 72 I=1,N 72 X(I+K1N)=X(I+K1N)-SN*Y(I) c C Normalize vector K: ANORM=spv(X,X,N,K1N,K1N) c C Check on convergence on gradient=norm of second trial vector: IF(K.EQ.2)THEN CONV=ANORM IF(ABS(CONV).LT.ALCONV)GOTO 8888 ENDIF ANORM=1.0d0/DSQRT(ANORM) DO 5 I=1,N 5 X(I+K1N)=X(I+K1N)*ANORM 3 CONTINUE rewind 49 M1N=(M-1)*N DO 204 I=1,N read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) DO 204 ii=1,nn J=j33(ii) AXV(I+M1N)=AXV(I+M1N)+a33(ii)*X(J+M1N) 204 AXV(J+M1N)=AXV(J+M1N)+a33(ii)*X(I+M1N) C Now looking for new eigenvector in the subspace of {Xk} C Form the sub-space matrix F=X.A.X: DO 21 II=1,M II1N=(II-1)*N DO 21 JJ=1,II JJ1N=(JJ-1)*N SUM=spv(AXV,X,N,II1N,JJ1N) F(JJ,II)=SUM 21 F(II,JJ)=SUM c C Diagonalize the sub-space directly: call TRED12(NSS,M,F,Z,E,2,IERR,EWORK) C Take the lowest/highest eigenvalue: if(LUP)then iu1=1 iu2=2 else iu1=M iu2=max(1,M-1) endif c C New approximation to this and following eigenvector: DO 24 I=1,N SUM=0.0d0 SUM2=0.0d0 DO 13 J=1,M SUM2=SUM2+Z(J,iu2)*X(I+(J-1)*N) 13 SUM=SUM+Z(J,iu1)*X(I+(J-1)*N) YNEW2(I)=SUM 24 YNEW(I)=SUM c call normv(YNEW,1,N) c c Get actual eigenvalue call vz(AXV,0.0d0,N) rewind 49 DO 206 I=1,N read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) DO 206 ii=1,nn J=j33(ii) AXV(I)=AXV(I)+a33(ii)*YNEW(J) 206 AXV(J)=AXV(J)+a33(ii)*YNEW(I) EL=spv(AXV,YNEW,N,0,0) wcm=sqrt(abs(EL))*c1 if(EL.lt.0.0d0)wcm=-wcm if(ipr.gt.0)WRITE(6,6007)ITER,wcm,CONV,M 6007 FORMAT(' It=',I5,' Lambda=',F20.10,' Grad=',F20.10,' M=',I3) DELL=ABS(ELAST-wcm) if(DELL.lt.ECONV)goto 8888 ELAST=wcm c c rewrite to first vector do 250 I=1,N 250 X(I)=YNEW(I) c IF(M.EQ.1)THEN WRITE(6,*)'Last eigenvector left' GOTO 8888 ENDIF c MNEW=0 if(lvar)then C Check if the sub-space is not redundant: MNEW=M DO 50 J=M,4,-1 FF=Z(J,iu1)**2 IF(FF.GT.OLIMIT)GOTO 4444 50 MNEW=J-1 4444 M=MAX(MNEW,3) C C Check if the sub-space should be increased: FF=Z(M,iu1)**2 IF(FF.GT.OLIMIT)M=MIN(M+1,NSS) IF(ITER.GT.1)THEN IF(CONV.GE.GRADOLD)M=MIN(M+1,NSS) ENDIF endif C GRADOLD=CONV GOTO 9999 8888 wm=dsqrt(dabs(EL))*c1 if(EL.lt.0.0d0)wm=-wm WRITE(6,6000)IE,ITER-1,CONV,wm,ALCONV,DELL,MNEW 6000 Format( 1' *************** The ',I3,'. eigenvalue ******************',/, 2' Iter',I6,' Gradient : ',F20.10,' Lambda : ',F20.10,/, 3' ',6X,' Limit grad.: ',F20.10,' Del Lam.: ',F20.10,/, 4' Final sub-space dim. : ',I20,/,/) write(6,6001)mclock()-itime0 6001 format(' Time: ',I10,' msec.') MV=IE D(IE)=EL c C Orthonormalize to previous eigenvectors Y1..YIE-1: REWIND 45 DO 81 IN=1,IE-1 READ(45)(Y(I),I=1,N) SN=spv(Y,YNEW,N,0,0) DO 81 I=1,N 81 YNEW(I)=YNEW(I)-SN*Y(I) call normv(YNEW,1,N) c C Save values & vectors into scratch-files: REWIND 45 DO 1003 J=1,IE-1 1003 READ(45) WRITE(45)(YNEW(I),I=1,N) OPEN(4,FILE='E.SCR',FORM='UNFORMATTED') WRITE(4)IE DO 2 IEE=1,IE 2 WRITE(4)D(IEE) CLOSE(4) IF(HW)THEN WRITE(6,*)' The eigenvector:' WRITE(6,3004)(YNEW(J),J=1,N) 3004 FORMAT(7F10.5) ENDIF c j1a=1 y1a=abs(YNEW(1)) do 4 J=2,N if(abs(YNEW(J)).gt.y1a)then j1a=J y1a=abs(YNEW(J)) endif 4 continue j2a=1 y2a=abs(YNEW(1)) if(j2a.eq.j1a)then j2a=2 y2a=abs(YNEW(2)) endif do 31 J=2,N if(abs(YNEW(J)).gt.y2a.and.J.ne.j1a)then j2a=J y2a=abs(YNEW(J)) endif 31 continue write(6,9006)YNEW(j1a),j1a,YNEW(j2a),j2a 9006 format(' Maximum contributions: ',2(F8.4,' of ',i6,'; ')) c if((LUP.and.wm.LT.ELMAX).or.(.NOT.LUP.AND.wm.GT.ELMAX))goto 6666 5555 CONTINUE 6666 M=MSAVE c REWIND 45 c OPEN(4,FILE='E.SCR',FORM='UNFORMATTED') read(4) ii=0 DO 10031 IE=1,MV READ(4)EL wcm=sqrt(abs(EL))*c1 if(EL.lt.0.0d0)wcm=-wcm READ(45)(YNEW(I),I=1,N) if(IE.ge.ievs1.and.IE.le.ievs2)then ii=ii+1 do 1 I=1,N 1 SMAT(I,ii)=real(YNEW(I)) c write(6,6688)real(ii),(SMAT(I,ii),I=1,N) 6688 format(10f8.4) EMAT(ii)=wcm endif 10031 continue ievs2=ievs1+ii-1 write(6,*)ii,'vectors recorded into S-matrix' CLOSE(4) CLOSE(45) close(49) c write(6,9009)MV 9009 format(I6,' eigenvectors found') RETURN END c SUBROUTINE VZ(V,Z,N) REAL*8 V(*),Z INTEGER*4 I,N DO 1 I=1,N 1 V(I)=Z RETURN END subroutine report(s) character*(*) s write(6,*)s stop end function spv(A,B,N,N1,M1) real*8 spv,A(*),B(*),S integer*4 I,N S=0.0d0 DO 1 I=1,N 1 S=S+A(I+N1)*B(I+M1) spv=S return end subroutine normv(A,N1,N2) real*8 A(*),ANORM integer*4 N1,N2,I ANORM=0.0d0 DO 1 I=N1,N2 1 ANORM=ANORM+A(I)**2 ANORM=1.0d0/DSQRT(ANORM) DO 2 I=N1,N2 2 A(I)=A(I)*ANORM return end SUBROUTINE MITINB(N,M,NE,ELMAX,IERR,MV,LUP,ipr,c1, 1GRADCONV,ievs1,ievs2,OLIMIT,ECONV,MAXIT,IVST,IPULAY,ibl) c c more like Davidson block diagonalization C C by Petr Bour, Prague 1994/2008 C C This subroutine finds few of lowest eigenvectors for a C very large matrix A, elements done on the fly C Diagonal terms are already divided by two in A C M is soft dimension of the subspace where direct diagonalization is done C NE is number of eigenvectors to be found C MV is the number of eigenvectors found C C For details look at C Alexander V. Mitin, J.Comp.Chem., Vol. 15, No.7, (1994), 747-751 C (Warning: The paper is full of missprints; the recommended procedure C is modified here.) C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX3=3*21000,MX=MX3/3,NV0=5000) PARAMETER (NSS=50,NM0=MX3) DIMENSION X(NM0*NSS),E(NSS),F(NSS,NSS),AIJV(NM0),DUM(1), 1AXV(NM0*NSS),Z(NSS,NSS),Y(NM0),YNEW(NM0),EL(NSS), 1j33(NM0),a33(NM0),EWORK(MX3),DD(MX3) LOGICAL HW,lup COMMON/OPTS/HW,IERM REAL*4 SMAT COMMON/ITERA/AMULT(MX3),SCFAC(NV0),ZM(MX3), 1SMAT(MX3,NV0),D(MX3),ZNU(MX3),NC(NV0,MX),NA,NM(MX), 2AS(3,MX3),EMAT(NV0) c MV=0 HW=.false. ELAST=-1.0d6 IERR=0 ANORM=0.0d0 DIAG=0.0d0 OPEN(45,FILE='CFI.SCR',FORM='UNFORMATTED') AF=1.0d0/DSQRT(DBLE(N)) itime0=mclock() amax=0.0d0 open(49,file='A3.SCR',form='unformatted') DO 2229 I=1,N read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) call vz(AIJV,0.0d0,I) do 26 ii=1,nn 26 AIJV(j33(ii))=a33(ii) DIAG=DIAG+ABS(AIJV(I)) DD(I)=AIJV(I) DO 2239 J=1,I if(ABS(AIJV(J)).gt.amax)amax=AIJV(J) 2239 ANORM=ANORM+ABS(AIJV(J)) 2229 continue write(6,*)'maximum component of A : ',amax write(6,*)'passed dimension of subspace: ',M IF(N.LT.M.OR.NSS.LT.M)THEN M=MIN(N,NSS) WRITE(6,*)' Sub-space dimension decreased to = ',M ENDIF MSAVE=M ANORM=2.0d0*ANORM/(DBLE(N*(N-1))) DIAG=DIAG/DBLE(N)*2.0d0 ALCONV=GRADCONV*MAX(DIAG,ANORM)*DBLE(N) WRITE(6,3003)ANORM,DIAG, ALCONV,ECONV,OLIMIT,MAXIT,N,M,IVST,NE,ibl 3003 FORMAT(' Average magnitude of matrix el. :',F20.10,/, 1 ' Average diagonal element :',F20.10,/, 2 ' Convergence criterium grad :',F20.10,/, 2 ' Convergence criterium E :',F20.10,/, 3 ' Cut-off norm for subspace :',F20.10,/, 4 ' Maximum number of iterations :',I20,/, 5 ' Dimension of the matrix space :',I20,/, 6 ' Dimension of the iterative space:',I20,/, 8 ' Total eigenvalue :',I20,/, 8 ' Number of vectors to be found :',I20,/, 8 ' Size of Davidson block :',I20) if(LUP)then WRITE(6,4004)ELMAX 4004 format( ' Lower limit for eigenvalues :',F20.10,/,/) else WRITE(6,4001)ELMAX 4001 format( ' Upper limit for eigenvalues :',F20.10,/,/) endif itime0=mclock() OPEN(4,FILE='E.SCR',FORM='UNFORMATTED') NDONE=0 READ(4,END=444)NDONE DO 23 IEE=1,NDONE READ(4)D(IEE) 23 if(ipr.gt.0)write(6,*)D(IEE) close(4) write(6,*)NDONE, ' eigenvalues found on disk' 444 IERM=NDONE+1 if(IVST.GT.IERM)then write(6,6009) 6009 format('Requested starting eigenvalue must be reduced') else IERM=IVST endif WRITE(6,6008)IERM 6008 FORMAT(' Calculation started from eigenvalue number ',I4) C Starting dimension of the subspace: if(M*ibl.gt.NSS)then write(6,*)'M=',M,' ibl=',ibl call report('insufficient subspace dimension') endif c c Loop over eigenvectors: DO 5555 IE=IERM,MIN(N,NE) c c check, that there is still space IF(M+IE-1.GT.N)THEN WRITE(6,*)' Subspace-dimension has to be decreased' M=N-IE+1 IF(M.LT.1)THEN WRITE(6,*)' No subspace for the eigenvector left !' WRITE(6,*)' Try to increase the subspace-dimension.' GOTO 6666 ENDIF ENDIF c c ib1-ib2 interval of evs solved simultaneously ib1=IE ib2=MIN(IE+ibl-1,N,NE) if(ipr.gt.0)write(6,*)'block ',ib1,'-',ib2 c C Estimate starting eigenvalue and eigenvector as unit vector: c for advanced, take previous guess IF(IE.GT.IERM)then do 299 ib=1,ib2-ib1+1-1 EL(ib)=EL(ib+1) do 299 I=1,N 299 YNEW(I+(ib-1)*N)=YNEW(I+(ib+1-1)*N) EL(ib2-ib1+1)=0.0d0 do 2991 I=1,N 2991 YNEW(I+(ib2-ib1+1-1)*N)=AF else do 399 ib=1,ibl 399 EL(ib)=0.0d0 call VZ(YNEW,AF,N*ibl) endif c c C Orthonormalize to previous eigenvectors Y1..YIE-1: DO 51 ib=1,ib2-ib1+1 REWIND 45 DO 71 IN=1,IE-1 READ(45)(Y(I),I=1,N) SN=spv(Y,YNEW,N,0,(ib-1)*N) DO 71 I=1,N 71 YNEW(I+(ib-1)*N)=YNEW(I+(ib-1)*N)-SN*Y(I) call normv(YNEW,1+(ib-1)*N,N+(ib-1)*N) DO 51 I=1,N 51 X(I+(ib-1)*N*M)=YNEW(I+(ib-1)*N) C C Start of the iterative cycle: ITER=0 DELL=0.0d0 9999 ITER=ITER+1 IF(ITER.GT.MAXIT)THEN WRITE(6,*)' Maximum number of iterations exceeded' IERR=8888 GOTO 8888 ENDIF c C Form gradient vectors Xk,k=2..M, new basis: {X1,..,Xk} C Xk+1=AXk-ELXk (Xk,K>1 generated by X1) call vz(AXV,0.0d0,N*M*ibl) DO 3 ib=1,ib2-ib1+1 DO 3 K=2,M rewind 49 K2N=(K-2)*N+(ib-1)*N*M K1N=(K-1)*N+(ib-1)*N*M DO 20 I=1,N c read non-zero components read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) DO 20 ii=1,nn J=j33(ii) AXV(J+K2N)=AXV(J+K2N)+a33(ii)*X(I+K2N) 20 AXV(I+K2N)=AXV(I+K2N)+a33(ii)*X(J+K2N) c C Gradient vectors DO 2031 I=1,N 2031 X(I+K1N)=AXV(I+K2N)-EL(ib)*X(I+K2N) c C Check on convergence on gradient=norm of second trial vector: c of the first block if(ib.eq.1.and.K.EQ.2)then CONV=spv(X,X,N,K1N,K1N) IF(ABS(CONV).LT.ALCONV)GOTO 8888 ENDIF if(IPULAY.eq.1.and.K.eq.M)then DO 2131 I=1,N 2131 if(DABS(DD(I)-EL(ib)).gt.1.0d-6)X(I+K1N)=X(I+K1N)/(DD(I)-EL(ib)) endif c C Orthogonalize vector XK to vectors XIN, IN=1..K-1 (including trial eigv): do 7 ibp=1,ib kend=M if(ibp.eq.ib)kend=K-1 DO 7 IN=1,kend IN1=(IN-1)*N+(ibp-1)*N*M SN=spv(X,X,N,IN1,K1N) DO 7 I=1,N 7 X(I+K1N)=X(I+K1N)-SN*X(I+IN1) c C Orthogonalize vector XK to previous eigenvectors Xin, in=1..IE-1 REWIND 45 DO 72 IN=1,IE-1 READ(45)(Y(I),I=1,N) SN=spv(Y,X,N,0,K1N) DO 72 I=1,N 72 X(I+K1N)=X(I+K1N)-SN*Y(I) c C Normalize vector K: 3 call normv(X,1+K1N,N+K1N) c c make also the matrix product with the last vector do 204 ib=1,ib2-ib1+1 M1N=(M-1)*N+(ib-1)*N*M rewind 49 DO 204 I=1,N read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) DO 204 ii=1,nn J=j33(ii) AXV(I+M1N)=AXV(I+M1N)+a33(ii)*X(J+M1N) 204 AXV(J+M1N)=AXV(J+M1N)+a33(ii)*X(I+M1N) C Now looking for new eigenvector in the subspace of {Xk} C Form the sub-space matrix F=X.A.X: DO 21 ib=1,ib2-ib1+1 DO 21 II=1,M II1N=(II-1)*N+(ib-1)*N*M i1=II+(ib-1)*M DO 21 ibp=1,ib jend=M if(ib.eq.ibp)jend=II DO 21 JJ=1,jend i2=JJ+(ibp-1)*M JJ1N=(JJ-1)*N+(ibp-1)*N*M SUM=spv(AXV,X,N,II1N,JJ1N) F(i2,i1)=SUM 21 F(i1,i2)=SUM c C Diagonalize the sub-space directly: call TRED12(NSS,M,F,Z,E,2,IERR,EWORK) c C New approximation to the eigenvectors: DO 24 ib=1,ib2-ib1+1 if(LUP)then iu=ib else iu=max((ib2-ib1+1)*M-ib+1,1) endif DO 24 I=1,N SUM=0.0d0 DO 13 ibo=1,ib2-ib1+1 DO 13 J=1,M i1=J+(ibo-1)*M 13 SUM=SUM+Z(i1,iu)*X(I+(J-1)*N+(ibo-1)*N*M) 24 YNEW(I+(ib-1)*N)=SUM c call normv(YNEW,1+(ib-1)*N,N+(ib-1)*N) c c Get actual eigenvalue call vz(AXV,0.0d0,N) rewind 49 DO 206 I=1,N read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) DO 206 ii=1,nn J=j33(ii) AXV(I)=AXV(I)+a33(ii)*YNEW(J+(ib-1)*N) 206 AXV(J)=AXV(J)+a33(ii)*YNEW(I+(ib-1)*N) EL(ib)=spv(AXV,YNEW,N,0,(ib-1)*N) wcm=sqrt(abs(EL(ib)))*c1 if(EL(ib).lt.0.0d0)wcm=-wcm 812 if(ipr.gt.0)WRITE(6,6007)ITER,ib,wcm,CONV,M 6007 FORMAT(' It=',I5,' L(',i2,') = ',F20.10,' Grad=',F20.10,' M=',I3) wcm=sqrt(abs(EL(1)))*c1 if(EL(1).lt.0.0d0)wcm=-wcm DELL=ABS(ELAST-wcm) if(DELL.lt.ECONV)goto 8888 ELAST=wcm c c transcript as initial x-vectors: do 250 ib=1,ib2-ib1+1 do 250 I=1,N 250 X(I+(ib-1)*N*M)=YNEW(I+(ib-1)*N) c IF(M.EQ.1)THEN WRITE(6,*)'Last eigenvector left' GOTO 8888 ENDIF c GOTO 9999 8888 wm=dsqrt(dabs(EL(1)))*c1 if(EL(1).lt.0.0d0)wm=-wm WRITE(6,6000)IE,ITER-1,CONV,wm,ALCONV,DELL 6000 Format( 1' *************** The ',I5,'. eigenvalue ******************',/, 2' Iter',I6,' Gradient : ',F20.10,' Lambda : ',F20.10,/, 3' ',6X,' Limit grad.: ',F20.10,' Del Lam.: ',F20.10,/,/) write(6,6001)(mclock()-itime0)/100 6001 format(' Diagonalization time: ',I10,' sec.') MV=IE D(IE)=EL(1) c C Orthonormalize to previous eigenvectors Y1..YIE-1: REWIND 45 DO 81 IN=1,IE-1 READ(45)(Y(I),I=1,N) SN=spv(Y,YNEW,N,0,0) DO 81 I=1,N 81 YNEW(I)=YNEW(I)-SN*Y(I) call normv(YNEW,1,N) c C Save values & vectors into scratch-files: REWIND 45 DO 1003 J=1,IE-1 1003 READ(45) WRITE(45)(YNEW(I),I=1,N) OPEN(4,FILE='E.SCR',FORM='UNFORMATTED') WRITE(4)IE DO 2 IEE=1,IE 2 WRITE(4)D(IEE) CLOSE(4) IF(HW)THEN WRITE(6,*)' The eigenvector:' WRITE(6,3004)(YNEW(J),J=1,N) 3004 FORMAT(7F10.5) ENDIF c j1a=1 y1a=abs(YNEW(1)) do 4 J=2,N if(abs(YNEW(J)).gt.y1a)then j1a=J y1a=abs(YNEW(J)) endif 4 continue j2a=1 y2a=abs(YNEW(1)) if(j2a.eq.j1a)then j2a=2 y2a=abs(YNEW(2)) endif do 31 J=2,N if(abs(YNEW(J)).gt.y2a.and.J.ne.j1a)then j2a=J y2a=abs(YNEW(J)) endif 31 continue write(6,9006)YNEW(j1a),j1a,YNEW(j2a),j2a 9006 format(' Maximum contributions: ',2(F8.4,' of ',i6,'; ')) c if((LUP.and.wm.LT.ELMAX).or.(.NOT.LUP.AND.wm.GT.ELMAX))goto 6666 5555 CONTINUE 6666 M=MSAVE c REWIND 45 c OPEN(4,FILE='E.SCR',FORM='UNFORMATTED') read(4) ii=0 DO 10031 IE=1,MV READ(4)EL1 wcm=sqrt(abs(EL1))*c1 if(EL1.lt.0.0d0)wcm=-wcm READ(45)(YNEW(I),I=1,N) if(IE.ge.ievs1.and.IE.le.ievs2)then ii=ii+1 do 1 I=1,N 1 SMAT(I,ii)=real(YNEW(I)) EMAT(ii)=wcm endif 10031 continue ievs2=ievs1+ii-1 write(6,*)ii,'vectors recorded into S-matrix' CLOSE(4) CLOSE(45) close(49) c write(6,9009)MV 9009 format(I6,' eigenvectors found') RETURN END