PROGRAM Ala22cgGS5
C CCC for ALA 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=12,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
      Character*1 CharAtoms(22)
      Data CharAtoms/'h','c','h','c','h','o','n','c','h','c','h','c',
     1               'o','n','c','h','h','h','h','h','h','h'/
c
c 7: all points of chain for Mma input
      OPEN(7,FILE='ala22gs.weg') 
c 77: all energy points for Mma input
      OPEN(77,FILE='ala22keten.weg')       
c 9: all points of chain
      OPEN(9,FILE='ala22chain.dat') 
      OPEN(11,FILE='ala22param.dat') 
      OPEN(12,FILE='ala22proj.dat') 
c 13: actual point X of GS iteration
      OPEN(13,FILE='ala22point.dat') 
      OPEN(14,FILE='ala22ener.dat') 
      OPEN(15,FILE='ala22grad.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 
c
c Notbremse: emergency brake only 
      IF(itGS.gt. 3*N+3) 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
         READ(32,*) (Natoms(I),I=1,N3)
         rewind 33  
        Do 7, J=1,N3
         write(44,*) Natoms(J), ( Y(LA,3*(J-1)+I ),I=1,3)
 7       WRITE(33,*) Natoms(J), ( Y(LA,3*(J-1)+I ),I=1,3)
         WRITE(33,*) '  '
c
c        write(44,*) ' ar   ' , ( Y(LA,3*(J-1)+I ),I=1,3)
c 7      continue 
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(F15.9, 3H , ), F15.9, 3H }} )
 22   FORMAT(3H { , 21(F15.9, 3H , ), F15.9, 3H }, )    
 23   FORMAT(F15.9, 3H ,   ) 
 25   FORMAT(' Node    ',I4,'  ITall    ',I5)
 50   FORMAT(1X,3F20.10)
      close(44)
      call exit(istatus)
      STOP                                                       
      END