       PROGRAM FK3Dmodel
C 2018 for FK model potential, 3 atoms  
c  modified old programs for FK chains !!
c                NDIM = number of atoms,
c                as  = period of side potential
c                a0  = length of spring 
c                k    = spring constant
c ccccccccccccccccccccccccccccccccccccccccccccccccccccc
C   output files *.weg, usable for Mathematica input
C                   to visualize the calculated path 
C*cccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C W.Quapp October in 2007 for Newton Trajectory        c
C GS nach: Baron Peters et al.JCP 120(2004)7877-7886   c
C and: W.Quapp JCP 122,iss17(2005)174106               c
C    "A growing string method for the reaction pathway c 
C     defined by a Newton trajectory"                  c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C    eps:  tolerance
C    L :   chain length 
C          are to defined by the user
C          nodes of the NT: L-2, minimum1 is L=1, minimum2 is L
c          NOTE:  min2= min1 +N/2 *a_s for all atoms        
C    NDIM: dimension of the problem
C    Y(LA,NDIM) are the chain nodes of the NT, LA=1,..,L+1  
c    do not interchange the atoms of the chain !! 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      Integer NDIM,L
      PARAMETER(NDIM=3,L=153)
      REAL*8 Y(L+1,NDIM),Yb(NDIM),Ye(NDIM),X(NDIM),EPS
      REAL*8 rpl(L+2),Rs(NDIM),ProMat(NDIM,NDIM),ENall(L+1)
      real*8 as,a0,sprk,pi,detH,GrNormAll(L+1)
      real*8 energy,hessian(NDIM,NDIM),gradient(NDIM)
      real*8 Told(NDIM)
      Integer N,J,JJ,I,LA,ITall,ist,NFUN,IFLAG,K
       common /paramI/ JJ,LA,ITALL,NFUN,IFLAG,ist
       common /paramR/ EPS
       common /rrpl/ rpl 
       common /Bproj/ ProMat 
       common /ppoint/ X 
       common /yyy/ Y 
       common /eenerg/ energy,ENall,GrNormAll 
       common /ggrad/ gradient 
       common /hhess/ hessian 
       common /rrs/ Rs 
       common /olds/ Told
       common /BBP/ detH
cccccccccccccccccccccccccccccccccccccccc                                        
       OPEN(17,FILE='FKmstart.weg')
c 91: Newton-Raphson solution: often useless by reordering
       OPEN(91,FILE='FKmNRsol.mol')
cccccccccccccccccccccccccccccccccccccccc
c protocol file: all interesting output of the GS calculation 
      OPEN(44,FILE='FKmprotocol.txt',status='unknown',access='append')
c 7: all points of final chain of the NT for Mma input
      OPEN(7,FILE='FKmgs.weg') 
c 77 energies of the final chain points      
      OPEN(77,FILE='FKmketen.weg')
c 79 gradient norms of the final chain points      
      OPEN(79,FILE='FKmketgr.weg')      
      rewind  7
      Write(7,*)'{' 
      rewind 17
      rewind 44
C Constants           
      N=NDIM
      EPS=1.000000040816D-5
      ITall=0
      IFLAG=0
c standard search direction for the FK modoel
       Do 11 I=1,N
11     Rs(I)=1.0d0
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      Write(44,*)
      write(44,*)' linear Chain for FKm, atoms:',N 
      Write(44,*)
      PRINT*,' Chain for FKm, N atoms ' 
      write(44,*)'Nodes L=',L,' threshold eps = ',eps 
c  Start - chain between minimums, Yb start, Ye end
         Yb(1)=1.17262D0
         Yb(2)=6.28319D0
         Yb(3)=11.3937D0      
c min of 3-chain with given data
        pi= 3.14159265359D0
        as= 2.0D0*pi
        a0= 4.0D0/3.d0*pi
        sprk=1.0D0
      WRITE(44,*) 'Parameters of the FK model'
      WRITE(44,*) 'as=',as ,'     = 2 *pi*',as/(2.0D0*pi)
      WRITE(44,*) 'a0=',a0
      WRITE(44,*) 'springK=',sprk
      WRITE(44,*) 
