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
