        PROGRAM A22ntGS3  
C CCC for Ala potential NDIM=60 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C
C W.Quapp June 2007 for Newton Trajectory by GS following
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
C  W.Q.: JCP 122, iss17, (2005) 174106 (11 pages)
C
c script version for GamessUS (or other PES function)  
c GS cycle for RP: NewPoint calculation for ProjectedGradient=0
c corrector step is done by second order methods of NT theory
c program parts from Michael Hirsch (Dissertation 2003, p.81)
C
C  the old NT point for predictor step is Xold
C  the actual point, to start the correction, is Y(LA,..)
C  put tangent = X-Yold = X - Y(LA-1) 
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc         
      Integer NDIM,L,N3  
      PARAMETER(NDIM=58,L=30,N3=22)
c                    steplength-suche mit 2308 
      Integer LA,N,Nm1
      Integer istatus,ITALL,Natoms(N3)
      REAL*8 energy,gnorm,tnorm,EPS,vec_ddot,D1
c      REAL*8 torad,btoa
      REAL*8 VV(NDIM),told(NDIM),X(NDIM),Proj(NDIM)  
      REAL*8 Y(L+1,NDIM),Yold(NDIM),Xold(NDIM)
      REAL*8 Rs(NDIM),ProMat(NDIM,NDIM),hessian(NDIM,NDIM)
c           ,gold(NDIM) ,hold(NDIM,NDIM)
      Integer JJ,I,J,Np,Np1,IterMax,IFLAG,NFUN
      Character*2 CharAtoms(N3)  
      Character*5 Zcoor(NDIM+2)
c CHASS:
      Data Zcoor/'ch2  ','cc3  ','cch3 ','oc4  ','occ4 ','dih4 ',
     1   'hc5  ','hcc5 ','dih5 ','hc6  ','hcc6 ','dih6 ','nc7  ',
     2   'ncc7 ','dih7 ','cn8  ','cnc8 ','dih8 ','cc9  ','ccn9 ',
     3   'dih9 ','oc10 ','occ10','dih10','hc11 ','hcn11','dih11',
     7   'hn12 ','hnc12','dih12','cc13 ','ccn13','dih13','hc14 ',
     8   'hcc14','dih14','hc15 ','hcc15','dih15','hc16 ','hcc16',
     7   'dih16','nc17 ','ncc17','dih17','cn18 ','cnc18','dih18',
     8   'hc19 ','hcn19','dih19','hn20 ','hnc20','dih20','hc21 ',
     7   'hcn21','dih21','hc22 ','hcn22','dih22'/
      Data CharAtoms/'h','c','c','o','h','h','n','c','c','o','h','h',
     1               'c','h','h','h','n','c','h','h','h','h'/
      Data Natoms/1,6,6,8,1,1,7,6,6,8,1,1,6,1,1,1,7,6,1,1,1,1/
c
       OPEN(9,FILE='ala2chain.dat')  
       OPEN(10,FILE='ala2Rs.dat')
       OPEN(11,FILE='ala2param.dat') 
       OPEN(12,FILE='ala2proj.dat') 
       OPEN(13,FILE='ala2point.dat') 
       OPEN(14,FILE='ala2ener.dat') 
       OPEN(15,FILE='ala2grad.dat')
       OPEN(16,FILE='ala2engs.weg')            
       OPEN(18,FILE='ala2gmatr.dat')
       OPEN(19,FILE='ala2prest.dat')
       OPEN(20,FILE='ala2hesse.dat')
       OPEN(21,FILE='ala2iflag.dat')
       OPEN(25,FILE='ala2told.dat') 
       OPEN(27,FILE='ala2xold.dat')
       OPEN(28,FILE='ala2yold.dat') 
       OPEN(33,FILE='ala2poig.dat')
c       OPEN(34,FILE='ala2gold.dat')
c       OPEN(35,FILE='ala2hold.dat') 
       OPEN(44,FILE='ala2protocol.txt',status='unknown',access='append')
       OPEN(54,FILE='ala2dih917p.txt',status='unknown',access='append')
       OPEN(55,FILE='ala2dih917a.txt',status='unknown',access='append')
c   Read problem input information
       N  = NDIM
       Nm1= NDIM -1 
C      N3 = (N+6)/3
C N3:   number of atoms
       ITERMAX= 10
       istatus= 1
c  (bohr,radian)->(angstroem,degree) for internal coordinates
c      torad=dacos(-1.0d0)/180.0d0
c      btoa =0.52917706d0
c
       read(21,*) IFLAG
       write( 6,*) '        IFLAG is   ',  IFLAG
       write(44,*) '        IFLAG is   ',  IFLAG
c
c     IFLAG=0 indicates an initial entry into corrector
c
       rewind 11
       read(11,*) JJ, LA, eps ,ITALL, NFUN
cc       JJ=L  is constant   !!
       rewind 11
       write( *,*) L, LA, eps ,ITALL, NFUN
       write(11,*) L, LA, eps ,ITALL, NFUN
       write(44,*)' L=',L,' LA= ',LA,' ITall = ',ITALL,'+',NFUN
       IF(LA.eq.L+1) then
        istatus=10
        goto 999
       Endif 
c  Print parameters
       if ( LA.lt.3 .AND. IFLAG.eq.0)then
         write(*,820)
         write(*,840) N
         write(44,820)
         write(44,840) N
       endif
c current point X
       rewind 13
       READ(13,*) (X(I),I=1,N)
        write(44,*) ' X is  '
        write(44,50) (X(I),I=1,5)
        write(44,50) (X(I),I=6,N)
c string Y 
       rewind 9
       Do 4, J=1,L+1
       Read(9,50) (Y(J,I),I=1,N)
4      Read(9,50)
c current string Y goes up to LA
       if(IFLAG.eq.0)  then
         do 45 I=1,N
45       Yold(I) = Y(LA-1,I)
         rewind 28
         write(28,50) (Yold(I),I=1,N)
c         write(*,*) ' Yold is '
c         write(*,50) (Yold(I),I=1,N)
         Do 211 I=1,N
211      Xold(I)= Y(LA,I)
         write(44,*) 'IFLAG is ', IFLAG 
         rewind 27
         write(27,50) (Xold(I),I=1,N)
       else
         read(27,50) (Xold(J),J=1,N)
         read(28,50) (Yold(J),J=1,N)