c atoms at Yb, in starrer Kette 
c      Do 91, I=1,N
c91    Yb(I)=a0*(I-1)
c            
      write(44,*)' Yb is '
      WRITE(44,50)(Yb(I),I=1,N)
c  because FK model is periodic per definition, use the periodicity 
      Do 993, I=1,N
993   Ye(I)=Yb(I)+as
      WRITE(44,*) ' Ye is'
      WRITE(44,50) (Ye(I),I=1,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Do 1 I=1,N
      X(I)=Yb(I)
      Do 1 J=1,L+1
c   first straigt chain between 2 Minima
      Y(J,I)= ((L+1-J)*Yb(I)+(J-1)*Ye(I))/L
 1    continue
c actual node is LA, node number is LA-1
c so, min1=Y_1, min_fin=Y_L+1 
c  
       WRITE(17,*) '{'
       Do 6, J1=1,L+1
       WRITE(17,*) (Y(J1,I),I=1,N)
 6     continue       
       WRITE(17,*) '}' 
c       
       LA=1
       goto 7771
c TEST  energy of start chain
       do 7772 I=1,L+1   
       do 7773 J=1,N
7773   X(J)=Y(I,J)        
        call hge
        If(LA .eq.1) then
         WRITE(44,*) ' gradient of start point '
         WRITE(44,50) (gradient(K), K=1,N)
        endif
        WRITE(44,*) ' energy of start chain ', ENall(LA) 
        WRITE(*,*)  ' energy of start chain ', ENall(LA)
7772   continue    
        rewind 77
       Write(77,*)' start chain '
       WRITE(77,*)' { ', ENall(1), ' , '
        Do 78, I=2,LA-1 
78      WRITE(77,*) ENall(I)  , ' , '
        WRITE(77,*) ENall(L)  , ' } '
       Stop
7771   continue                                                     
c
      LA=1
      WRITE(6,*)  L, LA, ITall, N 
      WRITE(44,*) L, LA, ITall, N
      I=1
      J=0
      JJ=L
      ITALL=I 
      NFUN=J
c now it starts really  !!!!!!!!!!!!!!!!!!!!!
c     PES of the problem 
      call hge
c projector
      call ntGS2                        
c predictor     
      call ntGS5
      ist=0
      do 7775 J=1,N
7775   X(J)=Y(LA,J)
ccccccccccccccccccccccccccccccccccc      
555    continue
       call hge
c  corrector
       call ntGS3
       If(ist .eq. 0) goto 555
       call hge
c  predictor   
       call ntGS5
       do 7774 J=1,N
7774   X(J)=Y(LA,J)
c               
       If(ist .lt. 10) goto 555
ccccccccccccccccccccccccccccccccc
c      close(44)
22    FORMAT(3H { , F13.9, 3H , F13.9, 3H }  )
21    FORMAT(3H { , 2(F13.9, 3H , ) )
50    FORMAT(1X,2F20.10)
      Stop
      End 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c
      subroutine ntGS2
c Projector 
C CCC for FKmodel  potential N=3 vorerst  CCCCCCCCCCCCC
C W.Quapp October 2007  for Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
C und   WQ, JCP 122, iss17, (2005) 174106 (11 pages)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer NDIM,L
      PARAMETER(NDIM=3,L=153)
      REAL*8 Rs(NDIM),RsN,EPS,vec_ddot
      REAL*8 PP(NDIM,NDIM),RR(NDIM,NDIM)
      REAL*8 Y(L+1,NDIM),ProMat(NDIM,NDIM)
      Integer N,I,J,JJ,LA,Lha,Lnum,ITALL,NFUN,IFLAG
       common /paramI/ JJ,LA,ITALL,NFUN,IFLAG,ist
       common /paramR/ EPS
       common /Bproj/ ProMat 
       common /yyy/ Y
       common /rrs/ Rs
cc
cc  Rs: search direction 
c   in classical FK model (1,1,...,1,1)
       N=NDIM 
       JJ=L
cccccccccccccccccccccccccccccccccccccc
        Do 7, J=1,N
 7       write(44,*) J, Y(LA,J)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  projector direction: 
        Lnum= 1
c Search direction of NT: generally any direction in N-space
c       READ(83,*)(Rs(I),I=1,N)
c classical FK model: Rs constant
c       Do 313 I=1,N
c313    Rs(I)=1.0d0
c  Normalization
      RsN=dsqrt(vec_ddot(Rs,Rs,N))
c      If(RsN.eq.0.0d0) then
c       write(6,*) ' '
c       write(6,*) 'STOP ,  search direction error '
c       write(44,*) 'STOP ,  search direction error '
c       write(6,*) ' '
c       write(6,*) ' '
c       STOP 
c      endif 
      call vec_dmult_constant(Rs,N,1.0d0/RsN,Rs)
cccccccccccccccccccccccccccccccccccccccccccccc
       call PROJECTOR(N,Rs,ProMat)
cccccccccccccccccccccccccccccccccccccccccccccc
50    FORMAT(1X,3F20.10)
52    FORMAT(1X,2A,1X,3F20.10)
      return
      End
c ccccccccccccccccccccccccccccccccccccccccccccc
       SUBROUTINE PROJECTOR(N,Rs,ProMat) 
c  projector to vector Rs of dim N in the
c  Projection Matrix ProMat of dim (N,N) 
      INTEGER N
      REAL*8 Rs(N),ProMat(N,N),TN,vec_ddot
      INTEGER J,JJ
      TN=dsqrt(vec_ddot(Rs,Rs,N))
c  normalization  
      call vec_dmult_constant(Rs,N,1.0d0/TN,Rs)
C  Projektionsoperator
      Do 34 J =1,N  
      Do 34 JJ=1,N  
      ProMat(J,JJ)=0.0d0
34    continue      
      Do 35 J=1,N   
      ProMat(J,J)=1.0d0
      Do 36 JJ=1,N     
      ProMat(J,JJ)= ProMat(J,JJ) - Rs(J)*Rs(JJ)
36    continue
35    continue
      RETURN  
      END     
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       SUBROUTINE PROJECTION(N,VEC,ProMat)
c  projection of vector VEC of dim N
c  by Projection Matrix ProMat of dim(N,N)
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine matmult(m1,d1,m2,d2,mout,dim)
c  matrix multiplication:
c    m1[d1,dim]*m2[dim,d2]=mout[d1,d2]
      real*8 m1,m2,mout 
      integer d1,d2,dim,i,j,k
      dimension m1(d1,dim),m2(dim,d2),mout(d1,d2)
c
      do 10,i=1,d1
          do 20,j=1,d2
            mout(i,j)=0.0d0
            do 30,k=1,dim  
             mout(i,j)=mout(i,j)+m1(i,k)*m2(k,j)
   30       continue
  20     continue   
 10   continue      
      return        
      end       
cccccccccccccccccccccccccccccccccccccccccccccccccc
c 
      subroutine ntGS3
c corrector       
C 2018 new for FKm potential  
C W.Quapp October 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 GS cycle for RP: NewPoint calculation for ProjectedGradient=0
c corrector step is done by second order methods of NT theory
c program parts by Michael Hirsch (Dissertation 2003, p.81)
Ccccccccccccc
C  the old NT point for predictor step is Y(LA-1,..)
C  the actual point, to start the correction, is Y(LA,..)
C  put first tangent = Y(LA,..) - Y(LA-1,..) 
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc         
      Integer NDIM,L 
      PARAMETER(NDIM=3,L=153)
      REAL*8 energy,gnorm,EPS,vec_ddot,D1,fac
      REAL*8 VV(NDIM),Told(NDIM),Proj(NDIM),X(NDIM)
      REAL*8 Y(L+1,NDIM),ENall(L+1),GrNormAll(L+1)
      REAL*8 Rs(NDIM),ProMat(NDIM,NDIM), hessian(NDIM,NDIM)
      Integer JJ,I,J,Np,Np1,IterMax,IFLAG,NFUN
      Integer LA,N,Nm1,ist,ITALL
       common /paramI/ JJ,LA,ITALL,NFUN,IFLAG,ist
       common /paramR/ EPS
c       common /rrpl/rpl
       common /Bproj/ ProMat 
       common /ppoint/ X
       common /eenerg/ energy,ENall,GrNormAll
       common /ggrad/ VV
       common /hhess/ hessian 
       common /rrs/ Rs 
       common /yyy/ Y
       common /olds/ Told
c input information
        N  =NDIM
        Nm1=NDIM-1
c  max cycles in corrector; if not convergent in 11 cycles, it is useless         
        ITERMAX=25 
        ist=1
       write( 6,*) ' in corr IFLAG is  ',  IFLAG
       write(44,*) ' in corr IFLAG is  ',  IFLAG
c IFLAG=0 indicates an initial entry into corrector
       write(*,*) L, LA, ITALL, NFUN, EPS
       write(44,*)'LA= ',LA,' ITall=',ITALL,'+',NFUN
       IF(LA.eq.L+1) then
         ist=10
         goto 999
       Endif
c  Print parameters
       if( LA.lt.3 .AND. IFLAG.eq.0)then
         write(*,820)
         write(*,840) N, LA
         write(44,820)
         write(44,840) N, LA
         do 451 I=1,N
451      Told(I) = Y(LA,I)- Y(LA-1,I)
         write(44,*) '   Told is '
         write(44,50) (Told(I),I=1,N)
         write(44,*) 'IFLAG is ', IFLAG 
       endif
       if(NFUN .eq. 0) then 
        Do 452, I=1,N
452     X(I)=Y(LA,I)      
       endif
        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
          ist=0
          IFLAG=2
          ITALL = ITALL+NFUN
          NFUN=0 
          write(44,*)'LA=',LA,' ITall=',ITALL,' NFUN=',NFUN
          write(44,*)
          write(44,*) ' STOP corrector  '
          goto 111
        Endif
c  SearchDirection Rs and global ProjectionMatrix 
        IF(LA. le. 2) then
         WRITE(44,*) ' Rs  is  '  
         WRITE(44,50) (Rs(I),I=1,N)
        endif 
         gnorm=vec_ddot(Rs,Told,N)
         WRITE(44,*)' Rs * Told =', gnorm
c       Do 3 J=1,N
c       if(IFLAG.eq.0 .AND. LA.eq.2) WRITE(44,*)' J= ',J, ' ProMat:'
c       if(IFLAG.eq.0 .AND. LA.eq.2) WRITE(44,50) (ProMat(J,I),I=1,N)
c 3      continue
c function of energy, gradient, hessian, use VV --> gradient 
       call hge
       WRITE(44,*)'energy=',energy 
       write(44,*)'gradient '
       write(44,50) (VV(I),I=1,N)
c   norm of gradient
        gnorm= DSQRT(vec_DDOT(VV,VV,N))
        WRITE(6,*) '  gradient norm ',gnorm
        WRITE(44,*) '  gradient norm ',gnorm
cccccccccccccccccccccccccccccccccccccccccc
C   Newton trajectory with Projector to Rs   
C   Do the projektion of gradient    
       call vec_dcopy(VV,Proj,N)
       call PROJECTION(N,Proj,ProMat)
cccccccccccccccccccccccccccccccccccccccccc
       D1= DSQRT(vec_DDOT(Proj,Proj,N))
       WRITE(6,*) ' projected   grad   norm ', D1
       WRITE(44,*)'  projected  grad   norm ', D1
      If(D1 .lt. EPS) THEN
        WRITE(44,*)'at IteGS=',NFUN,' ok proj norm less then', 
     1  ' tolerance EPS ', EPS
        WRITE(44,*)' successful convergence (no ERROR) '
        WRITE(6, *)' IteGS=',NFUN,' = ok - proj norm',D1
        WRITE(7, *) '{' 
        WRITE(7,23) (X(I), I=1,N-1)
        WRITE(7, *)  X(N), '},'  
 23     FORMAT(F15.9,3H , )
        ITALL = ITALL+NFUN
        NFUN=0 
        IFLAG=2 
        ist=1
        GOTO 111
      ENDIF
c     Call the corrector code
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      CALL newtcorr(NDIM,IFLAG,LA,EPS,X,VV,hessian,Rs,Told)     
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       IFLAG= 1 : Corrector step inside the corr - loop
c              2 : return to predictor
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        NFUN=NFUN+1 
        write(44,*) 'NFUN=', NFUN,' at IFLAG',IFLAG
        write(44,*) 'LA  =', LA
        write(6,*)' IFLAG after corr ',  IFLAG
        write(44,*)'IFLAG after corr ',  IFLAG
c         write(44,*)
         tnorm=dsqrt(vec_ddot(Told,Told,N))
         call vec_dmult_constant(Told,N,1.0d0/tnorm,Told)
          I=0
 41       I=I+1
          IF(I.GT.N) THEN
            write(44,*) '      *** corrector finish ***   '
            ITALL=ITALL+NFUN
            NFUN=0
            IFLAG=2 
            ist=1
            GO TO 501
          ENDIF
          IF(DABS(Proj(I)).GT.EPS) THEN
          write(44,*)'Large Comp of ProjGrad(',I,')',Proj(I)
          write(44,*) 
            IFLAG=1
            ist=0
            GO TO 501
          ELSE
            GO TO 41
          ENDIF
 501  continue
      IF(IFLAG.EQ.2) THEN
         write(44,890) energy
         write(44,*) ' point in corrector'
         write(44,50)  (Y(LA,i),i=1,N)
      ENDIF 
      IF(LA.le.L .and. IFLAG.eq.2) write(44,*)
     +  ' New step of Newton Trajectory '
cccccccccccccccccccccccccccccccccccccccccccccc
111   continue
c Istatus 'ist' controls the corrector loop 
      IF(IFLAG.EQ.2) THEN
        ist=1
      ELSE
        ist=0 
      ENDIF
      write(44,*) 'istatus im corrector',ist
999   continue
c      call exit(ist)
       write(44,*) 'ende corrector'
 50    FORMAT(1X,3F20.10)
 51    FORMAT(1X, F25.16)
 52    FORMAT(1X,3F15.10)
 820  format ( ' NT Routine for ProjGrad=0 ')
 840  format ( ' Ndim   =',i6,'   LA = ',i6 )
 890  format (/,' Energy at pes(x*) = ', 1pd16.8)
      return
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
       subroutine newtcorr(nn,IFLAG,LA,EPS,point,grad,hessian,rsel,
     *  told)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC W.Quapp, April 2007 and July 2018 
CC for GS-NT corrector by NT second order method,
CC replace tangent of NT methods by predictor step of GS
CC with parts of the program of Michael Hirsch
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      integer i,j,nn,N,LA
      real*8 EPS,dsqrt,gnorm,dnst,
     * tt,sdet,gesmass,xm,ym,zm,dih2,dih1,vec_mdot,detHess,as,
     * detH,vec_ddot,rgnorm,detk,det,tgt,tg,dmnorm,r12,r23,phi
c
      real*8 rsel(nn),dk(nn,nn+1),hessian(nn,nn),grad(nn),
     * rgeq(nn-1,nn),dstep(nn),projNm1(nn-1,nn),stepmax,
     * sw(nn),rgrad(nn-1),point(nn),unit(nn,nn)
c
      real*8 tang(nn),told(nn),gmat(nn,nn)
      common /BBP/ detH
c
        N = nn
c set all to unit matrix
        call mat_diag(unit,nn,1.0d0)
        call mat_diag(gmat,nn,1.0d0)
        call vec_dinit(dstep,nn,0.0d0)
ccc
        gnorm=dsqrt(vec_ddot(grad,grad,nn)) 
c
C  rsel search direction to (N-1)xN - ProjMat 
       call ortcompl(rsel,projNm1,unit,nn)
        write(44,*) ' hessian matrix '
        Do 7717, J=1,nn
7717    WRITE(44,50) (hessian(J,I),I=1,nn)
c 
        do 7727, I=1,nn
        do 7727, J=1,nn
7727    gmat(I,J)=hessian(I,J)
        detHess=det(gmat,nn)
        write(44,*) '  det Hessian', detHess
        if (detHess*detH .lt. 0.0d0) then 
          write(44,*) ' ? ? BBP point crossed ? ?  '
        endif 
        detH=detHess
        call mat_diag(gmat,nn,1.0d0)
c
       call tri_mat(projNm1,nn-1,nn,unit,nn,hessian,nn,rgeq)
       call vec_dinit(rgrad,nn-1,0.0d0)
       call matmult(projNm1,nn-1,grad,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))
c Tangent 
       call mat_diag(gmat,nn,1.0d0)
       call gettang(rgeq,gmat,tang,nn)
       write(44,*) ' tang is '
       write(44,50) (tang(I),I=1,nn)
cccccccccccccccccccccccccccccccccccccccccccccc 
        call buildk(dk,rgeq,tang,rgrad,nn)
cc
         detk=det(dk,nn)
         write(44,*) '        detK is    ', detk
         tgt=vec_ddot(tang,told,nn)
       write(44,*) ' tang * told is          ', tgt
cccccccccccccccccccccccccccccccccccccccccccccc
c
       if (rgnorm.gt.EPS) then
c  Corrector
        call linsolve(dk,dstep,nn)
        IFLAG=1
        goto 111
       else
c   goto Predictor
        IFLAG=2
        write(44,*)
        write(44,*)'ConvergentStep (no ERROR) go to next Pred, LA=',LA
       endif
        call vec_dcopy(tang,told,nn)
        if(IFLAG.eq.2) return
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c wq normalization of dstep 
c
111     continue
c
        tt=vec_ddot(dstep,dstep,nn)
       write(44,*) ' step-lenth of corrector ', tt
c Schritt 
       write(44,*) ' dstep is '
       stepmax=0.5000002308000047000d0
       if(tt .gt. stepmax) then
       write(44,*)' !! ERROR corr-step shorten to ', stepmax 
        call vec_dmult_add(point,dstep,nn,stepmax/tt,point)
        write(44,50) (stepmax/tt*dstep(I),I=1,nn)
        IFLAG=0
       else
        call vec_dmult_add(point,dstep,nn,1.0d0,point)
        write(44,50) (dstep(I),I=1,nn)
       endif
c
       write(44,*) ' point in corrector '
       write(44,50) (point(I),I=1,nn)
        
c hypothetical Newton-Raphson step to next stationary point
       call mat_diag(gmat,nn,1.0d0)
       call newton(hessian,grad,dstep,dnst,gmat,NN)
       write(44,*) ' dNewtonStep length whould be: ', dnst
       write(91,*) ' dNewtonStep length whould be: ', dnst
       write(91,*) ' det Hessian', detH
       gnorm=dsqrt(vec_ddot(grad,grad,nn))
       WRITE(91,*) '  gradient norm before step ',gnorm
        Do 44, I=1,nn
44      sw(I)=point(I)+dstep(I)
        write(91,*) 'New min by Newton step' 
        write(91,505) (sw(I),i=1,nn)       
505    FORMAT(1X, F18.8) 
        call vec_dcopy(tang,told,nn)
       return 
 50    FORMAT(1X,3F20.10)
       end       
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c                    
      subroutine buildk(dk,req,tang,rgr,nn)
c build the K-Matrix
      real*8 dk,req,tang,rgr
      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 zero durch kleinen Schritt
      dk(nn,nn+1)=0.0d0
      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)
