PROGRAM HCNvri3 
! Programm-Teil 3 : organize and execute the VRI search
! Programm bestimmt Suchrichtung der Newtontrajektorie, die den
! VRI-Punkt enthält
!
! W.Quapp, Benjamin Schmidt              2009
! 
!     nn = Dimension
!     le = Kettenlänge
!    eps = Genauigkeit des RGF-Verfahrens
!  kette = Kette zwischen den Punkten xstart und xend
!      
!     this part:
! (3) Calculate VRI approximation,
!     one version: use its gradient for new "dir"
!     other version: use only Ag for a predictor step
!     lather parts:
! (4) Organize points for Gamess-US input for (2) or (3)
! HCNvriReadOUT: translates output into usable input files
! (5) Organize Iterations up to convergence of VRI
!     by repeating of (2)-(4)
!
      INTEGER nn,le,N3,I5
      DOUBLE PRECISION eps
      PARAMETER (nn=3,N3=3,le=50,I5=5)
!
! Bestimmung des VRI-nahen Punktes
! "vridir" ist Gradient des VRI-nahen Punktes
!  Verbindungen zwischen allen Punkten von "kette" und VRI
!      
      DOUBLE PRECISION lkette(le+1,nn),kette(le+1,nn),
     1  dir(nn),xstart(nn),xend(nn),vri(nn),
     2  vridir(nn),Agminx,AgMin,xx,
     3  grad(nn),H(nn,nn),Adj(nn,nn),Ag(nn),
     4  vriold(nn),dmnorm,norm,tmin,
     5  gradj(I5+1,nn),Hj(I5+1,nn,nn),
     6  ginvj(I5+1,nn,nn),gmatj(I5+1,nn,nn),
     7  ginv(nn,nn),gmat(nn,nn)
      INTEGER indexJ,istatus,iter,j1,k1,le22,jj
ccccccccccccccccccccccccccccccccccccccccc
      OPEN(8,FILE='hcnVkette.dat')
ccccccccccccccccccccccccccccccccccccccccc
      OPEN(9,FILE='hcnVketNrj.dat')
      OPEN(10,FILE='hcnVket3.dat')
      OPEN(11,FILE='hcnVdir.dat')
      OPEN(14,FILE='hcnVgrad.dat')
      OPEN(15,FILE='hcnVhesse.dat')
      OPEN(16,FILE='hcnVeps.dat')
      OPEN(17,FILE='hcnViter.dat')
      OPEN(18,FILE='hcnVistat.dat')
      OPEN(19,FILE='hcnVicalc.dat')
      OPEN(28,FILE='hcnVketNrj1.dat')
      OPEN(30,FILE='hcnVvri.dat')
      OPEN(33,FILE='hcnVmin.dat')
      OPEN(48,FILE='hcnVgmatr.dat')
c protocol file: all interesting output of the vri-search
      OPEN(44,FILE='hcnVprotocol.txt',status='unknown',access='append')
cccccccccccccccccccccccccccccccccccccccc
       rewind 16
       read(16,*) eps
       j=1
       rewind 9
       write(9,*) j  
       rewind 28
       read(28,*) j1,k1
       rewind 19
       read(19,*) icalc,sgn
       rewind 18       
       read(18,*) istatus
       rewind 17
       read(17,*) iter
       rewind 33
       read(33,*) AgMin, Agminx, xx
c start vri:
       rewind 19
       if(icalc.eq.0) then
        icalc=1
        write(19,*) icalc,sgn
        goto 77
       endif
c loop  vri:
       if(icalc.eq.1) then
        write(44,*)'cycle j1 in vri3  is =', j1
        icalc=2
        rewind 19
        write(19,*) icalc,sgn
        goto 77
       else
        icalc=1
        rewind 19
        write(19,*) icalc,sgn
        goto 88
       endif 
cccc       
77     continue
cccc
       j1=j1+1
       if(j1.eq.le+2) then
        istatus=2
        rewind 18
        write(18,*) istatus
        goto 50
       else 
        istatus=4
        rewind 18
        write(18,*) istatus
       endif
cccc
       write(*,*) '              j1 in  iteration=  (',j1,iter,')'
       write(44,*) '             j1 in  iteration=  (',j1,iter,')'
       rewind 28
       write(28,*) j1,k1
       rewind 8
       Do 23 k=1,le+1
        read(8,*) (lkette(k,i),i=1,nn)