c         read(34,50) (gold(J),J=1,N)
c         write(44,*) ' Xold is '
c         write(44,50) (Xold(I),I=1,N)
       endif 
c
        IF((LA.lt.3) .and. (IFLAG.eq.0)) then
         do 451 I=1,N
451      told(I) = X(I)- Yold(I)
c (bohr,radian)<<---(angstroem,degree) for internal coordinates
         call angtobohr(told,N)
         gnorm=DSQRT(vec_ddot(told,told,N))
         call vec_dmult_constant(told,N,1.0d0/gnorm,told)
         write(44,*) ' first told is '
         write(44,50) (told(I),I=1,5)
         write(44,50) (told(I),I=6,N)
        else
         read(25,50) (told(J),J=1,N)
         write(44,*) ' told is read '
         write(44,50) (told(I),I=1,5)
         write(44,50) (told(I),I=6,N)
        endif
cc 
        If(NFUN .gt. Itermax) then
         write(6,*) ' Cancel Loop: ERROR , no Convergence at LA',LA
         write(6,*) ' LA=  ' ,LA, ' NFUN  =  ', NFUN
         write(44,*)' Cancel Loop: ERROR , no Convergence at LA ',LA
         istatus=1
         IFLAG=2
         ITALL = ITALL+NFUN
         NFUN=0 
         write(11,*) L, LA, eps ,ITALL, NFUN
         write(44,*)' LA= ',LA,' ITall= ',ITALL,' NFUN= ',NFUN
         write(44,*)
         write(44,*) ' STOP corrector  '
         goto 111
        Endif
c read the SearchDirection Rs and (old) global ProjectionMatrix here
       rewind 10
       READ(10,50) (Rs(I),I=1,N)
       if(IFLAG.eq.0) WRITE(44,*) 'Rs '
       if(IFLAG.eq.0) WRITE(44,50) (Rs(I),I=1,5)
       if(IFLAG.eq.0) WRITE(44,50) (Rs(I),I=6,N)
       gnorm=vec_ddot(Rs,told,N)
       WRITE(44,*)'               Cartesian  Rs * told =', gnorm
c
       rewind 12
       Do 3 J=1,N
       READ(12,50) (ProMat(J,I),I=1,N)
cc       if(IFLAG.eq.0) WRITE(44,*) ' J= ',J
cc       if(IFLAG.eq.0) WRITE(44,50) (ProMat(J,I),I=1,N)
3      READ(12,50) 
c read function,  chain, and  gradient values here
       rewind 14
       READ(14,*)  energy
       WRITE(44,*)'energy=',energy 
       WRITE(*,*) 'energy=',energy
c file with Grad
       rewind 15
       READ(15,*) (VV(I),I=1,6)
c note: override first dih4 of first CH3 group troughout !! 
       READ(15,*) (VV(I),I=6,N)
       write(44,*) ' gradient is '
       write(44,50) (VV(I),I=1,5)
       write(44,50) (VV(I),I=6,N)
c  cccc   IF( NFUN .eq.0 ) then
          rewind 20
          Do 120, J=1,6
          READ(20,52) ( hessian(J,I),I=1,6 )
120       READ(20,52) ( hessian(J,I),I=6,N )
          Do 121, J=6,N
          READ(20,52) ( hessian(J,I),I=1,6 )
121       READ(20,52) ( hessian(J,I),I=6,N )
c test
          write(44,*) '  Diag of Hessian is '
          write(44,50) (hessian(I,I),I=1,5)
          write(44,50) (hessian(I,I),I=6,N)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Hess Update 
c          call vec_dcopy(hessian,hold,N*N)
c         write(44,*) ' === start analytical Hessian === '
c  cccc       ENDIF
c        IF( NFUN.gt.0 ) then
c          rewind 35 
c          Do 122, J=1,N
c122       READ(35,50) (hold(J,I),I=1,N)
c          write(44,*) ' === use updated Hessian from former LA === '
c        ENDIF
c         write(44,*) ' NFUN = ', NFUN
c         write(44,*) ' hessian matrix befor Corr'
c         Do 7770, J=1,N
c7770      WRITE(44,50) (hold(J,I),I=1,N)
c
c    Cartesian  norm  of  gradient
        gnorm= DSQRT(vec_DDOT(VV,VV,N))
        WRITE(*,*) ' PES Cartesian grad norm ',gnorm
        WRITE(44,*)' PES Cartesian grad norm ',gnorm
ccccccccccccccccccccccccccccccccccccccccccALT:
C   Newton trajectory with Projector to Rs   
C   Old: do the "Cartesian" Projektion of Gradient    
       call vec_dcopy(VV,Proj,N)
       call PROJECTION(N,Proj,ProMat)
       D1= DSQRT(vec_DDOT(Proj,Proj,N))
c     WRITE(6,*) ' projected   grad   norm ', D1
      WRITE(44,*)' projected Cart. grad norm', D1
      If(D1 .lt. 1.0d0/30.0d0*eps) THEN
        WRITE(44,*)'at IteGS=',NFUN,' ok Cart proj norm less then', 
     1  ' 1/30 * tolerance EPS ', eps/30.0d0 
        WRITE(44,*)' successful convergence (no ERROR) '
        WRITE(6, *)' IteGS=',NFUN,' == ok - Cart proj norm',D1
        ITALL = ITALL+NFUN
        NFUN=0 
        IFLAG=2 
        istatus=1
        GOTO 111
      ENDIF
c
c     Call the new Newton corrector code
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      CALL newtcorr(NDIM,IFLAG,LA,Eps,X,VV,hessian,Rs,Xold,told)     
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       IFLAG=
c              1 : Corrector step inside the corr - loop
c              2 : return to predictor
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        NFUN=NFUN+1 
        rewind 11
        write(11,*) L, LA, eps ,ITALL, NFUN
c actual point approximation:
        rewind 13
        WRITE(13,50) (X(I),I=1,N)
ccc   for gamessUS input file
        rewind 33       
        Do 7, J=1,5
 7      WRITE(33,*) Zcoor(J), X(J)
        Do  J=6,N
        WRITE(33,*) Zcoor(J+1), X(J)
        enddo
        WRITE(33,*) '  ' 
        rewind 21
        write(21,*) IFLAG