c
       call mat_diag(gmat,nn,1.0d0)
c                       
       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 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 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
        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 .. S.51 Eqn.(21)/(22)
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
      real*8 gesmass,xm,ym,zm
      integer nn,i,j
      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)
       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
      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
      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
      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
C for FKmodel potential NDim=3 CCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp  April 2007 for Newton Trajectory
C GS nach: Baron Peters et al., JCP 120 (2004) 7877
C  and  WQ, JCP 122, iss 17, (2005) 174106 
CCCCCCCCCCCCCCCCCCCCCCCC
c  GS5 makes a predictor step
       subroutine ntGS5 
      Integer N,L
      PARAMETER(N=3,L=153)
      REAL*8 Y(L+1,N),energy,ENall(L+1),ww,grad(N),EPS,X(N),Rs(N)
      REAL*8 rpl(L+2),gmat(N,N),w3(N),ProMat(N,N),hessian(N,N),Told(N)
      REAL*8  GrNormAll(L+1)
      Integer J,JJ,I,LA,ITall,Istatus,itGS,IFLAG,NFUN,ist
       common /paramI/ JJ,LA,ITALL,NFUN,IFLAG,ist
       common /paramR/ EPS
       common /rrpl/ rpl
       common /Bproj/ ProMat 
       common /ppoint/ X
       common /yyy/ Y 
       common /eenerg/ energy,ENall,GrNormAll
       common /ggrad/ grad 
       common /hhess/ hessian 
       common /rrs/ Rs                                                                
       common /olds/ Told
