PROGRAM lj7sc5
C CCC for LJ_7 potential N=21 CCCCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp  30.06.2004  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c script version - modular structure 
c predictor step
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL*8  Y(52,21),D1,en ,ENall(52) 
      Integer L,N,J,I,LA,ITall,Istatus,itGS
c 7: all points of chain for Mma input
       OPEN(7,FILE='lj7gs.weg') 
c 77: all energy points for Mma input
       OPEN(77,FILE='lj7keten.weg')       
c 9: all points of chain
       OPEN(9,FILE='lj7chain.dat') 
       OPEN(11,FILE='lj7param.dat') 
       OPEN(12,FILE='lj7proj.dat') 
c 13: actual point X of GS iteration
       OPEN(13,FILE='lj7point.dat') 
       OPEN(14,FILE='lj7ener.dat') 
       OPEN(15,FILE='lj7grad.dat') 
       OPEN(44,FILE='protocol.txt',status='unknown',access='append')
C Constants
       N=21 
c       ETA=0.12d0
c       eps=0.075d0
c       L= 25
c  itGS iteration counter in the GS loop 
       itGS=0
       istatus=0
       rewind 11
       READ(11,*)  L, LA, ETA, EPS, ITALL, itGS
       WRITE(44,*) ' predictor step  '
       WRITE(44,*)  L, LA, ETA, EPS, ITALL, itGS
       rewind 9
       Do 11, J=1,L
 11     READ(9,*) (Y(J,I),I=1,N)
c GS gives actual point at node LA-1, on chain to final minimum
c
       rewind 14
       READ(14,*) ENall(LA)
       IF((LA.gt.1) .and. (LA.le.L)) then
         rewind 77
         READ(77,*) (ENall(I),I=1,LA-1)
       ENDIF
       IF(LA.le.L) then
        rewind 77
        WRITE(77,*)(ENall(I),I=1,LA)              
       ENDIF
c
       rewind 13
       READ(13,*)(Y(LA,I),I=1,N)
c
      LA=LA+1
      rewind 11
      WRITE(11,*) L, LA, ETA, EPS ,ITALL , itGS 
      WRITE(44,*)  L, LA, ETA, EPS, ITALL, itGS

      If(LA.eq.L) then
        rewind 13
        write(13,*)(Y(L,I),I=1,N)
        ITall=ITall+itGS
        WRITE(6,25)LA-2,ITall
        itGS=0
        rewind 11
        WRITE(11,*) L, LA, ETA, EPS ,ITALL , itGS
        WRITE(44,*)  L, LA, ETA, 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,N)
 15   continue        
      WRITE(7, 21)(Y(L,I),I=1,N)
c
      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
      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
cc       Do 12, J=LA,L-1
cc       Do 13, I=1,N
cc 13     Y(J,I)=( Y(J-1,I)*(L-LA-1)+ Y(L,I) )/(L-LA)
cc 12    continue 
       rewind 9
       Do 14, J=1,L
 14     WRITE(9,*) (Y(J,I),I=1,N)
c 
       rewind 13
       write(13,*)(Y(LA,I),I=1,N)
       WRITE(44,*) ' guess point of predictor '
       write(44,*)(Y(LA,I),I=1,N) 
c 
       ITall=ITall+itGS
       WRITE(6,25)LA-2,ITall
       WRITE(44,25)LA-2,ITall
       itGS=0
       rewind 11
       WRITE(11,*) L, LA, ETA, EPS ,ITALL , itGS 
       WRITE(44,*) L, LA, ETA, EPS ,ITALL , itGS

CC bis hier Iterations-Zyklus
   
999   CONTINUE              
                
 23   FORMAT(F13.9, 3H ,   ) 
 22   FORMAT(3H { , F13.9, 3H , , F13.9, 3H }, )    
 21   FORMAT(3H { , F13.9, 3H , , F13.9, 3H }}  )
 25   FORMAT(' Node    ',I4,'  ITall    ',I5)
      close(44)
      call exit(istatus)
      STOP                                                       
      END