23     continue 
         rewind 30
         read(30,2) (vri(i),i=1,nn)
         read(30,2) (vriold(i),i=1,nn)
        if(j1.eq.1) WRITE(*,*) '   vri      ', (vri(i),i=1,nn)
        if(j1.eq.1) WRITE(*,*) '   vriold   ', (vriold(i),i=1,nn)
c use only xstart on chain, xend will be the former vri point
         DO 10 i=1,nn
           xstart(i)= lkette(j1,i)
cc           xend(i)= lkette(le+1,i)
10       CONTINUE
       if(j1.eq.1) 
     1 WRITE(*,*) '  xstart    ', (xstart(i),i=1,nn)
       DO 173 i=1,nn
       xend(i)  =vriold(i)
173    CONTINUE
2     FORMAT(3F20.14)
      CALL straight_chain(kette,xstart,xend,nn,le)
cccccccccccccccccccccccccccccccccccccccccccccccccccc        
c Use only last I5 nodes of kette up to the VRI
C      I5=5  is given above                    
c Muss in Teil (4) auch eingestellt werden !!  
      rewind 10
      DO 313 k=le-I5,le+1
       WRITE(10,2) (kette(k,i),i=1,nn)
313   CONTINUE
      goto 50  
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
88    continue
cccc
       rewind 30
       read(30,2) (vri(i),i=1,nn)
       read(30,2) (vriold(i),i=1,nn)
cc  I5 siehe oben !!
       rewind 10
       DO 32 jj=1,I5+1
         READ(10,*) (kette(jj,i),i=1,nn)
32     CONTINUE
       rewind 14
       DO 111 jj=1,I5+1
111      read (14,*) (gradj(jj,i),i=1,nn)
       rewind 15
       DO 110 jj=1,I5+1
         Do 112 k=1,nn
112        read (15,*) (Hj(jj,k,i),i=1,nn)
110    continue
       rewind 48
       DO 113 jj=1,I5+1
         Do 766, k=1,nn
766         READ(48,*) (ginvj(jj,k,I),I=1,nn)
            READ(48,*)
            Do 767, k=1,nn
767           READ(48,*) (gmatj(jj,k,I),I=1,nn)
              READ(48,*)
113    continue
cccccccccccccccccccccccccccccccccccccccccccccccc
       indexJ=0  
       DO 101 jj=1,I5 + 1 
c  use current values 
       DO 115 i=1,nn
115       grad(i)= gradj(jj,i)
       Do 116 k=1,nn
       DO 116 i=1,nn 
          gmat(k,i)= gmatj(jj,k,i)
          ginv(k,i)= ginvj(jj,k,i)
          H(k,i)   =    Hj(jj,k,i)
116    continue       
       IF(dmnorm(grad,ginv,nn).GT.0.02d0) THEN
       CALL mat_adj(H,nn,Adj)
       CALL matmult(Adj,nn,grad,1,Ag,nn)
c
c H, g are in bohr, radian (kette is in angstroem, degree)
C Metric norm: , Ag is 1 x contravariant:
       tmin=dmnorm(Ag,gmat,nn) 
       write(44,*) 'tmin with    gmat-norm         ',tmin
       IF(tmin .LT. AgMin ) THEN
           AgMin=tmin
           rewind 33
           write(33,2) tmin, AgMinx, xx
        write(*,*)'          new Ag  -  minimum  = ',tmin 
        write(44,*)'         new Ag  -  minimum  = ',tmin
           do 102 i=1,nn
               vri(i)=kette(jj,i)
102            vridir(i)=gradj(jj,i)
           indexJ=jj  
           ENDIF
        ENDIF
cccc       
101    CONTINUE
cccc
       if( (indexJ.ne.0) ) then
         write(*,*) 'vri   ', (vri(i),i=1,nn)
         WRITE(*,*) 'vriold', (vriold(i),i=1,nn)
         write(*,*) 'vridir', (vridir(i),i=1,nn)
         write(44,*) 'vri   ', (vri(i),i=1,nn)
         WRITE(44,*) 'vriold', (vriold(i),i=1,nn)
         write(44,*) 'vridir', (vridir(i),i=1,nn)
         rewind 30
         write(30,2) (vri(i),i=1,nn)
         write(30,2) (vriold(i),i=1,nn)
         rewind 11
         write(11,2) (vridir(i),i=1,nn) 
       endif