c 7: all points of chain for Mma input
c      OPEN(7,FILE='FKmgs.weg') 
c 77: all energy points for Mma input
c      OPEN(77,FILE='FKmketen.weg')       
c      OPEN(44,FILE='FKmprotocol.txt',status='old',access='append')
c  GS has actual point at node LA on chain to final minimum
      WRITE(44,*) ' Predictor '       
c  Constants
      IFLAG=0
      ist=0
c  itGS iteration counter in the GS loop
      itGS=NFUN
      J=L
      IF(LA.eq.L+1) then
        ist=10
        goto 999
      Endif
cccccccccccccccccccccccccccccccc
c Notbremse: emergency brake only 
c      IF(itGS .gt. 2*N+1) ist=10
ccccccccccccccccccccccccccccccccc
      WRITE(44,*) ' energy before predictor ', ENall(LA)
       WRITE(*,*) ' energy before predictor ', ENall(LA)
        rewind 77
        Write(77,*)' LA is = ',LA
        WRITE(77,*)' { ', ENall(1), ' , '
        Do 78, I=2,LA-1
78      WRITE(77,*) ENall(I)  , ' , '                
        WRITE(77,*) ENall(L)  , ' } ' 
        rewind 79
        Write(79,*)' LA is = ',LA
        WRITE(79,*)' { ', GrNormAll(1), ' , '
        WRITE(44,*) 'NormGrad', GrNormAll(LA) 
        IF(GrNormAll(LA) .lt. 3.3D-03) Write(44,*) '  SP ?? nearby'
        Do 80, I=2,LA-1 
