        PROGRAM a162cgGS5 
C CCC for poly ALA15 potential N=486 CCCCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp  30.11.2004 /07.07.2006  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c script version - modular structure
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
c predictor step with z-Mat coordinates 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer N,L,N3
      PARAMETER(N=486,L=12,N3=162,N6=480)
      REAL*8  Y(L,N),Ye(N),Yint(L,N6),ENall(L)
      Integer J,I,LA,ITall,Istatus,itGS,iflag,Natoms(N3)
      LOGICAL FINISH
      Character*2  CharAtoms(N3)
      Character*50 zmatLinesLA1(N3+4),zmatLinesYe(N3+4)
      Character*7  znames(N6) 
      Data CharAtoms/
     1 'C','C','N','C','C','O','C','N','C','C','O','C','N','C','C','O',
     1 'C','N','C','C','O','C','N','C','C','O','C','N','C','C','O','C',
     1 'N','C','C','O','C','N','C','C','O','C','N','C','C','O','C','N',
     1 'C','C','O','C','N','C','C','O','C','N','C','C','O','C','N','C',
     1 'C','O','C','N','C','C','O','C','N','C','C','O','C','N','C','O',
     1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H',
     1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H',
     1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H',
     1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H',
     1 'H','H','H','H','H','H','H','H','H','H','H','H','H','H','H','H',
     1 'H','H'/
c
      Data Natoms/6,6,7,6,6,8,6,7,6,6,8,6,7,6,6,8,
     * 6,7,6,6,8,6,7,6,6,8,6,7,6,6,8,6,
     * 7,6,6,8,6,7,6,6,8,6,7,6,6,8,6,7,
     * 6,6,8,6,7,6,6,8,6,7,6,6,8,6,7,6,
     * 6,8,6,7,6,6,8,6,7,6,6,8,6,7,6,8,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1/ 
cc
c 7: all points of chain for Mma input
      OPEN(7,FILE='a162gs.weg') 
c 77: all energy points for Mma input
      OPEN(77,FILE='a162keten.weg')       
c 9: all points of chain
      OPEN(9,FILE='a162chain.dat') 
      OPEN(11,FILE='a162param.dat') 
      OPEN(12,FILE='a162proj.dat') 
c 13: actual point X of GS iteration
      OPEN(13,FILE='a162point.dat') 
      OPEN(14,FILE='a162ener.dat') 
c       OPEN(15,FILE='a162grad.dat')
      OPEN(21,FILE='a162iflag.dat')
      OPEN(22,FILE='a162cgplFi.dat') 
      OPEN(23,FILE='a162cggsW.dat')
cc      OPEN(32,FILE='a162atoms.dat',status='old')
      OPEN(33,FILE='a162poig.dat')
      OPEN(44,FILE='a162proto.txt',status='unknown',access='append')
cc  z-mat coordinates
        OPEN(83,FILE='starta162.zmat')
        OPEN(84,FILE='endea162.zmat')
        OPEN(85,FILE='a162current.zmat')
        OPEN(89,FILE='a162chainint.zmat')
cccc
c GS has actual point at node LA on chain to final minimum
C Constants
c      N=66  
c      N3=N/3
c      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
      WRITE(44,*) L, LA, EPS, ITALL, itGS
      WRITE(44,*) '  predictor step to LA+1 = ',LA+1 
c
c Notbremse: emergency brake only 
c      IF(itGS.gt. 3*N+3) istatus=10
c
      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 GS searches actual point at node LA+1 on chain to final minimum
ccccccccccccccccccccccccccccccccccccccccccccc
c !                                         c 
      LA=LA+1
c !                                         c
ccccccccccccccccccccccccccccccccccccccccccccc
      rewind 11
      WRITE(11,*) L, LA, EPS ,ITALL , itGS 
      WRITE(44,*)  L, LA, EPS, ITALL, itGS
      WRITE(44,*) ' Predictor: go with LA one step further '
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c new guess point in Cartesian coordinates
C extra besides the script organization:
      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
      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,*) '  '
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
        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)
c   harmonize the sign of the dihedrals:
      Do 311 I=6,N6,3
      IF(Yint(LA-1,I)*Ye(I) .lt. 0.0d0
     1  .and. Dabs(Yint(LA-1,I)).gt. 130.0d0
     2  .and. Dabs(       Ye(I)).gt. 130.0d0  )  then
       write(44,*)'dih - ERROR in Yint(',LA-1,I,') corrected ' 
          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
c
c    new guess points
      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

      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,*) 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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,*)

       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
