PROGRAM A22cgGS5 
C CCC for ALA 22 potential N=66 CCCCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp  30.11.2004/09.06.2006 for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c script version - modular structure  
c predictor step with z-Mat coordinates  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer N,L,N3
      PARAMETER(N=66,L=32,N3=22,N6=60)
      REAL*8  Y(L,N),Ye(N),Yint(L,N6),ENall(L),eps,ww
      REAL*8  DiStart
      Integer J,I,LA,ITall,Istatus,itGS,iflag,Natoms(N3)
      LOGICAL FINISH, DiLog
      Character*2  CharAtoms(N3)
      Character*50 zmatLinesLA1(N3+4),zmatLinesYe(N3+4)
      Character*7  znames(N6) 
      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'/
      Data Natoms/1,6,1,6,1,8,7,6,1,6,1,6,8,7,6,1,1,1,1,1,1,1/
cc
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') 
cccccccccccccccccccccccccccccccccccc
      OPEN(14,FILE='ala22ener.dat') 
      OPEN(21,FILE='ala22iflag.dat')
      OPEN(22,FILE='ala22cgplFi.dat') 
      OPEN(23,FILE='ala22cggsW.dat')
      OPEN(33,FILE='ala22poig.dat')
      OPEN(35,FILE='ala22Distart.dat')
      OPEN(44,FILE='ala22proto.txt',status='unknown',access='append')
cc  z-mat coordinates
       OPEN(84,FILE='endeWEG.zmat')
       OPEN(85,FILE='ala22current.zmat')
       OPEN(89,FILE='ala22chainint.zmat')
cc/cc
c GS has actual point at node LA on chain to final minimum
C Constants
c                N=66     N3=N/3      N6=N-6
      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, IFLAG, IFLAG 
      istatus=0
c  itGS iteration counter in the GS loop
      rewind 11
      READ(11,*)  J, LA, EPS, ITALL, itGS
      J=L
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     special case of strong uphill pathway V41 H11 zu N14 L=32
      If(LA.gt.24.and.LA.lt.30)then
       EPS=EPS*1.1d0
      endif
      If(LA.gt.29)then
       EPS=0.005d0
      endif
CCCCCCCCCCCCCCCCCCCC
C CCC wq special case of steep slope on PES: v20/v21  ccccccccccc
c       IF((LA.eq.(J+1)/3).or.(LA.eq.(J+2)/3).or.(LA.eq.(J+3)/3))
c     1  Eps=1.11*Eps                        
C CCC use only on Ala2 path from Min to SP, not from Min to Min 
C CCC ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc       
      WRITE(44,*) L, LA, EPS, ITALL, itGS
      WRITE(44,*) ' predictor step to LA+1 = ',LA+1,' being node =',LA  
c
c Notbremse: emergency brake only 
c     IF(itGS .gt. 3*N) istatus=10
      rewind 9
      Do 11, J=1,L
      READ(9,*) J1
      READ(9,*)
      Do 11, JJ=1,N3
 11   READ(9,555) CharAtoms(JJ),(Y(J,3*(JJ-1)+I),I=1,3)
555   Format(1X,A3, 3F16.8)
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) then
         rewind 77
         READ(77,*) (ENall(I),I=1,LA-1)
         rewind 77
         WRITE(77,*)(ENall(I),I=1,LA)              
       ENDIF
       rewind 13
       READ(13,*)(Y(LA,I),I=1,N)
c
c GS5 calculates guess point for node LA+1 on chain to final minimum
ccccccccccccccccccccccccccccccccccccccccccccc
      LA=LA+1
ccccccccccccccccccccccccccccccccccccccccccccc
      rewind 11
      WRITE(11,*) L, LA, EPS ,ITALL, itGS
      WRITE(44,*) ' Predictor: go with LA one step further '  
      WRITE(44,*) L, LA, EPS, ITALL, itGS
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c new guess point in Cartesian coordinates
      Do 712, J=1,L-LA+1
      Do 713, I=1,N   
713   Y(LA-1+J,I)=( Y(LA-1,I)*(L-LA-J+1)+ Y(L,I)*J )/(L-LA+1)
712   continue    
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      rewind 9
      Do 714, J=1,L 
       WRITE(9,*) '   ',N3
       WRITE(9,*) ' '  
      Do 714, JJ=1,N3