80      WRITE(79,*) GrNormAll(I)  , ' , '
        WRITE(79,*) GrNormAll(L)  , ' } '
c
c GS searches actual point at node LA+1 on chain to final minimum
c make Reaction Path Length up to point LA: rpl
c
       If(LA.eq.1) goto 1
c in FK model: no curvilinear metric is involved ! 
       call mat_diag(gmat,N,1.0d0)
c                        
       do 104 J=1,N
104    w3(J)=(Y(LA-1,J)-Y(LA,J))
       ww=0.0d0
       Do 9, I=1,N
       Do 9, J=1,N
       ww=ww+gmat(I,J)*w3(I)*w3(J)
9      continue
       ww=Dsqrt(ww)
        
       rpl(LA)=rpl(LA-1)+ww
       Write(44,*)' Reaction Path increases here: ',ww
       Write(44,*)' Reaction Path Length (rpl) up to here:',rpl(LA)
       Write(*,*) ' RPL', rpl(LA)
1      continue 
c go to next node of NT
       LA=LA+1
c
       If(LA.eq.L) then
        ITall=ITall+itGS
        WRITE(6,25) LA, ITall
        itGS=0
        Goto 999
       ENDIF
c
       If(LA.gt.L) then
        ist=10   
        WRITE(6, *) ' --- STOP:  Final Step ok   --- '
        WRITE(6, *) ' --- Sum of Gradient Calculations  -- >> ',ITall
        WRITE(44,*) ' --- Sum of Gradient Calculations  -- >> ',ITall
        goto 999
       ENDIF
