        PROGRAM Ala22ntGS1
C CCC for Alanine potential CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C*********************************************************
C   C   for Alanine-dipeptide
C   energy + gradient + Hessian + metric from GamessUS
c  cccccccccccccccccccccccccccccccccccccccccccccc
C   output in file **.weg, usable for Mathematica input
C   to visualize the calculated path 
C*cccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C W.Quapp October 2007  for Newton Trajectory             c
C GS nach: Baron Peters et al. JCP 120 (2004)7877-7886 c
C and: W.Quapp JCP 122, iss17, (2005) 174106 (11 pages)c
C    "A growing string method for the reaction pathway c 
C     defined by a Newton trajectory"                  c
c script version - modular structure of program        c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C    eps: tolerance
C    L :  chain length, maximum 52 are to defined by the user
C         nodes between minima: L-1, min1 is L=1, min2 is L=L+1 
C    N :  dimension of the problem
C         here N=60   
C    Y(LA,N) is the chain node, LA=1,...,L,L+1  
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  NDIM, L in ntGS? are to put in every program: 1,2,3,5 
C
      Integer NDIM,L,NA,IFLAG
      PARAMETER(NDIM=58,L=30,NA=22)
      REAL*8  Y(L+1,NDIM),Yb(NDIM),Ye(NDIM)
      REAL*8  xm,ym,zm,gesmass
      Integer N,J,I,LA,ITall,Natoms(NA) 
      Character*2  CharAtoms(NA)
      Character*5  Zcoor(NDIM+2)
      Character*70 zmatLines(35)
c
      Data Zcoor/'ch2  ','cc3  ','cch3 ','oc4  ','occ4 ','dih4 ',
     1   'hc5  ','hcc5 ','dih5 ','hc6  ','hcc6 ','dih6 ','nc7  ',
     2   'ncc7 ','dih7 ','cn8  ','cnc8 ','dih8 ','cc9  ','ccn9 ',
     3   'dih9 ','oc10 ','occ10','dih10','hc11 ','hcn11','dih11',
     7   'hn12 ','hnc12','dih12','cc13 ','ccn13','dih13','hc14 ',
     8   'hcc14','dih14','hc15 ','hcc15','dih15','hc16 ','hcc16',
     7   'dih16','nc17 ','ncc17','dih17','cn18 ','cnc18','dih18',
     8   'hc19 ','hcn19','dih19','hn20 ','hnc20','dih20','hc21 ',
     7   'hcn21','dih21','hc22 ','hcn22','dih22'/
      Data CharAtoms/'h','c','c','o','h','h','n','c','c','o','h','h',
     1               'c','h','h','h','n','c','h','h','h','h'/
      Data Natoms/1,6,6,8,1,1,7,6,6,8,1,1,6,1,1,1,7,6,1,1,1,1/
c other z-matrix                          
c      Data Zcoor/'ch2  ','hc3  ','hch3 ','cc4  ','cch4 ','dih4 ',
c     2   'hc5  ','hcc5 ','dih5 ','oc6  ','occ6 ','dih6 ',
c     4   'nc7  ','ncc7 ','dih7 ','cn8  ','cnc8 ','dih8 ',
c     7   'hn9  ','hnc9 ','dih9 ','cc10 ','ccn10','dih10',
c     8   'hc11 ','hcn11','dih11','cc12 ','ccn12','dih12',
c     1   'oc13 ','occ13','dih13','nc14 ','ncc14','dih14',
c     1   'cn15 ','cnc15','dih15','hn16 ','hnc16','dih16',
c     1   'hc17 ','hcn17','dih17','hc18 ','hcn18','dih18',
c     1   'hc19 ','hcn19','dih19','hc20 ','hcc20','dih20',
c     5   'hc21 ','hcc21','dih21','hc22 ','hcc22','dih22'/
c      Data CharAtoms/'h','c','h','c','h','o','n','c','h','c','h','c',
c     1               'o','n','c','h','h','h','h','h','h','h'/
c      Data Natoms/1,6,1,6,1,8,7,6,1,6,1,6,8,7,6,1,1,1,1,1,1,1/
c      Data Zcoor/'ch2  ','hc3  ','cc4  ','hc5  ','oc6  ',
c     2 'nc7  ','cn8  ','hn9  ','cc10 ','hc11 ','cc12 ', 
c     3 'oc13 ','nc14 ','cn15 ','hn16 ','hc17 ','hc18 ',
c     4 'hc19 ','hc20 ','hc21 ','hc22 ',        'hch3 ',
c     5 'cch4 ','hcc5 ','occ6 ','ncc7 ','cnc8 ','hnc9 ',
c     6 'ccn10','hcn11','ccn12','occ13','ncc14','cnc15',
c     7 'hnc16','hcn17','hcn18','hcn19','hcc20','hcc21',
c     8 'hcc22',
c     9 'dih4 ','dih5 ','dih6 ','dih7 ','dih8 ','dih9 ',
c     1 'dih10','dih11','dih12','dih13','dih14','dih15',
c     2 'dih16','dih17','dih18','dih19','dih20','dih21',
c     3 'dih22'/
c                
c     these data were taken from
c     pure and applied chemistry, 51, 1 (1979) unless noted otherwise
c                   angstroms per bohr
      data toang  /0.52917706d0/
