      PROGRAM HCNvri1
!      
!  W.Quapp, Benjamin Schmidt, Decembre 2009      
!      
! Program determines search direction of the singular Newtontrajectory 
! containing the  VRI point
! Example HCN, 3 dimensions 
! energy + gradient + Hessian + metric from GamessUS
! use scratchfile 'scriHCNntVRI' for overall run ! 
!
!      n = Dimension,  N3 = number of atoms
! length = Chainlength (even integer)
!    stl = Steplength
!    eps = Convergence in RGF corrector, 
!          is useless for predictors only
!  chain = Chain between points xstart and xend, 
!          xend is useless for predictors only 
!    dir = Searchdirection RGF method = rsel in RGF corrector
! 
! scriHCNntVRI: is the working code in Linux, combining the      
! followings FORTRAN programms, which use following steps
! 
!     this part is 
! (1) Start: Initial settings, with
!  (i) Fix start- and end point of the chain
!  (ii) Fix searchdirection "dir", and initial guess of "VRI" 
!  (iii) Put straigt chain of "length" between both points
!     later parts:
! (2) To "dir" calculate corresponding NT
! (3) Calculate VRI approximation, 
!     one version: use its gradient for new "dir"
!     other version: use only Ag for a predictor step 
! (4) Organize points for Gamess-US input for (2) or (3)
! HCNvriReadOUT: translates output into usable input files
! (5) Organize Iterations up to convergence of VRI
!     by repeating of (2)-(4) 
!
      INTEGER n,N3,length
      DOUBLE PRECISION eps
c
      PARAMETER (n=3,N3=3,length=50,eps=0.0001d0)
      INTEGER Istatus,iter,Natoms(N3)
      DOUBLE PRECISION chain(length+1,n),vri(n),chainP(length/2+1,n)
      DOUBLE PRECISION xstart(n),xend(n),dir(n),torad,btoa
      Character*2  CharAtoms(N3),Zcoor(n)
      Character*70 zmatLines(20)
      Data Zcoor/'hc','nc','w '/
      Data CharAtoms/'H ','C ','N '/
      Data Natoms/1,6,7/
ccccccccccccccccccccccccccccccccccccccccc
      OPEN(8,FILE='hcnVkette.dat')
ccccccccccccccccccccccccccccccccccccccccc
      OPEN(9,FILE='hcnVketNrj.dat')
      OPEN(11,FILE='hcnVdir.dat')
      OPEN(12,FILE='hcnVproj.dat')
      OPEN(16,FILE='hcnVeps.dat')
      OPEN(17,FILE='hcnViter.dat')
      OPEN(18,FILE='hcnVistat.dat')
      OPEN(19,FILE='hcnVicalc.dat')     
      OPEN(20,FILE='hcnVketst.dat')
      OPEN(21,FILE='hcnVketen.dat')
      OPEN(25,FILE='hcnVtold.dat')
      OPEN(26,FILE='hcnVt1old.dat')
      OPEN(30,FILE='hcnVvri.dat')
      OPEN(32,FILE='hcnVnode.dat')
      OPEN(33,FILE='hcnVmin.dat')
c protocol file: all interesting output of the GS
      OPEN(44,FILE='hcnVprotocol.txt',status='unknown',access='append')
cccccccccccccccccccccccccccccccccccccccc
c 80,81: initial point input:  min1, min2, or SP
c 80 start  
c 81 end - not importend by predictor only
      OPEN(88,FILE='obenVhcn')
      OPEN(80,FILE='HCNbranch1.dat')
      OPEN(81,FILE='HCNbranch2.dat')
c       OPEN(80,FILE='HCNbranch3.dat')
c       OPEN(81,FILE='HCNbranch4.dat')
ccccccccccccccccccccccccccccccccccccccccccc
c     iter: iteration over |Ag|--> min, 
c     iterNT: count convergence in part 2
      iter=0 
      write(17,*) iter, iter