c        write(44,*)' IFLAG after corr ',  IFLAG
        write(44,*)
        rewind 25
        write(25,50) (told(J),J=1,N)
c
c         call vec_dcopy(VV,gold,N)
c         rewind 35
c         Do 77, J=1,N
c77       WRITE(35,50) (hold(J,I),I=1,N)
c         write(44,*) ' hessian-matrix nach Corr'
c         Do 770, J=1,N
c770      WRITE(44,50) (hold(J,I),I=1,N)
c         Do 771, J=1,N
c771      WRITE(20,50) (hold(J,I),I=1,N)
c         WRITE(20,50) 
c         rewind 34
c         write(34,50) (gold(J),J=1,N)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      write (*,*) ' energy   ', energy 
      IF(IFLAG.EQ.2) THEN
         rewind 27
         write(27,50) (Xold(J),J=1,N)
         write(44,890) energy
         write(44,*)
         Do 17, J=1,5
17       WRITE(44,*) Zcoor(J),'  ',X(J)
         Do  J=6,N
         WRITE(44,*) Zcoor(J+1),'  ',X(J)
         enddo 
c output for Mma Figure:
        Do J=1,N 
         If(Zcoor(J).eq.'dih9') WRITE(54,*) '{  ',X(J-1),' ,'
         If(Zcoor(J).eq.'dih17') WRITE(54,*) X(J-1),' }, '
         If(Zcoor(J).eq.'dih9') WRITE(55,*) '{  ',X(J-1),' ,'
         If(Zcoor(J).eq.'dih17') WRITE(55,*) X(J-1),' }, '
        Enddo 
         IF(LA.le.L) write(44,*)' new step of NewtonTrajectory: '
      ENDIF
ccccccccccccccc
111   continue
ccccccccccccccc
      rewind 11
      write(11,*) L, LA, eps ,ITALL, NFUN
      rewind 21
      write(21,*) IFLAG
c   Istatus controls the corrector loop 
      IF(IFLAG.EQ.2) THEN
        istatus=1
      ELSE
        istatus=0 
      ENDIF
c  new chain point becomes:
      Do 117 I=1,N
117   Y(LA,I)= X(I)
      rewind 9
      Do 118, J=1,L+1
       write(9,50) (Y(J,I),I=1,N)
118    write(9,50)
      WRITE(16,*) '{ ' 
      IF(LA.le.L)  WRITE(16,*) '  ',  energy,' , '
      IF(LA.eq.L+1)WRITE(16,*) '  ',  energy,' } '
999   continue 
      close(44)
      call exit(istatus)
      STOP        
 50    FORMAT(1X,3F20.10)
 51    FORMAT(1X, F25.16)
 52    FORMAT(1X,3F15.10)
 820  format ( ' *** NT Routine for ProjGrad=zero *** ')
 840  format ( ' Ndim   =',i6 )
 890  format (/,' Energy at pes(x*) = ', 1pd16.8)
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
       SUBROUTINE PROJECTION(N,VEC,ProMat)
c  projection of vector VEC of dim N 
c  by Projection Matrix ProMat of dim(N,N)
c   ! old version in Cartesian metric !
      REAL*8 VEC(N),ProMat(N,N),Proj(N)
      INTEGER N,J,JJ
      Do 1 J=1,N
      Proj(J)=0.0d0
      Do 2 JJ=1,N   
2     Proj(J)=Proj(J) + Vec(JJ)*ProMat(J,JJ)
1     Continue
      Do 3 J=1,N 
3     Vec(J) = Proj(J)     
      RETURN
      END  
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
      subroutine newtcorr(nn,IFLAG,LA,EPS,point,grad,hessian,rsel,
     *  xold,told)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC W.Quapp, February 2008
CC for GS-NT corrector by NT second order method
CC with parts of a program written by Michael Hirsch
CC use metric of internal z matrix coordinates 
CC
CC Pred: replace tangent of NT methods by predictor step of GS
CC ( next program part ***5.f )
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      integer i,j,nn,n3
      real*8 EPS,EPSneu,hessian,gradcont,
c            torad,btoa,dsqrt,gold,
     1 rsel,dk,gmat,unit,rgeq,rgrad,point,vec_mdot,
     2 grad,tang,tcov,told,gnorm,dnst,tgt,tg,dmnorm
      real*8 xold,projNm1,tt,vec_ddot,dstep,rgnorm,ginv,pla,
     3       prestep  
      dimension rsel(nn),dk(nn,nn+1),hessian(nn,nn),grad(nn),
     1 rgeq(nn-1,nn),dstep(nn),projNm1(nn-1,nn),
     2 rgrad(nn-1),point(nn),xold(nn),gradcont(nn),
c       gold(nn),hold(nn,nn),dx(nn),dg(nn),
     3  tcov(nn),tang(nn),told(nn),gmat(nn,nn),unit(nn,nn),
     4  ginv(nn,nn),prestep(nn)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c        write(44,*) '( point ) in angstroem/degree in Corrector '
c        write(44,50) point(1)
c        write(44,50) point(2), point(3)
c        write(44,50) (point(I),I=4,NN)

        n3=(nn+ 6 +2 )/3
         
c test holds only for full dimensional case:
c         tt=(point(2)**2)*Dsin(point(3))
c         Do 327 K=2,NAT/3
c         tt=tt*(point(3*(K-1)+1)**2)*Dsin(point(3*(K-1)+2))
c 327     continue
c         write(44,*) 'theor.Sqrt of |g_ij| in angstroem ',tt
c         tt=tt*(1.0d0/0.52917706d0)**21
c         write(44,*) 'theor.Sqrt of |g_ij| in bohr      ',tt
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c (bohr,radian)<<--(angstroem,degree) for internal coordinates
c note: angtobohr depends from the z-matrix used !! 
c       it is to adapt at every change           !!
c       torad=dacos(-1.0d0)/180.0d0
c       btoa =0.52917706d0
        call angtobohr(point,NN)
        write(44,*) '( point ) in bohr,radian '
        write(44,50) (point(I),I=1,5)
        write(44,50) (point(I),I=6,NN)
        write(44,*)
        call mat_diag(unit,nn,1.0d0)