c new guess point (without a tangent by GS: growing string )
      Do 12, J=1,L-LA
      Do 13, I=1,N   
 13    Y(LA-1+J,I)=( Y(LA-1,I)*(L-LA-J+1)+ Y(L,I)*J )/(L-LA+1)
 12    continue    
c
       WRITE(44,*) ' guess point of predictor '
       Do 7, J=1,N
        X(J)=Y(LA,J)
7      continue
       write(44,50) (X(J),J=1,N)
       ITall=ITall+itGS
       WRITE(*,25) LA,ITall
       WRITE(44,25)LA,ITall
       WRITE(44,26) itGS
       itGS=0
C                                end iteration cycle 
999   CONTINUE              
c      WRITE(44,*)' Exit-Status Predictor ', ist 
 21   FORMAT(3H { , F15.9, 3H ,  ,F15.9, 4H }}  )
 22   FORMAT(3H { , F15.9, 3H ,  ,F15.9, 4H },  )    
 23   FORMAT(F15.9, 3H , ) 
 25   FORMAT(' Node ',I4,' ITall ',I5)
 26   FORMAT(' itGS',I5)
 50   FORMAT(1X,3F20.10)
 51   FORMAT(2(1X,I3,3F20.10))
c      call exit(ist)
      IF(ist.eq.10) STOP
      return                                                       
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c WQ:     23.05.2018
c         Testpotential   -- former GS4 program part 
c  3 dim- PES of FK-model in x_i for 'chain' of 3 atoms
C  for N=3, but it should work for any N 
C  pes[x_i]= V(x,y)= \sum_{i=1}^{N} [1-cos(2 pi x_i / a_s)]+
C                    \sum_{i=1}^{N-1} k/ 2 [x_{i+1}-x_i-a_o]^2
      subroutine hge