cccccccccccccccccccccccccccccccccccccccccccccccccccc
cc ende jj cycle for any j1  !
cccccccccccccccccccccccccccccccccccccccccccccccccccc
50     CONTINUE
        rewind 18
        write(18,*) istatus
        call exit(istatus)
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccc
!
!  Kette der Länge le zwischen Punkten x und y
!        
      SUBROUTINE straight_chain(kette,x,y,nn,le)
      DOUBLE PRECISION kette(le+1,nn),x(nn),y(nn)
      DO 20 j=1,le+1
        DO 10 i=1,nn
                kette(j,i)=x(i)+((j-1)*(y(i)-x(i)))/le 
10      CONTINUE
20    CONTINUE
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function dmnorm(vec,gm,nn)
c   metric norm of VEC with resp. to metric GM
      real*8 vec,vec_mdot,gm
      integer nn
      dimension vec(nn),gm(nn,nn)
          dmnorm=dsqrt(vec_mdot(vec,gm,vec,nn))
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccc      
! Berechnet adjunkte Matrix A zu M
!     in 3D only ! 
!
      SUBROUTINE adjoint(A,M,nn)
      DOUBLE PRECISION A(nn,nn),M(nn,nn)
      A(1,1)=M(2,2)*M(3,3)-M(2,3)*M(3,2)
      A(1,2)=M(1,3)*M(3,2)-M(1,2)*M(3,3)
      A(1,3)=M(1,2)*M(2,3)-M(1,3)*M(2,2)
      A(2,1)=M(2,3)*M(3,1)-M(2,1)*M(3,3)
      A(2,2)=M(1,1)*M(3,3)-M(1,3)*M(3,1)
      A(2,3)=M(1,3)*M(2,1)-M(1,1)*M(2,3)
      A(3,1)=M(2,1)*M(3,2)-M(2,2)*M(3,1)
      A(3,2)=M(1,2)*M(3,1)-M(1,1)*M(3,2)
      A(3,3)=M(1,1)*M(2,2)-M(1,2)*M(2,1)
      END
! 
      subroutine mat_adj(matin,ndim,matout)
c adjoint matrix of matin
      real*8 matin,matout,matscr,det
      integer ndim,i,j,k,l,ii,jj
      dimension matin(ndim,ndim),matout(ndim,ndim),
     +  matscr(ndim-1,ndim-1)
      call vec_dinit(matscr,(ndim-1)*(ndim-1),0.0d0)
      do 10, i=1 , ndim
        do 20, j=1 , ndim
             do 30, ii=1 , ndim-1
                do 40, jj=1 , ndim-1
                        if (ii.lt.i) then
                         k=ii
                        else
                         k=ii+1
                        endif
                        if (jj.lt.j) then
                         l=jj
                        else
                         l=jj+1
                        endif
                        matscr(ii,jj)=matin(k,l)
c                       if(dabs(matscr(ii,jj)).lt.1.0d-30) then
c                        matscr(ii,jj)=0.0d0
c                       endif
 40               continue
 30             continue
                matout(i,j)=det(matscr,ndim-1)
                if(dabs(matout(i,j)).lt.1.0d-30) goto 20
                matout(i,j)= (-1.0d0)**(i+j)*matout(i,j)
 20     continue
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function det(matein,nn)
c   calculation of the determinant of MATEIN
       integer i,nn
       real*8 matein,mataus,w
       dimension matein(nn,nn),mataus(nn,nn),w(nn)
       call householder(matein,mataus,w,nn,nn)
       det=1.0d0
       do 10,i=1,nn
       if (dabs(mataus(i,i)).lt.1.0d-25) then
          det=0.0d0
          goto 20
       endif
       if (dabs(det) .lt.1.0d-25) then
          det=0.0d0
          goto 20
       endif
 10    det=det*mataus(i,i)
       if (iand(nn,1).eq.0) det=-det
 20    continue
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine householder(matein,mataus,w,n,sp)
c    householder transformation
       integer k,l,m,n,sp
       real*8 alpha,rho,sk,s,w,matein,mataus
       dimension matein(n,sp),mataus(n,sp),w(n)
       call vec_dcopy(matein,mataus,n*sp)
       do 10,k=1,n-1
        s=0.0d0
        do 20,l=k,n
         s=s+mataus(l,k)*mataus(l,k)
 20     continue
        if (mataus(k,k).lt.0) then
         alpha=dsqrt(s)
        else
         alpha=-dsqrt(s)
        endif
        rho=dsqrt(s-mataus(k,k)*mataus(k,k)+
     *       (mataus(k,k)-alpha)*(mataus(k,k)-alpha))
        if(rho .lt. 1.0d-63) then
         rho= 1.0d-63
