      PROGRAM A22Vri1
!
!  W.Quapp, Decembre 2009
!  
!    CCC for Alanine potential CCCCCCCCCCCCCCCCCC     
! energy + gradient + Hessian + metric from GamessUS
! Programm bestimmt Suchrichtung der Newtontrajektorie, die den
! VRI-Punkt enthält
!
!      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)
! A22vriReadOUT: 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      2 dof reduced !! (of 60) thus n=58
      PARAMETER (n=58,N3=22,length=12,eps=0.0001d0)
c                         length should be even 
      INTEGER Istatus,iter,Natoms(N3)
      DOUBLE PRECISION chain(length+1,n),vri(n),chainP(length/2+1,n)
      DOUBLE PRECISION xstart(n+1),xend(n+1),dir(n+1),torad,btoa
      Character*2  CharAtoms(N3)
      Character*5  Zcoor(n+2)
      Character*70 zmatLines(35)
c Chass et al. numbering is used
      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     pure and applied chemistry, 51, 1 (1979):
c          angstroms per bohr
      data toang  /0.52917706d0/
c 7: 2D points of initial chain for Mma output
      OPEN(7,FILE='chain.dat')
ccccccccccccccccccccccccccccccccccccccccc
      OPEN(8,FILE='a22Vkette.dat')
      OPEN(9,FILE='a22VketNrj.dat')
      OPEN(11,FILE='a22Vdir.dat')
      OPEN(16,FILE='a22Veps.dat')
      OPEN(17,FILE='a22Viter.dat')
      OPEN(18,FILE='a22Vistat.dat')
      OPEN(19,FILE='a22Vicalc.dat')     
      OPEN(25,FILE='a22Vtold.dat')
      OPEN(26,FILE='a22Vt1old.dat')
      OPEN(30,FILE='a22Vvri.dat')
      OPEN(32,FILE='a22Vnode.dat')
      OPEN(33,FILE='a22Vmin.dat')
c protocol file: all interesting output of the GS
      OPEN(44,FILE='a22Vprotocol.txt',status='unknown',access='append')
cccccccccccccccccccccccccccccccccccccccc
c 80,81: initial minima input:  min1, min2, or SP
c start: C7ax ,  end: C5
       OPEN(88,FILE='obenAla2CHASS')
       OPEN(81,FILE='Ala2Alt08Run60Node19.dat')
       OPEN(82,FILE='Ala2grad.dat')      
       OPEN(83,FILE='Ala2VRI.dat')
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     iter: iteration over |Ag|--> min, 
c     iterNT: count convergence in part 2
      iter=0 
      write(17,*) iter, iter
c 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: import value of last iteration over |Ag|--> min      
      xminimum=0.0004d0
      write(33,*) xminimum, xminimum, 1000.0d0*xminimum
ccccccccccccccccccccccccccccccccccccccccccc
      Write(44,*)
      write(44,*)' Search singular VRI-chain for ala2, 22 atoms'
      write(*,*) ' Search singular VRI-chain for ala2, 22 atoms'
      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,35
88    READ(88,881) zmatLines(I)
      Write(44,*) ' Ansatz in GamessUS:  '
      Do 880 I=1,35
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 
812   Format(F12.7)
c 
      READ(83,*)  (xstart(I),I=1,N)
      READ(81,*)  (xend(I),I=1,N)
c move test VRI for start in direction of SP:
       xstart(20)=  xstart(20) +1.0d0
       xstart(44)=  xstart(44) +0.20d0
       WRITE(44,50)(xstart(I),I=1,5)
       WRITE(44,50)(xstart(I),I=6,N)
       WRITE(*,50) (xstart(I),I=1,5)
       WRITE(*,50) (xstart(I),I=6,N)
       WRITE(44,*) 'ketten-point start phi & psi is '
       WRITE(44,*) '{', xstart(20), ',', xstart(44) ,' },'
       WRITE(7,*)  '{'
       WRITE(7,*)  '{', xstart(20), ',', xstart(44) ,' },'
c  store xstart:
      rewind 32
      WRITE(32,50)(xstart(I),I=1,5)
      WRITE(32,50)(xstart(I),I=6,N)
50    FORMAT(1X,3F20.10)      
c xend is only hypothetical, because it is a predictor-corrector-run
c      WRITE(44,*) ' xend is with 6'
c      WRITE(44,50)(xend(I),I=1,N+1)
c      WRITE(*,*) ' xend is with x6'
c      WRITE(*,50)(xend(I),I=1,N+1)
c read previous VRI run:
      rewind 83
      read(83,*)(vri(I),I=1,N)
c old: set phi and psi value directly, vri(6) is already overwritten:     
c      vri(20)= 115.0d0
c      vri(44)= -72.0d0
      rewind 30
        WRITE(30,50)(vri(I),I=1,N)
        WRITE(30,50)(vri(I),I=1,N)
      WRITE(44,*)' start-vri-approx run 22_7  '
        WRITE(44,50)(vri(I),I=1,5)
        WRITE(44,50) (vri(I),I=6,N)
      WRITE(*,*) 'start-vri-approx run 22_7 '
        WRITE(*,50) (vri(I),I=1,5)
        WRITE(*,50) (vri(I),I=6,N)
ccccccccccccccccccccccccccccccccccccccccccccccccc      
c       DO 5 i=1,n+1
c         dir(i)=xend(i)-xstart(i)
c 5     CONTINUE        
cc   Note: zwischen SPs ist die direkte Richtung oft schlechte Wahl! 
cc   Normierung: 
cc   (bohr,radian) <<-- (angstroem,degree) for internal coordinates
c        torad=dacos(-1.0d0)/180.0d0
c        btoa = 0.52917706d0
c        dir(1)=dir(1)/btoa   
c        dir(2)=dir(2)/btoa 
c        dir(3)=dir(3)*torad
c        dir(4)=dir(4)/btoa
c        dir(5)=dir(5)*torad
c       DO 6 k=2,19         
c        dir(3*k  )= dir(3*k  +1)/btoa
c        dir(3*k+1)= dir(3*k+1+1)*torad
c        dir(3*k+2)= dir(3*k+2+1)*torad
c 6     continue                 
ccccccccccccccccccccccccccccccccccccccccccccccccccccc            
c  58 data coordinates for selected grad - direction
      read(82,*) (dir(i),i=1,n)  !  alt  n+1
c in 11 data from previous run
c     read(11,*) (dir(i),i=1,n)      
      write(*,*) ' start dir   '
        write(*,50) (dir(i),i=1,5)
        write(*,50) (dir(i),i=6,n) 
      write(44,*)'start dir   '
        write(44,50)(dir(i),i=1,5)
        write(44,50)(dir(i),i=6,n)
      write(11,2) (dir(i),i=1,n)
c AnfangsTangenten for direction + - 
       write(25,2) (-1.0d0*dir(i),i=1,n)
       write(26,2) (-1.0d0*dir(i),i=1,n)
c note: the tangent direction can be turnd by multiplication with -1
c       this also applies to (interactive) calculations 
c       in part 2 - NT - along an iteration       
c       use programm helpTurnTang.f 
2     FORMAT(3F20.14)
c c make start chain between xstart to vri, vri to xend
      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