c        
c note: gmat and ginv interchanged agains use in Vibrational Theory
c here: mathematical gmat is metric g_ij, 
c g^ij is its contravariate inverse
c now gmat is obtained by svd of ginv, ginv is from the B-matrix
c 
        rewind 18
        Do  J=1,6
        READ(18,50) (ginv(J,I),I=1,6)
        READ(18,50) (ginv(J,I),I=7,nn)
        enddo
        Do  J=7,nn
        READ(18,50) (ginv(J,I),I=1,6)
        READ(18,50) (ginv(J,I),I=7,nn)
        enddo
c
        READ(18,*)
        Do 767, J=1,6
        READ(18,50) (gmat(J,I),I=1,6)
767     READ(18,50) (gmat(J,I),I=7,nn)
        Do  J=7,nn
        READ(18,50) (gmat(J,I),I=1,6)
        READ(18,50) (gmat(J,I),I=7,nn)
        enddo
cccccccc
        write(44,*) '  ginv I,I is '
        WRITE(44,50) (ginv(I,I),I=1,5)
        WRITE(44,50) (ginv(I,I),I=6,nn)
        write(44,*) 
c        write(44,*) '  gmat(21,I) is dih 9=phi '
c        WRITE(44,50) (gmat(21,I),I=1,nn)
c        write(44,*) '  ginv(45,I) is dih 17=psi'
c        WRITE(44,50) (ginv(45,I),I=1,nn)
c        write(44,*)
        write(44,*) '  gmat I,I is '
        WRITE(44,50) (gmat(I,I),I=1,5)
        WRITE(44,50) (gmat(I,I),I=6,nn)
c unit test 
        call matmult(ginv,nn,gmat,nn,unit,nn)
        write(44,*) '  ginv*gmat I,I is '
        WRITE(44,50) (unit(I,I),I=1,nn)
        write(44,*) '  ginv*gmat 17,I is '
        WRITE(44,50) (unit(17,I),I=1,nn)
        gnorm=dmnorm(grad,ginv,nn)
        WRITE(44,*) 'PES full grad-metric norm',gnorm
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc  test: it comes out successful
c        write(44,*)' only  DIAG METRIC is used in ginv'        
c        Do 7661, J=1,nn                                        
c        Do 7661, I=1,nn  
c7661    If (I.ne.J) ginv(J,I)=0.0d0                            
c        call mat_diag(gmat,nn,1.0d0)                           
c        write(44,*)'DIAG METRIC 1/diag(ginv) is used in gmat'  
c        Do 7671, J=1,nn                                        
c7671    gmat(J,J)=1.0d0/ginv(J,J)                              
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C   contravariant gradient
        call vec_dinit(gradcont,nn,0.0d0)
        call matmult(ginv,nn,grad,1,gradcont,nn)
c        write(44,*) ' contra gradient is '
c        WRITE(44,50) (gradcont(I),I=1,nn)
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Version of old rgf: direction of rsel with first gradient...
c       if (istep.eq.0 .and. .... ) then
c        call matmult(ginv,nn,gradient,1,rsel,nn)
c        call ortcompl(rsel,proj,unit,nn)
c       endif
c Note: rsel not used in metric form, up to now !!
       tt=dmnorm(rsel,gmat,nn)
       WRITE(44,*)'metric norm of rsel befor normalization ',tt
       call vec_dmult_constant(rsel,nn,1.0d0/tt,rsel)
c
C  rsel search direction to (N-1)xN - ProjMat 
cc test: gmat to make projNm1 lines contravariate
cc       call ortcompl(rsel,projNm1,gmat,nn)
ccc  but use again old version by MH:
       call ortcompl(rsel,projNm1,unit,nn)
c           write(44,*)'  projNm1 is   '
c           do 121 J=1,nn-1
c           write(44,50) (projNm1(J,I),I=1,nn)
c121        write(44,50) 
c  then                        here  unit is correct for Pr * H  
cc      call tri_mat(projNm1,nn-1,nn,unit,nn,hessian,nn,rgeq)
ccc  but use again old version by MH:
        call tri_mat(projNm1,nn-1,nn,ginv,nn,hessian,nn,rgeq)
        write(44,*)'  rgeq(I,I) is   '
        write(44,50) (rgeq(I,I),I=1,5)
        write(44,50) (rgeq(I,I),I=6,nn-1)
        call vec_dinit(rgrad,nn-1,0.0d0)
cc  projNm1 is contravariate, thus it is to apply to grad
c   gradcont is useless now        
cc      call matmult(projNm1,nn-1,grad,1,rgrad,nn)
ccc  but use again old version by MH:
       call matmult(projNm1,nn-1,gradcont,1,rgrad,nn)
       IF(IFLAG.eq.0) then 
          write(44,*)' rgrad is  '
          write(44,50) (rgrad(I),I=1,nn-1)
       ENDIF
       rgnorm=dsqrt(vec_ddot(rgrad,rgrad,nn-1))
       WRITE(6,*) ' projected rgnorm        ',rgnorm
       WRITE(44,*)' projected rgnorm        ',rgnorm
c Tangent : gmat is used for normalization of the tangent
c which is contravariate
         call gettang(rgeq,gmat,tang,nn)
         write(44,*) ' tang is '
         write(44,50) (tang(I),I=1,5)
         write(44,50) (tang(I),I=6,nn)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccctest         
c new (calculation without bifurcation test)
         tgt=vec_mdot(tang,gmat,told,nn)
         write(44,*) ' tang*gmat*told is ', tgt
      if (tgt.lt.0.0d0) call vec_dmult_constant(tang,nn,-1.0d0,tang)
         call vec_dcopy(tang,told,nn) 
         call vec_dinit(tcov,nn,0.0d0)
         call matmult(gmat,nn,tang,1,tcov,nn)
         write(44,*) ' tcov is '
         write(44,50) (tcov(I),I=1,5)
         write(44,50) (tcov(I),I=6,nn)
         pla=0.0d0
cccccccccccccccccccccccccccccccccccccccccccccc 
c        call buildk(dk,rgeq,tang,rgrad,nn,pla)
ccc  but use again old version by MH:
         call buildk(dk,rgeq,tcov,rgrad,nn,pla)