C            w(k)= 0.0d0
C            goto 22
        endif
        w(k)=(mataus(k,k)-alpha)/rho
 22     continue
        do 30,l=k+1,n
         w(l)=mataus(l,k)/rho
 30     continue
        do 40,l=k,sp
         sk=0.0d0
         do 55,m=k,n
          sk=sk+w(m)*mataus(m,l)
 55      continue
         do 60,m=k,n
           mataus(m,l)=mataus(m,l)-2*sk*w(m)
 60      continue
 40     continue
 10    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccc
!
! Gibt die Norm von v aus
!
      DOUBLE PRECISION FUNCTION norm(v,nn)
      DOUBLE PRECISION v(nn),s
      s=0.d0
      DO 10 j=1,nn
        s=s+v(j)**2
10    CONTINUE
      IF (s.LT.1.D-15) THEN
        norm=0.d0
      ELSE
        norm=DSQRT(s)
      ENDIF
      END
!
! Gibt die Distanz zweier Punkte v und w aus
!
      DOUBLE PRECISION FUNCTION distance(v,w,nn)
      DOUBLE PRECISION v(nn),w(nn),s      
      s=0.d0
      DO 10 j=1,nn
        s=s+(v(j)-w(j))**2
10    CONTINUE
      distance=DSQRT(s)
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine tri_mat(dm1,nr1,nc1,dm2,nc2,dm3,nc3,dres)
c     dm1[nr1,nc1]*dm2[nc1,nc2]*dm3[nc2,nc3]=dres[nr1,nc3]
      real*8 dm1,dm2,dm3,dres,r1
      integer nr1,nc1,nc2,nc3
      dimension dm1(nr1,nc1),dm2(nc1,nc2),dm3(nc2,nc3),
     +          r1(nr1,nc2),dres(nr1,nc3)
      call matmult(dm1,nr1,dm2,nc2,r1,nc1)
      call matmult(r1,nr1,dm3,nc3,dres,nc2)
      return
      end
ccccccccccccccc
       subroutine mat_dcopy(arra1,arra2,ndim,mdim)
c  Copies matrix ARRA1 to matrix ARRA2
       real*8  arra1,arra2
       integer i,j,ndim,mdim
       dimension arra1(ndim,mdim),arra2(ndim,mdim)
       do 10,i = 1,ndim
       do 20,j = 1,mdim
        arra2(i,j) = arra1(i,j)
 20    continue
 10    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function vec_mmdot
     &(nri,matin,din,nra,matout,dout,dim,dmet)
c     multiply the (nri) col of (matin) and the (nra) col of
c     (matout). (dim) is the number of the rows with metric dmet
      real*8 matin,matout,prod,v1,v2,dmet,vec_mdot
      integer nri,din,nra,dout,dim
      dimension matin(dim,din),matout(dim,dout),dmet(dim,dim),
     +   v1(dim),v2(dim)
ccc   if ((nri.gt.din).or.(nra.gt.dout)) call cerr(-100)
      call vec_mat_centry(v1,matin,nri,din,dim,'o')
      call vec_mat_centry(v2,matout,nra,dout,dim,'o')
      prod=vec_mdot(v1,dmet,v2,dim)
      vec_mmdot=prod
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine vec_dcopy(arra1,arra2,ndim)
c  Copies ARRA1 to ARRA2
      real*8 arra1,arra2
      integer i,ndim
      dimension arra1(ndim),arra2(ndim)
      if(ndim .le. 0) return
      do 10 i = 1,ndim
           arra2(i) = arra1(i)
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_mat_centry(vec,mat,npos,ndim,ncol,zkind)
c Copies VEC 'i'n or 'o'ut the NPOS column of matrix MAT
        real*8 vec,mat
        integer npos,ndim,ncol,i
        character zkind
        dimension vec(ndim),mat(ndim,ncol)
        if(npos.gt.ncol) return
        if (zkind.eq.'i') then
         do i=1,ndim
           mat(i,npos)=vec(i)
         end do
        endif
        if (zkind.eq.'o') then
         do i=1,ndim
           vec(i)=mat(i,npos)
         enddo
        endif
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function vec_mdot(vec1,metric,vec2,dim)
c  metric scalar product of VEC1 and VEC2 with respect to METRIC
       integer dim
       real*8 vec1,metric,vec2,cov2,vec_ddot
       dimension vec1(dim),vec2(dim),cov2(dim),metric(dim,dim)
       call matmult(metric,dim,vec2,1,cov2,dim)
       vec_mdot=vec_ddot(vec1,cov2,dim)
       return
       end
