PROGRAM A22cgGS5
C CCC for at_22 potential N=66 CCCCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp  30.09.2004  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c script version - modular structure for LJ22 test 
c predictor step
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer N,L,N3
      PARAMETER(N=66,L=25,N3=22)
      REAL*8  Y(L+1,N),ENall(L+1),ww,grad(N),eps 
      Integer J,I,LA,ITall,Istatus,itGS,iflag,Natoms(N3)
      LOGICAL FINISH
c 7: all points of chain for Mma input
      OPEN(7,FILE='at22gs.weg') 
c 77: all energy points for Mma input
      OPEN(77,FILE='at22keten.weg')       
c 9: all points of chain
      OPEN(9,FILE='at22chain.dat') 
      OPEN(11,FILE='at22param.dat') 
      OPEN(12,FILE='at22proj.dat') 
c 13: actual point X of GS iteration
      OPEN(13,FILE='at22point.dat') 
      OPEN(14,FILE='at22ener.dat') 
      OPEN(15,FILE='at22grad.dat')
      OPEN(21,FILE='iflag.dat')
      OPEN(22,FILE='cgplFi.dat') 
      OPEN(23,FILE='cggsW.dat')
      OPEN(32,FILE='alaatoms.dat',status='old')
      OPEN(33,FILE='alapoig.dat')
      OPEN(44,FILE='protocol.txt',status='unknown',access='append')
c GS has actual point at node LA on chain to final minimum
c
C Constants
c      N=66  
c      N3=N/3
      IFLAG=0

      ww=0.0d0
      rewind 21
      write(21,*) IFLAG 
      rewind 23
      write(23,*) (ww,I=1,N+1)
      rewind 22
      FINISH=.FALSE.
      write(22,*) FINISH
      istatus=0

c  itGS iteration counter in the GS loop
      rewind 11
      READ(11,*)  J, LA, EPS, ITALL, itGS
      J=L
      WRITE(44,*) L, LA, EPS, ITALL, itGS
      WRITE(44,*) '  predictor step to LA+1 = ',LA+1 
cc Notbremse: emergency brake only 
      IF(itGS.gt. 2*N+1) istatus=10
c
      rewind 9
      Do 11, J=1,L+1
 11   READ(9,*) (Y(J,I),I=1,N)
c
      rewind 14
      READ(14,*) ENall(LA)
      WRITE(44,*) '  energy before predictor ',  ENall(LA) 
       IF( LA.eq.1)    WRITE(77,*) ENall(LA)
       IF((LA.gt.1) .and. (LA.le.L)) then
         rewind 77
         READ(77,*) (ENall(I),I=1,LA-1)
         rewind 77
         WRITE(77,*)(ENall(I),I=1,LA)              
       ENDIF
c
       rewind 15
ccc       Do 1, J=1,N3
ccc 1     READ(15,*) J1, JJ, (grad( 3*(J-1)+I ),I=1,3)

        READ(15,*) (grad(I),I=1,N) 

c       call PERROR("I/O processing")
c       write(6,*) '               for gradient '  
c       IF(J1 .ne. N3 ) WRITE(44,*) ' ERROR in gradient reading !!! '
c using the special file with grad output from Gaussian03 
       ww=0.0D0 
       Do 3, I=1,N
 3     ww =ww + grad(I)*grad(I)
       ww =Dsqrt(ww)
c       write(6,*) '   gradnorm = ' ,ww, '      at   LA =' , LA
c       write(44,*) '  gradnorm = ' ,ww, '      at   LA =' , LA  
       rewind 13
       READ(13,50)(Y(LA,I),I=1,N)
c
c GS searches actual point at node LA+1 on chain to final minimum
c
      LA=LA+1
      rewind 11
      WRITE(11,*) L, LA, EPS ,ITALL , itGS 
c      WRITE(44,*)  L, LA, EPS, ITALL, itGS

      If(LA.eq.L) then
        rewind 13
        write(13,50)(Y(L,I),I=1,N)
ccc     N3=N/3
        READ(32,*) (Natoms(I),I=1,N3)
        rewind 33  
        Do 17, J=1,N3
 17     WRITE(33,*) Natoms(J), ( Y(L, 3*(J-1)+I ),I=1,3)
        WRITE(33,*) '  '

        ITall=ITall+itGS
        WRITE(6,25) LA, ITall
        itGS=0
        rewind 11
        WRITE(11,*) L, LA, EPS ,ITALL , itGS
c        WRITE(44,*) L, LA, EPS, ITALL, itGS
        Goto 999
      ENDIF
c
      If(LA.gt.L) then
        istatus=10 
        WRITE(6, *) ' ---  STOP: Final Minimum  --- '      
c   Mma - output of Reaction Path  
        rewind 7  
        WRITE(7,*) '{'
        Do 15, J=1,L-1
        WRITE(7,22)(Y(J,I),I=1,N3)
 15     continue        
        WRITE(7,21)(Y(L,I),I=1,N3)
c
        rewind 77
        READ(77,*) (ENall(I),I=1,L)
        rewind 77
        WRITE(77,*) ' { ' 
        WRITE(77,23)(ENall(I),I=1,L-1)    
        WRITE(77,*)  ENall(L), ' } ' 
c
        WRITE(6, *) ' --- Sum of Gradient Calculations  -- >> ',ITall
        WRITE(44,*) ' --- Sum of Gradient Calculations  -- >> ',ITall
        goto 999
      ENDIF
c
c new guess point 
      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
       rewind 9
       Do 14, J=1,L+1
 14     WRITE(9,50) (Y(J,I),I=1,N)
c 
       rewind 13
       write(13,50)(Y(LA,I),I=1,N)
       WRITE(44,*) ' guess point of predictor '
ccc       N3=N/3
c         READ(32,*) (Natoms(I),I=1,N3)
c         rewind 33  
        Do 7, J=1,N3
c         write(44,*) Natoms(J), ( Y(LA,3*(J-1)+I ),I=1,3)
c 7       WRITE(33,*) Natoms(J), ( Y(LA,3*(J-1)+I ),I=1,3)
        write(44,*) ' ar   ' , ( Y(LA,3*(J-1)+I ),I=1,3)
 7      continue 
c       WRITE(33,*) '  '
c 
       ITall=ITall+itGS
       WRITE(6,25)LA,ITall
       WRITE(44,25)LA,ITall
       itGS=0
       rewind 11
       WRITE(11,*) L, LA, EPS ,ITALL , itGS 
       WRITE(44,*) L, LA, EPS ,ITALL , itGS
C                                     end iteration cycle 
999   CONTINUE              
 21   FORMAT(3H { , 21(F13.9, 3H , ), F13.9, 3H }} )
 22   FORMAT(3H { , 21(F13.9, 3H , ), F13.9, 3H }, )    
 23   FORMAT(F13.9, 3H ,   ) 
 25   FORMAT(' Node    ',I4,'  ITall    ',I5)
 50   FORMAT(1X,3F20.10)
      close(44)
      call exit(istatus)
      STOP                                                       
      END