cccccccccccccccccccccccccccccccccccccccccccccc
c  t1 for bifurcation detection in old program:
c         detk=det(dk,nn)
c         write(44,*) '        detK is         ', detk
c        if (detk.lt.0.0d0) then
c           call vec_dmult_constant(tang,nn,-1.0d0,t1)
c        else
c           call vec_dcopy(tang,t1,nn)
c        endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c         tgt=vec_mdot(tang,gmat,told,nn)
c         write(44,*) ' tang*gmat*told is  ', tgt
c         tgt=vec_ddot(tang,told,nn)
c         write(44,*) ' tang 1 *  told is  ', tgt
cc        if (tgt.lt.0.0d0) then
cc           call vec_dmult_constant(tang,nn,-1.0d0,tang)
cc           call vec_dmult_constant(tcov,nn,-1.0d0,tcov)
cc           call buildk(dk,rgeq,tcov,rgrad,nn,pla)
cc         endif
cc         call vec_dcopy(tang,told,nn)
ccccccccccccccccccccccccccccccccccccccccccccc
c        tighter near stationary points the eps:
c         EPSneu=EPS
c         If(Eps.gt.0.1d0*gnorm .and. LA.gt.5) then
c           EPSneu=0.1d0* gnorm
c           write(44,*) 'new eps is 0.1*gnorm',EPSneu
c         else
c           write(44,*) ' eps is           ', EPS
c         Endif
cccccccccccccccccccccccccccccccccccccccccccccc
       if (rgnorm.gt.EPS) then
c  Corrector
        call linsolve(dk,dstep,nn)
         tt=dmnorm(dstep,gmat,nn)
         write(*,*) 'test dstep norm tdn ',tt
         if(tt .gt. 1.0d0) tt=1.0d0
         read(19,50) (prestep(I),I=1,nn)
         call angtobohr(prestep,nn)
         pla=vec_ddot(prestep,tcov,nn)
         write(*,*)  'prestep*tcov ',pla
         tgt=dmnorm(prestep,gmat,nn)
         write(*,*)  ' |predstep|  ',tgt
         pla=pla/tgt
         tgt=dmnorm(tcov,ginv,nn)
         write(*,*)  ' norm tcov   ',tgt
         pla=pla/tgt
         write(*,*)  'pla vor delta',pla
         pla=Dsqrt(Dabs(1.0d0-pla**2))*tt*0.90d0 
         write(*,*)  '  sin(alpha)*tdn*0.9  =  ',pla 
c        If(pla .gt. 0.25d0) pla=0.25d0
         write(44,*) ' delta von pred-step    ', pla
         call buildk(dk,rgeq,tcov,rgrad,nn,pla)
         call linsolve(dk,dstep,nn)
ccccccccccccccccccccccccccccccccccccccccccccccccccccc                                                                       
         IFLAG=1
         goto 111
       else
c   goto Predictor
        IFLAG=2
        write(44,*)
        write(44,*)'Convergence (no ERROR) - go to next Pred, LA=',LA
c (bohr,radian)-->>(angstroem,degree)
        call bohrtoang(point,nn)
        call vec_dcopy(point,xold,nn)
       endif
       return
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c  normalization of correctors dstep
111     continue
        tt=dmnorm(dstep,gmat,nn)
        ttr=tt 
        write(44,*) ' step-norm of corrector ', tt
        stepmax=1.000000000000000000002308d0
        if(tt .lt. stepmax) then
         tt=1.0d0 
        else
         write(44,*)' ! ERROR corr-step shorten to ',stepmax
         tt=stepmax/tt
        endif
c execute step
       call vec_dmult_add(point,dstep,nn,tt,point)
       write(44,*) ' dstep is '
       write(44,50) (tt*dstep(I),I=1,5)
       write(44,50) (tt*dstep(I),I=6,nn)
cc TestStep to next stationary point
       call newton(hessian,grad,dstep,dnst,gmat,NN)
       write(44,*) ' dNewtonStep norm whould be: ', dnst
c(bohr,radian)-->>(angstroem,degree)
        call bohrtoang(point,NN)
        call vec_dcopy(tang,told,nn)
        return 
50      FORMAT(1X,3F20.10)
        end       
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine buildk(dk,req,tang,rgr,nn,pla)
c build the K-Matrix
      real*8 dk,req,tang,rgr,pla
      integer nn,i,j
      dimension dk(nn,nn+1),req(nn-1,nn),tang(nn),rgr(nn-1)
      do 10,i=1,nn-1
        do 11,j=1,nn
        dk(i,j)=req(i,j)
 11   continue
 10   continue
      do 20,j=1,nn
         dk(nn,j)=tang(j)
 20   continue
      do 30,i=1,nn-1
         dk(i,nn+1)=-rgr(i)
 30   continue
ccccc eventuell: ersetze pla=zero durch kleinen Schritt
      dk(nn,nn+1)=pla
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine newton(hi,fi,step,stepl,gmat,nn)
c  calculates STEP = - HI^(-1) * FI and returns step and steplength
       real*8 hi,fi,step,stepl,gmat,gl,vec_mdot
       integer i,nn
       dimension hi(nn,nn),fi(nn),step(nn),gmat(nn,nn),gl(nn,nn+1)
       call mat_dcopy(hi,gl,nn,nn)
       do 10,i=1,nn
        gl(i,nn+1)=-fi(i)
  10   continue
       call linsolve(gl,step,nn)
       stepl=dsqrt(vec_mdot(step,gmat,step,nn))
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc       
       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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       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                    
       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
c         ccccccccccccccccccccccccccccccccccccccccccccccccc                                                              
       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
c                                                                                                           
       subroutine vec_mat_rentry(vec,dmat,npos,ndim,nrow,zkind)
c Copies VEC 'i'n or 'o'ut the NPOS row of matrix MAT
       real*8 vec,dmat
       integer npos,ndim,nrow,i
       character zkind
       dimension vec(ndim),dmat(nrow,ndim)
       if(npos.gt.nrow) return
       if (zkind.eq.'i') then
        do 10,i=1,ndim
 10      dmat(npos,i)=vec(i)
       endif
       if (zkind.eq.'o') then
        do 20,i=1,ndim
 20     vec(i)=dmat(npos,i)
       endif
       return
       end
c
       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
c
       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
c
       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
c                                                                       
       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
c       
       subroutine linsolve(matein,x,n)