ccccccccccccccccccccccccccccccccccccc
       real*8 function vec_ddot(arra1,arra2,ndim)
c   cartesian scalar product of ARRA1 and ARRA2
       real*8  arra1,arra2
       integer ndim,i
       dimension arra1(ndim),arra2(ndim)
       vec_ddot = 0
       if(ndim .le. 0)return
       do 70 i = 1,ndim
         vec_ddot = vec_ddot + arra1(i) * arra2(i)
  70   continue
       return
       end
ccccccccccccccccccccccccccccccccccccccc
       subroutine matmult(m1,d1,m2,d2,mout,dim)
c  matrix multiplication:  m1[d1,dim]*m2[dim,d2]=mout[d1,d2]
       real*8 m1,m2,mout,swap
       integer d1,d2,dim,i,j,k
       dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2)
       call vec_dinit(swap,d1*d2,0.0d0)
       do 10,i=1,d1
        do 20,j=1,d2
         do 30,k=1,dim
           swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j)
 30      continue
 20     continue
 10    continue
       call vec_dcopy(swap,mout,d1*d2)
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_dinit(arra1,ndim,wert)
c Set all elements of ARRA1 to the double value of WERT
       real*8   arra1,wert
       integer  ndim,i
       dimension arra1(ndim)
       if(ndim .le. 0) return
       do i=1,ndim
         arra1(i)=wert
       end do
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      subroutine mat_trans(matin,nrow,ncol,matout)
c     transpose MATOUT of matrix MATIN
      real*8 matin,matout,swap
      integer nrow,ncol,i,j
      dimension matin(nrow,ncol),matout(ncol,nrow),swap(ncol,nrow)
      do 10,i=1,nrow
        do 20,j=1,ncol
                swap(j,i)=matin(i,j)
   20   continue
   10 continue
      call mat_dcopy(swap,matout,ncol,nrow)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_mrowcopy(nri,matin,din,nra,matout,dout,dim)
c  copies the (nri) row of (matin) with (din) row and (dim) col
c  to the (nra) row of (matout) with (dout) row and (dim) col
      real*8 matin,matout
        integer nri,din,nra,dout,dim,i
        dimension matin(din,dim),matout(dout,dim)
ccc      if ((nri.gt.din).or.(nra.gt.dout)) call cerr(17)
      do 10,i=1,dim
        matout(nra,i)=matin(nri,i)
   10 continue
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_madd
     * (nr1,mat1,d1,fact1,nr2,mat2,d2,fact2,nrx,matx,dx,dim)
c     matx(nrx)=fact1*mat1(nr1)+fact2*mat2(nr2)
c     d1,d2,dx - columns
c     dim      - rows
       real*8 mat1,mat2,matx,fact1,fact2,vec
       integer nr1,nr2,nrx,d1,d2,dx,dim,i
       dimension mat1(dim,d1),mat2(dim,d2),matx(dim,dx),vec(dim,1)
       call vec_dinit(vec,dim,0.0d0)
       do 10,i=1,dim
          vec(i,1)=fact1*mat1(i,nr1)+fact2*mat2(i,nr2)
 10    continue
       call vec_mat_centry(vec,matx,nrx,dim,dx,'i')
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_dmult_constant(arra1,ndim,zahl,arra2)
c   ARRA2 = ZAHL * ARRA1
      real*8  arra1,arra2,zahl
      integer ndim,i
      dimension arra1(ndim),arra2(ndim)
      if(ndim .le. 0)return
      do i = 1,ndim
         arra2(i) = arra1(i)*zahl
      end do
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mat_diag(d,n,v)
c create a NxN diagonal matrix D with entries V
       real*8  d,v
       integer n,i
       dimension d(n,n)
              if(n .le. 0) return
               call vec_dinit(d,n*n,0.0d0)
               do i=1,n
                 d(i,i)=v
               end do
               return
              end
c