c use istatus=2 for rgf,   =3 for vri,   =10 for end of iteration
      istatus=2 
      write(18,*) istatus
c icalc: alternate use of parts 2/3 or 4 for external gradient, Hessian      
      icalc=0
      sgn=-1.0d0
      write(19,*) icalc,sgn
c xminimum: value of last iteration over |Ag|--> min      
      xminimum=   6.710d-5
      write(33,*) xminimum, xminimum, 1000.0d0*xminimum
ccccccccccccccccccccccccccccccccccccccccccc
      Write(44,*)
      write(44,*)' Search singular VRI-chain for HCN, 3 atoms'
      write(44,*)'          coor(1)=CH,  coor(2)=CN, coor(3)=angle '  
      write(*,*) ' Search singular VRI-chain for HCN, 3 atoms'
      write(*,*) '          coor(1)=CH,  coor(2)=CN, coor(3)=angle '
      Write(44,*)
      write(16,*) eps
      write(*,*) 'eps in rgf-part for NT convergence', eps
      Write(44,*)
      write(44,*)'eps in rgf-part for NT convergence', eps
      Do 88 I=1,14
88    READ(88,881) zmatLines(I)
      Write(44,*) ' Ansatz in the GamessUS pp: '
      Do 880 I=1,14
880   Write(44,881) zmatLines(I)
881   Format(A70)
      Write(44,*)
      write(44,*)' Nodes Length = ',length
      write(*,*)'  Nodes Length = ',length                                                                       
c  Start - chain between minimums, or SPs, xstart and xend 
      READ(80,*)  (xstart(I),I=1,N)
c best run:  Agmin=3.8*10^-7 of VRI_4
       vri(1)=1.37312087d0
       vri(2)=1.46142511d0
       vri(3)=47.6628534d0
         xstart(1)=vri(1)
         xstart(2)=vri(2)
         xstart(3)=vri(3)            
c suche branch 4 , start bei VRI
c         xstart(1)= 1.35871132d0
c         xstart(2)= 1.47504847d0
c         xstart(3)=47.3797572d0
c start on branch 4  
c         xstart(1)= 1.322764d0
c         xstart(2)= 1.4338d0
c         xstart(3)=50.285d0
c near SP
c        xstart(1)=1.155d0
c        xstart(2)=1.169d0
c        xstart(3)=72.0d0
c HNC
c        xstart(1)=2.155d0
c        xstart(2)=1.169d0
c        xstart(3)=2.0d0
c nehme weitest "linken" VRI_1 bzgl. CH
c Ag - min 8.35924075E-05
c      xstart(1)=3.29713077d0
c      xstart(2)=1.18646316d0
c      xstart(3)=70.9500558d0
c next 
c      xstart(1)=3.4075681833d0
c      xstart(2)=1.1871519276d0
c      xstart(3)=76.11793184d0
c next
c Ag min 3.76293277E-05
c       xstart(1)=3.6191325d0
c       xstart(2)=1.18907874d0
c       xstart(3)=70.1254106d0
      write(44,*)' xstart is '
      WRITE(44,50)(xstart(I),I=1,N)
      WRITE(32,50)(xstart(I),I=1,N) 
50    FORMAT(1X,3F20.10)      
      READ(81,*)  (xend(I),I=1,N)
c  branch b4
c      xend(1)=1.403873d0    
c      xend(2)=1.6982133d0
c      xend(3)=43.77998d0
c  branch b2, last point at right border of a run
      xend(1)=1.52932871821135d0
      xend(2)=1.54048346614484d0
      xend(3)=41.02299440723087d0
      WRITE(44,*) ' xend is'
      WRITE(44,50)(xend(I),I=1,N)
c erste Belegung vri point  cccccccccccccccc
c      DO 4 i=1,n
c      vri(i)=0.5d0*(xend(i)+xstart(i))
c4     CONTINUE
c run with  Agmin=3.8*10^-7
       vri(1)=1.37312087d0
       vri(2)=1.46142511d0
       vri(3)=47.6628534d0