c  solves MATEIN * X = 0
       real*8 matein,dre,x,w
       integer n
       dimension matein(n,n+1),dre(n,n+1),x(n),w(n)
       call householder(matein,dre,w,n,n+1)
       call loesung(dre,x,n)
       return
       end
c  
       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
c       
       subroutine loesung(matein,x,n)
c   intern linear equation routine
       real*8 matein,x
       integer i,j,n,k
       dimension matein(n,n+1),x(n)
       do 10,i=1,n
 10     x(i)=matein(i,n+1)
       x(n)=x(n)/matein(n,n)
       do 20,i=1,n-1
         k=n-i
         do 30,j=k+1,n
 30        x(k)=x(k)-matein(k,j)*x(j)
           x(k)=x(k)/matein(k,k)
 20    continue
       return
       end
c      
      subroutine dfp(einmat,q,p,nn,upmat)
ccc   DFP-Update of Hessian
ccc Heidrich,Kliesch,Quapp: Properties...
ccc S.51 Eqn.(21)/(22)
ccc
ccc einmat        Hessian to update      (H_0)
ccc q             difference of gradient (g_1 - g_0)
ccc p             difference of point    (x_1 - x_0)
ccc upmat         updated Hessian        (H_1)
       real*8 einmat,q,p,upmat,qp,unit,b1,b2,b3,
     1     vec_ddot,eps,sw
       integer nn
       dimension unit(nn,nn),b3(nn,nn),
     +  einmat(nn,nn),upmat(nn,nn),
     +  b1(nn,nn),b2(nn,nn),sw(nn,nn),
     +  q(nn), p(nn)
        data eps/1.0d-30/
        call vec_dinit(upmat,nn*nn,0.0d0)
        call mat_diag(unit,nn,1.0d0)
        qp=vec_ddot(q,p,nn)
ccc        if(dabs(qp).lt.eps) call cerr(-31)
        call matmult(q,nn,p,nn,sw,1)
        call matadd(unit,(-1.0d0)/qp,sw,b1,nn,nn)
        call matmult(p,nn,q,nn,sw,1)
        call matadd(unit,(-1.0d0)/qp,sw,b2,nn,nn)
        call matmult(q,nn,q,nn,b3,1)
        call tri_mat(b1,nn,nn,einmat,nn,b2,nn,sw)
        call matadd(sw,1.0d0/qp,b3,upmat,nn,nn)
        return
        end
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
c
       subroutine ortcompl(vec,compl,dmet,ndim)
c   calculate an Orthonormalcomplement COMPL to vector VEC
       real*8 basis,vec,compl,swap,det,dmgs,dmet
       integer i,j,ndim
       dimension vec(ndim),compl(ndim-1,ndim),basis(ndim,ndim),
     *   swap(ndim),dmet(ndim,ndim)
       call vec_dinit(basis,ndim*ndim,0.0d0)
       do 10,i=1,ndim
         basis(i,i)=1.0d0
 10    continue
       do 20,i=1,ndim
         if (vec(i).ne.0.0d0) then
           do 30,j=1,ndim
             swap(j)=basis(j,1)
             basis(j,1)=vec(j)
             if (i.ne.1) basis(j,i)=swap(j)
 30        continue
           goto 40
         endif
 20    continue
 40    continue
       det=dmgs(basis,dmet,ndim)
       call vec_dinit(compl,(ndim-1)*ndim,0.0d0)
       call mat_trans(basis,ndim,ndim,basis)
       do 60,i=2,ndim
          call vec_mrowcopy(i,basis,ndim,i-1,compl,ndim-1,ndim)
 60    continue
       return
       end
c               
       real*8 function dmgs(a,dmet,dim)
c  modified Gram-Schmidt algorithm return det(a)
       real*8 det,rr,a,q,r,vec_mmdot,dmet,dzero
       integer dim,i,k
       dimension a(dim,dim),q(dim,dim),r(dim,dim),dmet(dim,dim)
       data dzero/0.0d0/
       call vec_dinit(q,dim*dim,0.0d0)
       call vec_dinit(r,dim*dim,0.0d0)
       det=1.0d0
       call mat_dcopy(a,q,dim,dim)
       do 10,k=1,dim
         if (k.eq.1) goto 21
         do 20,i=1,k-1
            r(i,k)=vec_mmdot(i,q,dim,k,q,dim,dim,dmet)
            call vec_madd(k,q,dim,1.0d0,i,q,dim,-r(i,k),k,q,dim,dim)
 20      continue
 21      continue
         rr=vec_mmdot(k,q,dim,k,q,dim,dim,dmet)
         r(k,k)=dsqrt(rr)
         det=det*r(k,k)
         if(rr.lt.dzero) goto 33
         call vec_madd(k,q,dim,1.0d0/r(k,k),k,q,dim,0.0d0,k,q,dim,dim)
 10    continue
       call mat_dcopy(q,a,dim,dim)
       dmgs=det
       return
 33    CONTINUE
cc       call cerr(-23)
       return
       end
c       
       subroutine bof(einmat,q,p,nn,upmat)
cc      Bofill Update of Hessian, see
cc     JCC 15 (1994) 1-11, eqs.(10)-(13)
cccccccccccccccccccccccccccccccccccccccccccccccccccc
cc einmat        Hessian to update      (H_0)       c
cc q (=gamma)   difference of gradient (g_1 - g_0)  c
cc p (=delta)   difference of point    (x_1 - x_0)  c
cc upmat         updated Hessian        (H_1)       c
cc zi            q - (H_0) p                        c
cc Hms           Murtag-Sargent update              ccc
cc Hpo           Powell update                        c
cc phi           optimal combination between Hms/Hpo  c
cc                                                    c
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 einmat,q,p,upmat,zi,unit,b1,v2,
     1     vec_ddot,eps,sw,pzi,phi,pp,zizi,Hms,Hpo
      integer nn
      dimension unit(nn,nn),
     +  einmat(nn,nn),upmat(nn,nn),
     +  b1(nn,nn),sw(nn,nn),v2(nn),
     +  q(nn),p(nn),zi(nn),Hms(nn,nn),Hpo(nn,nn)
      data eps/1.0d-30/
      call vec_dinit(upmat,nn*nn,0.0d0)
      call vec_dinit(Hms,nn*nn,0.0d0)  
      call vec_dinit(Hpo,nn*nn,0.0d0)  
      call mat_diag(unit,nn,1.0d0)     

      call matmult(einmat,nn,p,1,v2,nn)
      call vec_dmult_add(q,v2,nn,-1.0d0,zi)