c      torad=3.141592653589793d0/180.0d0 
c                   kilograms per atomic mass unit
c      data tokg   /1.6605655d-27/
c                   electrostatic units (esu)
c               per electron charge pure appl. chem., 2, 717 (1973)
c      data toe    /4.803242d-10/
c                   planck constant, joule-seconds
c      data planck /6.626176d-34/
c                   avogadro constant
c      data avog   /6.022045d+23/
c                   joules per calorie
c      data rjpcal /4.184d0/
c                   metres per bohr
c      data tomet  /5.2917706d-11/
c                   joules per hartre
c      data hartre /4.359814d-18/
c                   boltzman constant, joules per kelvin
c      data boltz  /1.380662d-23/
cccccccccccccccccccccccccccccccccccccccc
c 8,81: minima input min1, min2 
c       start is ala2 c7ax 
c       ende  is      c5
c       SP is between 
      OPEN(88,FILE='obenAla2CHASS')
c          OPEN(8,FILE='Ala22SP631G.dat')
      OPEN(8,FILE='startAla2CHASS.dat') 
c          OPEN(81,FILE='endeAla2testGE.dat')  ! for valley to south 
      OPEN(81,FILE='endeAla2CHASS.dat')
cccccccccccccccccccccccccccccccccccccccc                                        
c 9: all points of chain  
       OPEN(9,FILE='ala2chain.dat')
      OPEN(91,FILE='ala2kette.mol')
      OPEN(95,FILE='ala2kette.gau')
ccccccccccccccccccccccccc ccccccccccccccc
      OPEN(11,FILE='ala2param.dat') 
c 13 actual point X of GS iteration
      OPEN(13,FILE='ala2point.dat') 
      OPEN(21,FILE='ala2iflag.dat')
c for gamess / gaussian data: 
      OPEN(33,FILE='ala2poig.dat')
c for molden file:
      OPEN(38,FILE='ala2rpl.dat')
      OPEN(44,FILE='ala2protocol.txt',status='unknown',access='append')
C Constants           
        N=NDIM
        ITall=0
        IFLAG=0
        write(21,*) IFLAG 
        eps=1.000000100010000005000001d-4
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c        WRITE(6,*) 'give eps !! now '                  c
c        Read(5,*) eps                                  c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      Write(44,*)
      write(44,*)' Chain for ala2, 22 atoms ' 
      Write(44,*)
      Do 88 I=1,35
 88    READ(88,881)  zmatLines(I)
      Write(44,*) ' Ansatz in the gamessUS: '
      Do I=1,35
       write(44,881) zmatLines(I)
      enddo
 881  Format(A70)
      Write(44,*)
      PRINT*,' Chain for ala2, 22 atoms ' 
      write(44,*)' Nodes L = ' ,L,'  threshold eps = ',eps 
c  Start - Kette  zwischen den Minima, Yb start
c  note troughout mimic to fix dih4, and dih22 !!!
      READ(8,812) (Yb(I),I=1,5)
      READ(8,812)  Yb(6) 
      READ(8,812) (Yb(I),I=6,N)
811   Format(1X,F13.7)
      write(44,*)' Yb is '
      WRITE(44,50)(Yb(I),I=1,5)
      WRITE(44,50)(Yb(I),I=6,N)
      READ(81,*)  (Ye(I),I=1,5)
      READ(81,*)   Ye(6)
      READ(81,*)  (Ye(I),I=6,N)
812   Format(F12.7)
      WRITE(44,*) ' Ye is'
      WRITE(44,50)(Ye(I),I=1,5)
      WRITE(44,50)(Ye(I),I=6,N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Do 1 J=1,L+1
      Do 1 I=1,N
c   first straigt chain between 2 Minima
      Y(J,I)= ((L+1-J)*Yb(I)+(J-1)*Ye(I))/L
 1    continue
      rewind 91
       Do 92, J=1,L+1
       write(91,*) '     ',N
       write(91,*)
       Do 92, K=1,5
 92    WRITE(91,505) Zcoor(K), Y(J,K) 
       Do  K=6,N
       WRITE(91,505) Zcoor(K+1), Y(J,K)
       enddo

505   FORMAT(1X,A5, F18.8)
      rewind 95
       Do 925, J=1,L+1
       write(95,*) '     ',N
       write(95,*)
       Do 925, K=1,5
 925   WRITE(95,*) Zcoor(K), Y(J,K)
       Do  K=6,N
       WRITE(95,*) Zcoor(K+1), Y(J,K)
       enddo 
      rewind 9
       Do 2, J=1,L+1
       WRITE(9,50) (Y(J,I),I=1,N)
 2     WRITE(9,50)
c actual node is LA, node number is LA-1
c so, min1=Y_1=Yb,   min_fin=Y_(L+1)=Ye  
      LA=1
      rewind 13 
      WRITE(13,50)(Y(LA,I),I=1,N)

      WRITE(44,*) ' first point'
      rewind 33 
      Do 7 J=1,5 
      WRITE(44,*) Zcoor(J), Y(LA,J)
 7    WRITE(33,*) Zcoor(J), Y(LA,J)
      Do  J=6,N
      WRITE(44,*) Zcoor(J+1), Y(LA,J)
      WRITE(33,*) Zcoor(J+1), Y(LA,J)
      enddo            
 
      WRITE(33,*) '  ' 
      WRITE(6,*)  L, LA,  EPS, ITall, N 
      WRITE(44,*) L, LA,  EPS, ITall, N
      rewind 11
      I=1
      J=0
      WRITE(11,*) L, LA,  EPS, I, J
      xm=0.0d0 
      WRITE(38,*) xm
      WRITE(38,*) xm
      WRITE(38,*) xm  
      close(44)
 22   FORMAT(3H { , 21(F13.9, 3H , ), F13.9, 3H }, )    
 21   FORMAT(3H { , 21(F13.9, 3H , ), F13.9, 3H }} )
 50   FORMAT(1X,3F20.10)
      Stop
      End 