c for N=3 atoms
c former anal_xxx (XX,energy,gradient,hessian)
      Integer N,L 
      PARAMETER(N=3,L=153)
      Integer JJ,LA,ITALL,NFUN,IFLAG,I,J,ist
      Real*8 energy,as,a0,sprk,pi,ENall(L+1),GrNormAll(L+1)
      Real*8 hessian(N,N),gradient(N),XX(N),gnorm,vec_DDOT
      common /paramI/ JJ,LA,ITALL,NFUN,IFLAG,ist
      common /ppoint/ XX
      common /eenerg/ energy,ENall,GrNormAll
      common /ggrad/ gradient 
      common /hhess/ hessian
c min of 3-chain with given data
        pi=3.14159265359D0
        pi2=2.0D0*pi
        as = pi2       
        a0 = 4.0D0/3.d0*pi              
        sprk=1.0D0 
c        
      energy=0.0D0
      Do 1, I=1,N-1
1     energy=energy+1.0D0-Dcos(pi2/as*XX(I))+
     *              sprk/2.0D0*(XX(I+1)-XX(I)-a0)**2
      energy=energy+1.0D0-Dcos(pi2/as*XX(N))
      ENall(LA)=energy
c
      gradient(1)=-sprk*(XX(2)-XX(1)-a0)+pi2/as*Dsin(pi2/as*XX(1))
      Do 2, I=2,N-1
2     gradient(I)= sprk*(2.0D0*XX(I)-XX(I-1)-XX(I+1))+
     *             pi2/as*Dsin(pi2/as*XX(I))
      gradient(N)= sprk*(XX(N)-XX(N-1)-a0)+pi2/as*Dsin(pi2/as*XX(N))
      gnorm= DSQRT(vec_DDOT(gradient,gradient,N))      
      GrNormAll(LA)= gnorm
c
      Do 3, I=1,N
      Do 3, J=1,N
      hessian(I,J)=0.0D0
3     continue
      Do 4, I=1,N-1
      hessian(I,I+1)=-sprk
4     hessian(I+1,I)=-sprk
      hessian(1,1)=sprk +( (pi2/as)**2 )*Dcos(pi2/as*XX(1))
      Do 5, J=2,N-1 
      hessian(J,J)=2.0D0*sprk +( (pi2/as)**2 )*Dcos(pi2/as*XX(J))
5     continue     
      hessian(N,N)=sprk +( (pi2/as)**2 )*Dcos(pi2/as*XX(N)) 
      return
      end
cccccc and all !!!!!!!!!!!!!!  June 2018  !!!!!!!!!!!!!!!!!