Cc  Hms  Murtag-Sargent update -- eq.(11)
      pzi=vec_ddot(p,zi,nn)     
ccc      if(dabs(pzi).lt.eps) call cerr(-31)
      call matmult(zi,nn,zi,nn,b1,1)
      call matadd(einmat,1.0d0/pzi,b1,Hms,nn,nn)
Cc  Hpo  Powell update -- eq.(12)
      pp=vec_ddot(p,p,nn)
cccc      if(dabs(pp).lt.eps) call cerr(-31)
      zizi=vec_ddot(zi,zi,nn)
cccc      if(dabs(zizi).lt.eps) call cerr(-31)
      call matmult(p,nn,p,nn,sw,1)
      call matadd(einmat,(-1.0d0)*pzi*pzi/pp/pp,sw,Hpo,nn,nn)
      call matmult(zi,nn,p,nn,b1,1)
      call matadd(Hpo,1.0d0/pp,b1,sw,nn,nn)
      call matmult(p,nn,zi,nn,b1,1)
      call matadd(sw,1.0d0/pp,b1,Hpo,nn,nn)
c c   combination of Hms and Hpo by phi -- eqs.(10,13)
      phi= 1.0d0 - pzi*pzi/pp/zizi 
      call vec_dinit(sw,nn*nn,0.0d0)
      call matadd(sw,(1.0d0-phi),Hms,b1,nn,nn)
      call matadd(b1,phi,Hpo,upmat,nn,nn)
    2 continue
      return  
      end     

      subroutine gettang(curveq,gmat,tangent,nn)
c   This routine calculate the tangent on a curve
c   whereas the tangent is unique determined by:
c
c   H*t=0
c   |t|=1
c   det[H,t]>0
c
      real*8 tangent,gmat,curveq,q,ss1,det,zero,dmnorm
      integer nn,i,j
      dimension tangent(nn),curveq(nn-1,nn),q(nn*nn),ss1(nn,nn)
c      dimension tangent(nn),curveq(nn-1,nn),q(nn,nn),ss1(nn,nn)
      dimension gmat(nn,nn)
      data zero /0.0d0/
      call qr(curveq,q,nn-1)
      call vec_dcopy(q(nn*(nn-1)+1),tangent(1),nn)
c      do 1,i=1,nn
c  1   tangent(i)=q(nn,i)
c wq 
       call vec_dmult_constant
     +   (tangent,nn,1/dmnorm(tangent,gmat,nn),tangent)
       call mat_dcopy(curveq,ss1,nn-1,nn)
       do 10,i=1,nn
   10    ss1(nn,i)=tangent(i)
       if (det(ss1,nn).lt.zero) then
        call vec_dmult_constant(tangent,nn,-1.0d0,tangent)
      endif
      return
      end
c
       subroutine qr(gex,q,n)
C      QR-Zerlegung der Jacobimatrix von ge (gex)
       integer i,j,k,n
       real*8 s1,s2,s,gex,q,r,qq,rr
       dimension gex(n,n+1),q(n+1,n+1),r(n+1,n),
     * qq(n+1,n+1),rr(n+1,n)
       do 10,i=1,n
          do 20,j=1,n+1
             r(j,i)=gex(i,j)
 20       continue
 10    continue
       do 30,i=1,n+1
          do 40,j=1,n+1
             q(i,j)=0.0d0
 40       continue
          q(i,i)=1.0d0
 30    continue
       do 50,i=1,n
          do 60,j=i+1,n+1
             s1=r(i,i)
             s2=r(j,i)
             s=dsqrt(s1*s1+s2*s2)
             if (s.gt.0.0d0) then
                s1=s1/s
                s2=s2/s
                else
                s1=0.0d0
                s2=0.0d0
                endif
             call vec_dcopy(r,rr,(n+1)*n)
             do 70,k=1,n
                rr(i,k)=s1*r(i,k)+s2*r(j,k)
                rr(j,k)=s1*r(j,k)-s2*r(i,k)
 70          continue
             call vec_dcopy(rr,r,(n+1)*n)
             call vec_dcopy(q,qq,(n+1)*(n+1))
             do 80,k=1,n+1
                qq(i,k)=s1*q(i,k)+s2*q(j,k)
                qq(j,k)=s1*q(j,k)-s2*q(i,k)
 80          continue
             call vec_dcopy(qq,q,(n+1)*(n+1))
 60       continue
 50    continue
       do 90,i=1,n+1
          do 100,j=1,n+1
             q(i,j)=qq(j,i)
 100      continue
 90    continue
       return
       end
c
c      subroutine recorrect(corr,cold,gmat,ocorr,nn)
c   manipulation of the corrector for a better convergence
c      integer nn
c      real*8 corr,cold,dmnorm,gmat,vec_mdot
c      logical olen,odir,ocorr,odone
c      dimension corr(nn),cold(nn),gmat(nn,nn)
c      ocorr=.false.
c      if (dmnorm(corr,gmat,nn).gt.dmnorm(cold,gmat,nn)) then
c           olen=.true.
c      else
c           olen=.false.
c      endif
c      if (vec_mdot(corr,gmat,cold,nn).lt.0.0d0) then
c           odir=.true.
c      else
c           odir=.false.
c      endif
c      odone=.false.
c      if (odir.and.olen.and.ocorr) then
c        call vec_dmult_constant(corr,nn,0.6d0,corr)
c        odone=.true.
c      endif
c      ocorr=odone
c      return
c      end
c
      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
c
       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
c
      real*8 function adet(mat,sdet,dim)