c  other run  
c        vri(1)= 1.35871132d0
c        vri(2)= 1.47504847d0
c        vri(3)=47.3797572d0
c C-H disso  Ag minimum  6.7050953E-05
c         vri(1)=3.36082553d0
c         vri(2)=1.18691529d0
c         vri(3)=70.7027663d0
      WRITE(30,50)(vri(I),I=1,N)
      WRITE(30,50)(vri(I),I=1,N)
      WRITE(44,*)' start-vri-approx '
      WRITE(44,50)(vri(I),I=1,N)
      WRITE(*,*) 'start-vri-approx'
      WRITE(*,50) (vri(I),I=1,N)
ccccccccccccccccccccccccccccccccccccccccccccccccc      
c Note: between skew SPs, this direct direction is often bad choice!
      DO 5 i=1,n
        dir(i)=xend(i)-xstart(i)
5     CONTINUE        
c   Normierung, if dir is not from a gradient !!
c   (bohr,radian) <<-- (angstroem,degree) for internal coordinates
        torad=dacos(-1.0d0)/180.0d0
        btoa = 0.52917706d0
        dir(1)=dir(1)/btoa   
        dir(2)=dir(2)/btoa 
        dir(3)=dir(3)*torad
c run  Agmin=3.8*10^-7
        dir(1)=-0.0264043d0
        dir(2)= 0.2925649d0
        dir(3)= 0.0471219d0
c other
c        dir(1)=-0.0309788d0
c        dir(2)= 0.2944857d0
c        dir(3)= 0.0366888d0
c  C-H disso 
c  Ag  -  minimum  =   8.35924075E-05
c vri     3.29713077  1.18646316  70.9500558
c vridir 
c           dir(1)= 0.0382976d0
c           dir(2)=-0.0239331d0
c           dir(3)= 0.0496706d0
      write(*,*) ' start dir   '
      write(*,50) (dir(i),i=1,n) 
      write(44,*)'start dir'
      write(44,50)(dir(i),i=1,n)
      write(11,2) (dir(i),i=1,n)
c initial tangenten direction + - 
      write(25,2) (1.0d0*dir(i),i=1,n)
      write(26,2) (1.0d0*dir(i),i=1,n)
c note: the  direction can be turnd by multiplication with -1
c       this also applies to (interactive) calculations 
c       for file 25: hcnVtold.dat 
c       in part 2 - NT - along a chain       
2     FORMAT(3F20.14)
cccccccccccccccc make start chain between xstart to vri, vri to xend
c                only important if the corrector is used
      CALL straight_chain(chainP,vri,xend,n,length/2)
      Do 20 k=1,length/2+1
      Do 20 i=1,n
20     chain(length/2+k,i)=chainP(k,i)
       if((length/2)*2.ne.length) then
        Do 21 i=1,n
21      chain(length+1,i)=chain(length,i)
       endif   
      CALL straight_chain(chainP,xstart,vri,n,length/2)
      Do 22 k=1,length/2+1
      Do 22 i=1,n
22     chain(k,i)=chainP(k,i)
      rewind 8
      DO 31 k=1,length+1
       WRITE(8,2) (chain(k,i),i=1,n)
31    CONTINUE
      j=1
      rewind 9
      write(9,*) j
      call exit(istatus)
      END
!
! Legt gerade "kette" der Laenge "le" zwischen die Punkte "x" und  "y".
!
      SUBROUTINE straight_chain(kette,x,y,nn,le)
      DOUBLE PRECISION kette(le+1,nn),x(nn),y(nn)
       DO 20 j=1,le+1
        DO 10 i=1,nn
         kette(j,i)=x(i)+((j-1)*(y(i)-x(i)))/le
10      CONTINUE
20     CONTINUE
      END
cccccccccccccccccccccccccccccccccccccccccccccccc