714   WRITE(9,555) CharAtoms(JJ),(Y(J,3*(JJ-1)+I),I=1,3)
      rewind 13
      write(13,50)(Y(LA,I),I=1,N)
      rewind 33
      Do 17, J=1,N3
 17    WRITE(33,*) Natoms(J), (Y(LA,3*(J-1)+I),I=1,3)
      WRITE(33,*) '  '
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C     For step normalization in the corrector
      If(LA.le.L) then
       read(35,*) Distart,  DiLog
       Distart=0.0d0
       Do 717, I=1,N
717    Distart=Distart+(Y(LA-1,I)-Y(LA,I))**2
       Distart=Dsqrt(Distart)
       rewind 35
       Write(35,*) Distart, Dilog
       Write(44,*)'Distance of Predictor Point  ', Distart
      Endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      If(LA.eq.L) then
        rewind 13
        write(13,50)(Y(L,I),I=1,N)
        ITall=ITall+itGS
        WRITE(6,25) LA, ITall
        itGS=0
        rewind 11
        WRITE(11,*) L, LA, EPS ,ITALL , itGS
       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)
        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), ' } ' 
        WRITE(6, *) ' --- Sum of Gradient Calculations  -- >> ',ITall
        WRITE(44,*) ' --- Sum of Gradient Calculations  -- >> ',ITall
        goto 999
      ENDIF
c
      rewind 85
      Do 3 I=1,N3+1
3     READ(85,833)  zmatLinesLA1(I)
833   Format(50A)
      Do 333 I=1,N6
333   READ(85,*) znames(I),Yint(LA-1,I)
      Do 4 I=1,N3+1
      READ(84,833) zmatLinesYe(I)
cccc control order of dist, angles, dihedrals ccccccccccccccccc
      If(zmatLinesLA1(I). ne. zmatLinesYe(I)) then
        write(6,*) 'ERROR in z-mat structure is to correct by hand !!'
        write(6,*) '           LINE =  ', I
        write(44,*)'ERROR in z-mat structure is to correct by hand !!'
        write(44,*)'           LINE =  ', I
      endif
4     continue
      Do 41 I=1,N6
41    READ(84,*) znames(I),Ye(I)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccc harmonize the sign of dihedrals cccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
      Do 311 I=6,N6,3
      IF(Yint(LA-1,I)*Ye(I) .lt. 0.0d0
     1  .and. Dabs(Yint(LA-1,I)).gt. 135.0d0 
     2  .and. Dabs(       Ye(I)).gt. 110.0d0 )  then
       write(44,*) 'dih - ERROR in Yint(',LA-1,I,') corrrected' 
          if (Yint(LA-1,I).gt.0.0d0 ) then
              Yint(LA-1,I)=Yint(LA-1,I)-360.0d0  
            else
              Yint(LA-1,I)=Yint(LA-1,I)+360.0d0
          endif
       ENDIF   
311    continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c    new guess point for LA < L in z-matrix format 
      Do 12 J=1,L-La
      Do 111 I=1,N6
111   Yint(LA-1+J,I)= ((L-LA-J+1)*Yint(LA-1,I)+J*Ye(I))/(L-LA+1)
12    continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      rewind 89
      Do 2, J=LA,L-1
      Do 32 I=1,N3+1
32    Write(89,833)  zmatLinesLA1(I)
      Do 321 I=1,N6
321   WRITE(89,*) znames(I), Yint(J,I)
      WRITE(89,*)  
2     continue
      Do 33 I=1,N3+1
33    Write(89,833)  zmatLinesYe(I)
      Do 322 I=1,N6
322   WRITE(89,*) znames(I), Ye(I)
      WRITE(89,*) 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      rewind 85
      Do 3332 I=1,N3+1
3332    Write(85,833)  zmatLinesLA1(I)
      Do 3321 I=1,N6
3321    WRITE(85,*) znames(I), Yint(LA,I)
        WRITE(85,*)
       WRITE(44,*) ' guess point of predictor '
        rewind 33  
        Do 7, J=1,N3
7         write(44,*) Natoms(J), ( Y(LA,3*(J-1)+I ),I=1,3)
       ITall=ITall+itGS
       itGS=0
       WRITE(6,25)LA,ITall
       WRITE(44,25)LA,ITall
       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