c  calculate the Determinant of MAT
      real*8 mat,amat,det,fact,sdet,sz
      integer dim,itest,sp,nexp,k,dreieck,ll
      dimension mat(dim,dim),amat(dim,dim),sp(dim-1),fact(dim)
      sz=1.0d0/16.0d0
      call mat_dcopy(mat,amat,dim,dim)
      itest=dreieck(amat,sp,dim)
      do 5,k=1,dim-1
   5  continue
      if (itest.eq.-1) then
        adet=0.0d0
        sdet=0.0d0
      else
        nexp=0
        det=1
        sdet=1.0d0
        do 10,k=1,dim
                det=det*amat(k,k)
                fact(k)=amat(k,k)
                if (k.lt.dim) then
                        if (k.lt.sp(k)) det=-det
                endif
  20            continue
                if (dabs(det).lt.sz) then
                        det=det*16.0d0
                        nexp=nexp-1
                        goto 20
                endif
  30            continue
                if (dabs(det).ge.1.0d0) then
                        det=det*sz
                        nexp=nexp+1
                        goto 30
                endif
                if (k.eq.dim-1) then
                        sdet=det
                        ll=nexp
                endif
  10    continue
        if (nexp.lt.0) then
                nexp=-nexp
                do 40,k=1,nexp
                        det=det*sz
  40            continue
        elseif (nexp.gt.0) then
                do 50,k=1,nexp
                        det=det*16.0d0
  50            continue
        endif
        if (ll.lt.0) then
                ll=-ll
                do 410,k=1,ll
                        sdet=sdet*sz
  410           continue
        elseif (ll.gt.0) then
                do 510,k=1,ll
                        sdet=sdet*16.0d0
  510           continue
        endif
        adet=det
      endif
      return
      end
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
cccccccccccccccccccccccccccccccccccccccccccccc
c      subroutine toang(vec1,vec2,nn)
cc   (bohr,radian)->(angstroem,degree) for internal coordinates
c      real*8 vec1,vec2,torad,btoa
c      integer i,nr,nn
c      dimension vec1(nn),vec2(nn)
c      data btoa  /0.52917706d0/
c      torad=dacos(-1.0d0)/180.0d0
c      nr=(nn+6)/3-1
c      do i=1,nn
c      if (i.le.nr) then
c         vec2(i)=vec1(i)*btoa
c      else
c         vec2(i)=vec1(i)/torad
c      endif
c      end do
c      return
c      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine angtobohr(point,NN)
        real*8 torad,btoa,point(nn)
        integer i,nn,n3
c (bohr,radian) <<-- (angstroem,degree) for internal coordinates
c        torad=dacos(-1.0d0)/180.0d0
        torad=3.141592653589793d0/180.0d0
        btoa =0.52917706d0
        n3=(nn+6 +2 )/3
c
c        point(1)=point(1)/btoa
c        point(2)=point(2)/btoa
c        point(3)=point(3)*torad
c        If(nn.gt.3) then
c         Do 1 I=1,N3-3
c         point(3*I+1)=point(3*I+1)/btoa
c         point(3*I+2)=point(3*I+2)*torad
c1        point(3*I+3)=point(3*I+3)*torad
c        Endif
c note: fix dih4 , reorganize the list of coordinates
        point(1)=point(1)/btoa
        point(2)=point(2)/btoa
        point(3)=point(3)*torad
        point(4)=point(4)/btoa
        point(5)=point(5)*torad
        If(nn.gt.5) then
         Do 1 I=2,N3-4
         point(3*I)=point(3*I)/btoa
         point(3*I+1)=point(3*I+1)*torad
1        point(3*I+2)=point(3*I+2)*torad
        Endif
         point(nn-1)=point(nn-1)/btoa
         point(nn)=point(nn)*torad
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine bohrtoang(point,NN)
        real*8 torad,btoa,point(nn)
        integer i,nn,n3
c (bohr,radian)-->>(angstroem,degree) for internal coordinates
c        torad=dacos(-1.0d0)/180.0d0
        torad=3.141592653589793d0/180.0d0
        btoa =0.52917706d0
        n3=(nn+6 +2 )/3
c  note dih6 is fixed, also dih22 
        point(1)=point(1)*btoa
        point(2)=point(2)*btoa
        point(3)=point(3)/torad
        point(4)=point(4)*btoa
        point(5)=point(5)/torad
        If(nn.gt.5) then
         Do 1 I=2,N3-4
         point(3*I)=point(3*I)*btoa
         point(3*I+1)=point(3*I+1)/torad
1        point(3*I+2)=point(3*I+2)/torad
        Endif
         point(nn-1)=point(nn-1)*btoa
         point(nn)=point(nn)/torad
        return
        end
ccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_dmult_add(arra1,arra2,ndim,zahl,arra4)
c  For each Komponent: ARRA4 = ARRA1 + ZAHL*ARRA2
      real*8  arra1,arra2,arra4,zahl
      integer ndim,i
      dimension arra1(ndim),arra2(ndim),arra4(ndim)
      if(ndim .le. 0)return
      do 145 i = 1,ndim
       arra4(i) = arra1(i)+arra2(i)*zahl
 145  continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccc
      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
ccccccccccccccccccccccccccccccccccccccccccccccccc
      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 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
      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
      subroutine matadd(dm1,fact,dm2,out,n,m)
c    OUT(N,M) = DM1(N,M) + FACT * DM1(N,M)
      real*8 dm1,fact,dm2,out
      integer m,n
      dimension dm1(n,m),dm2(n,m),out(n,m)
      call vec_dmult_add(dm1,dm2,n*m,fact,out)
      return
      end
cccccccccccccccccccccccccccccccccccccccccccc
      integer function dreieck(mat,sp,dim)
c     Dreiecksfaktorisierung of MAT ; dim(SP) = DIM-1
      real*8 mat,z
      integer sp,dim,it,k,i,s,j
      dimension mat(dim,dim),sp(dim-1)
      it=0
      do 10,k=1,dim-1
        z=0.0d0
        do 20,i=k,dim
                if (dabs(mat(i,k)).gt.z) then
                        z=dabs(mat(i,k))
                        s=i
                endif
  20    continue
        if (z.eq.0.0d0) then
                dreieck=-1
                return
        endif
        sp(k)=s
        if (k.lt.s) then
                do 30,j=1,dim
                        z=mat(k,j)
                        mat(k,j)=mat(s,j)
                        mat(s,j)=z
  30            continue
        endif
        do 40,i=k+1,dim
                mat(i,k)=mat(i,k)/mat(k,k)
  40    continue
        do 50,j=k+1,dim
                do 60,i=k+1,dim
                        mat(i,j)=mat(i,j)-mat(i,k)*mat(k,j)
  60            continue
  50    continue
  10  continue
      if (mat(dim,dim).eq.0.0d0) it=-1
      dreieck=it
      return
      end
cccccccccccccccccccccccc
