C CCC CCCCCCCCCCCCCCCCCCCC C W.Quapp 24.10.2007 c Energy, GRADIENT, Hessian, Metric to scratchfile c from GamessUS output file c gmat by svd of ginv CCCCCCCCCCCCCCCCCCCCCCCCCCC Integer NDIM,L,N3 PARAMETER(NDIM=59,L=30,N3=22) Real*8 En,Um,Gr(NDIM),He(NDIM,NDIM),Bmat(NDIM,3*N3), 1 gmat(NDIM,NDIM),ginv(NDIM,NDIM) 2 ,XX(NDIM-1,NDIM-1),YY(NDIM-1,NDIM-1),ZZ(NDIM-1) Integer NAT,I,J,K,Ivar(NDIM),Jnum(6),Istatus,Natoms(N3), 2 mass(3*N3),LA Logical omat_inv Character*5 Z5 Character*6 Z6 Character*7 Z7 Character*2 Z2 Character*18 Z18 Character*80 Z80 Data Natoms/1,6,6,8,1,1,7,6,6,8,1,1,6,1,1,1,7,6,1,1,1,1/ cc OPEN(71,FILE='Energyfile') OPEN(72,FILE='Gradfile') OPEN(73,FILE='Hessfile') OPEN(74,FILE='Bmatfile') OPEN(11,FILE='ala2param.dat') OPEN(14,FILE='ala2ener.dat',status='unknown') OPEN(15,FILE='ala2grad.dat',status='unknown') OPEN(18,FILE='ala2gmatr.dat',status='unknown') OPEN(20,FILE='ala2hesse.dat',status='unknown') OPEN(51,FILE='ala2bmatr.dat',status='unknown') OPEN(44,FILE='ala2protocol.txt',status='unknown',access='append') cc NAT=NDIM istatus=0 Read(11,*) I, LA Write(44,*) ' +++++++++++++++++++++++++++++++++++++++ ' Write(44,*) ' Step LA = ',LA ccccccccccccccccccccccccccccccccccccccccccccccccccccc c Energy read(71,*,END=91) Z5, Z6, Z2, En 81 FORMAT(1X,A18,F18.10) rewind 14 write(14,523) En 523 FORMAT(1X,F18.10) cccccccccccccccccccccccccccccccccccccccccccccccccccccc c Gradient: Do I=1,3 read(72,*) enddo Do I=1,NAT read(72,82,END=92) Ivar(I), Z7, (Jnum(J),J=1,6), En, Gr(I) c write(*,82) Ivar(I), Z7, (Jnum(J),J=1,6), En, Gr(I) enddo 82 FORMAT(1X,I5,1X,A7,1X,6I5,2F16.7) rewind 15 write(15,52) (Gr(I),I=1,NAT) 52 FORMAT(1X, 3F15.11) c cccccccccccccccccccccccccccccccccccccccccccccccccccccc c Hesse matrix Do I=1,10 read(73,*) enddo Do 23 k=1,12 Do I=1,NAT read(73,83,END=93) Ivar(I), (He(I,J),J=(K-1)*5+1,(K-1)*5+5) enddo c Do I=1,3 ! beachte dass hessian output bis 60 geht Do I=1,4 read(73,*) enddo 23 continue 83 FORMAT(1X,I4,1X,5(1X,F11.7)) rewind 20 Do I=1,NAT WRITE(20,53) (He(I,J),J=1,NAT) enddo 53 FORMAT(1X,3F15.10) ccccccccccccccccccccccccccccccccccccccccccccccccccccccc WRITE(44,*) ' svd of hessian ' call SVDCMP(He,NAT,NAT,NAT,NAT,Gr,gmat) Um=10.1d0**10 Do J=1,NAT IF(Gr(J).lt.Um) Um=Gr(J) enddo WRITE(44,*) WRITE(44,*) ' smallest SVD of hessian is ', Um Um=1.0d0 Do J=1,NAT-1 Um=Um*Gr(J) enddo WRITE(44,*) ' Product of SVD of hessian is ', Um WRITE(44,*) ' vector of singular values of hessian ' WRITE(44,50) (Gr(J),J=1,NAT) ccccccccccccccccccccccccccccccccccccccccccccccccccccccc c B matrix Do I=1,5 read(74,*) enddo c 13= [ (Nat+6) : 5 ] Do 22 k=1,13 Do I=1,NAT read(74,84,END=94) Ivar(I), (Bmat(I,J),J=(K-1)*5+1,(K-1)*5+5) write(51,84) Ivar(I), (Bmat(I,J),J=(K-1)*5+1,(K-1)*5+5) enddo read(74,*) read(74,*) ! beachte dim-reduction write(51,*) read(74,80) Z80 write(51,80) Z80 read(74,*) write(51,*) 22 continue 80 FORMAT(A80) Do I=1,NAT read(74,844,END=94) Ivar(I), Bmat(I,66) write(51,844) Ivar(I), Bmat(I,66) enddo write(51,*) write(51,*) ' Ende B Matrix !! ' 84 FORMAT(1X,I4,1X,5(1X,F11.7)) 844 FORMAT(1X,I4,1X,1(1X,F11.7)) cccccc for mass-weighting (not used) Do 221 K=1,22 Do 231 J=1,3 mass((K-1)*3+J)=1 c If(Natoms(K).eq.1) then c mass((K-1)*3+J)=1 c else c mass((K-1)*3+J)=Natoms(K)*2 c endif 231 continue 221 continue c write(44,*) ' ginv is mass-weighted !!! ' cccccc Make Metric G^ij=B*B^T, this is ginv matrix DO 16 I=1,NAT DO 16 J=1,NAT UM=0.0d0 DO K=1,NAT+6 UM=UM+Bmat(I,K)/mass(K)*Bmat(J,K) enddo 16 Ginv(I,J)=UM ccc also dih6 fixed : do 371, I=1,5 do 371, J=1,5 371 XX(I,J)=Ginv(I,J) Do 377, I=6,NAT-1 Do 378, J=6,NAT-1 378 XX(I,J)=Ginv(I+1,J+1) 377 continue rewind 18 Do 37, I=1,NAT-1 37 WRITE(18,50) (XX(I,J),J=1,NAT-1) WRITE(18,*) 50 FORMAT(1X,3F20.10) c goto 1122 ccccccccccccccccccccccccccccccccccccccccccccc c inverse metric calculation WRITE(44,*) ' gmat calculated by ginvers of ginv ' call ginvers(XX,NAT-1,YY,omat_inv) If(.not.omat_inv)write(44,*)' ERROR in gmat calculation ' if(omat_inv) then Do 3381, J=1,NAT-1 3381 WRITE(18,50) (YY(J,I),I=1,NAT-1) endif goto 11 cccccccccccccccccccccccccccccccccccccccccccccc 1122 WRITE(44,*) ' gmat calculated by svd of ginv ' call SVDCMP(XX,NAT-1,NAT-1,NAT-1,NAT-1,ZZ,YY) Um=10.1d0**10 Do J=1,NAT-1 IF(ZZ(J).lt.Um) Um=ZZ(J) enddo WRITE(44,*) WRITE(44,*) ' smallest SVD of ginv(i,j) is ', Um Um=1.0d0 Do J=1,NAT-1 Um=Um*ZZ(J) enddo WRITE(44,*) ' Product of SVD of ginv(i,i) is ', Um WRITE(44,*) ' vector of singular values ' WRITE(44,50) (ZZ(J),J=1,NAT-1) c gmat by svd: write(44,*) ' NO svd threshold for gmat !!! ' Do 55 K=1,NAT-1 ZZ(K)=1.0d0/ZZ(K) 55 continue Do 111 I=1,NAT-1 Do K=1,NAT-1 XX(I,K)=XX(I,K)*ZZ(K) enddo 111 continue c Do 112 I=1,NAT-1 Do 110 J=1,NAT-1 gmat(I,J)=0.0d0 Do 109 K=1,NAT-1 109 gmat(I,J)= gmat(I,J)+XX(I,K)*YY(J,K) 110 continue 112 continue Do 3382, J=1,NAT-1 3382 WRITE(18,50) (gmat(J,I),I=1,NAT-1) WRITE(18,*) goto 11 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc 91 continue write(6,*) 'Stop ohne Energy ' write(44,*) ' ERROR !!! Energy is missing !!! ' goto1 92 continue write(6,*) 'Stop ohne Gradient ' write(44,*) ' ERROR !!! Gradient is missing !!! ' goto1 93 continue write(6,*) 'Stop ohne Hessian ' write(44,*) ' ERROR !!! Hessian is missing !!! ' goto1 94 continue write(6,*) 'Stop ohne Metric ' write(44,*) ' ERROR !!! Metric is missing !!! ' 1 continue istatus=10 WRITE(*,*)' Exit-Status DATEN ', istatus WRITE(44,*)' Exit-Status DATEN ', istatus 11 continue call exit(istatus) STOP End c ------------------------------------------------------------------------ subroutine ginvers(in,NAT,out,omat_inv) c calculates the inverse matrix of IN real*8 in,out logical omat_inv integer NAT,i,j,k dimension in(NAT,NAT),out(NAT,NAT) do 1,k=1,NAT do 1,j=1,NAT 1 out(j,k)=in(j,k) c do 10,k=1,NAT if (out(k,k).eq.0.0d0) then omat_inv=.false. return endif out(k,k)=1.0d0/out(k,k) do 20,i=1,NAT if (i.ne.k) out(i,k)=-(out(i,k)*out(k,k)) 20 continue do 30,j=1,NAT if (j.ne.k) then do 40,i=1,NAT if (i.ne.k) out(i,j)=out(i,j)+out(i,k)*out(k,j) 40 continue out(k,j)=out(k,j)*out(k,k) endif 30 continue 10 continue omat_inv=.true. return end C ccccccccccccccccccccccccccccccccccccccccc SUBROUTINE SVDCMP(A,M,N,Mp,Np,W,V) IMPLICIT NONE INTEGER M, N, Mp, Np DOUBLE PRECISION A(Mp, Np), W(Np), V(Np,Np) c Routine to do singular value decomposition of the matrix A. c Based on the Press et.al routine (Numerical Recipes) INTEGER NMAX PARAMETER (NMAX=100) DOUBLE PRECISION rv1(NMAX) DOUBLE PRECISION anorm, c, f, g, h, s, scale, x, y, z INTEGER i, its, j, k, l, nm, jj, l1 g = 0.0d0 scale = 0.0d0 anorm = 0.0d0 DO i = 1 , N l = i + 1 rv1(i) = scale*g g = 0.0d0 s = 0.0d0 scale = 0.0d0 IF ( i.LE.M ) THEN DO k = i , M scale = scale + ABS(A(k,i)) ENDDO IF ( scale .NE. 0.0d0 ) THEN DO k = i , M A(k,i) = A(k,i)/scale s = s + A(k,i)*A(k,i) ENDDO f = A(i,i) g = -SIGN(SQRT(s),f) h = f*g - s A(i,i) = f - g IF ( i.NE.N ) THEN DO j = l , N s = 0.0d0 DO k = i , M s = s + A(k,i)*A(k,j) ENDDO f = s/h DO k = i , M A(k,j) = A(k,j) + f*A(k,i) ENDDO ENDDO ENDIF DO k = i , M A(k,i) = scale*A(k,i) ENDDO ENDIF ENDIF W(i) = scale*g g = 0.0d0 s = 0.0d0 scale = 0.0d0 IF ( (i.LE.M) .AND. (i.NE.N) ) THEN DO k = l , N scale = scale + ABS(A(i,k)) ENDDO IF ( scale .NE. 0.0d0 ) THEN DO k = l , N A(i,k) = A(i,k)/scale s = s + A(i,k)*A(i,k) ENDDO f = A(i,l) g = -SIGN(SQRT(s),f) h = f*g - s A(i,l) = f - g DO k = l , N rv1(k) = A(i,k)/h ENDDO IF ( i.NE.M ) THEN DO j = l , M s = 0.0d0 DO k = l , N s = s + A(j,k)*A(i,k) ENDDO DO k = l , N A(j,k) = A(j,k) + s*rv1(k) ENDDO ENDDO ENDIF DO k = l , N A(i,k) = scale*A(i,k) ENDDO ENDIF ENDIF anorm = MAX(anorm,(ABS(W(i))+ABS(rv1(i)))) ENDDO DO i = N , 1 , -1 IF ( i.LT.N ) THEN IF ( g .NE. 0.0d0 ) THEN DO j = l , N V(j,i) = (A(i,j)/A(i,l))/g ENDDO DO j = l , N s = 0.0d0 DO k = l , N s = s + A(i,k)*V(k,j) ENDDO DO k = l , N V(k,j) = V(k,j) + s*V(k,i) ENDDO ENDDO ENDIF DO j = l , N V(i,j) = 0.0d0 V(j,i) = 0.0d0 ENDDO ENDIF V(i,i) = 1.0 g = rv1(i) l = i ENDDO DO i = N , 1 , -1 l = i + 1 g = W(i) IF ( i.LT.N ) THEN DO j = l , N A(i,j) = 0.0d0 ENDDO ENDIF IF ( g .NE. 0.0d0 ) THEN g = 1.0/g IF ( i.NE.N ) THEN DO j = l , N s = 0.0d0 DO k = l , M s = s + A(k,i)*A(k,j) ENDDO f = (s/A(i,i))*g DO k = i , M A(k,j) = A(k,j) + f*A(k,i) ENDDO ENDDO ENDIF DO j = i , M A(j,i) = A(j,i)*g ENDDO ELSE DO j = i , M A(j,i) = 0.0d0 ENDDO ENDIF A(i,i) = A(i,i) + 1.0 ENDDO DO k = N , 1 , -1 DO its = 1 , 30 DO l = k , 1 , -1 l1 = l nm = l - 1 IF ( (ABS(rv1(l))+anorm).EQ.anorm ) GOTO 380 IF ( (ABS(W(nm))+anorm).EQ.anorm ) GOTO 340 ENDDO 340 CONTINUE l = l1 c = 0.0d0 s = 1.0d0 DO i = l , k f = s*rv1(i) rv1(i) = c*rv1(i) IF ( (ABS(f)+anorm).NE.anorm ) THEN g = W(i) h = SQRT(f*f+g*g) W(i) = h h = 1.0d0/h c = (g*h) s = -(f*h) DO j = 1 , M y = A(j,nm) z = A(j,i) A(j,nm) = (y*c) + (z*s) A(j,i) = -(y*s) + (z*c) ENDDO ENDIF 360 ENDDO 380 CONTINUE z = W(k) IF ( l.EQ.k ) THEN IF ( z .LT. 0.0d0 ) THEN W(k) = -z DO j = 1 , N V(j,k) = -V(j,k) ENDDO ENDIF GOTO 500 ENDIF IF ( its .EQ. 30 ) write(*,*) & 'SVDCMP: No convergence in 30 iterations' x = W(l) nm = k - 1 y = W(nm) g = rv1(nm) h = rv1(k) f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0d0*h*y) g = SQRT(f*f+1.0d0) f = ((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x c = 1.0d0 s = 1.0d0 DO j = l , nm i = j + 1 g = rv1(i) y = W(i) h = s*g g = c*g z = SQRT(f*f+h*h) rv1(j) = z c = f/z s = h/z f = (x*c) + (g*s) g = -(x*s) + (g*c) h = y*s y = y*c DO jj = 1 , N x = V(jj,j) z = V(jj,i) V(jj,j) = (x*c) + (z*s) V(jj,i) = -(x*s) + (z*c) ENDDO z = SQRT(f*f+h*h) W(j) = z IF ( z .NE. 0.0d0 ) THEN z = 1.0d0/z c = f*z s = h*z ENDIF f = (c*g) + (s*y) x = -(s*g) + (c*y) DO jj = 1 , M y = A(jj,j) z = A(jj,i) A(jj,j) = (y*c) + (z*s) A(jj,i) = -(y*s) + (z*c) ENDDO ENDDO rv1(l) = 0.0d0 rv1(k) = f W(k) = x ENDDO 500 CONTINUE ENDDO RETURN END