PROGRAM A22ntGS5
C CCC for Ala potential NDim=60 CCCCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp June 2007 for Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
C  and   WQ, JCP 122 , iss17, (2005) 174106 (11 pages)
c script version - modular structure 
c predictor step
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer N,L,N3
      PARAMETER(N=58,L=30,N3=22)
      REAL*8  Y(L+1,N),ENall(L+1),ww,grad(N),eps,kcal
     + ,rpl(L+2),gmat(N,N),gmatold(N,N),w3(N),
     +  predstep(N),told(N),pla
      Integer J,I,LA,ITall,Istatus,itGS,iflag,Natoms(N3)
      LOGICAL FINISH, omass
      Character*2 CharAtoms(N3)
      Character*5 Zcoor(N+2)
c CHASS:
      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
c
c      Data Zcoor/'ch2  ','hc3  ','hch3 ','cc4  ','cch4 ','dih4 ',
c     2   'hc5  ','hcc5 ','dih5 ','oc6  ','occ6 ','dih6 ',
c     3   'nc7  ','ncc7 ','dih7 ','cn8  ','cnc8 ','dih8 ',
c     4   'hn9  ','hnc9 ','dih9 ','cc10 ','ccn10','dih10',
c     5   'hc11 ','hcn11','dih11','cc12 ','ccn12','dih12',
c     6   'oc13 ','occ13','dih13','nc14 ','ncc14','dih14',
c     7   'cn15 ','cnc15','dih15','hn16 ','hnc16','dih16',
c     8   'hc17 ','hcn17','dih17','hc18 ','hcn18','dih18',
c     9   'hc19 ','hcn19','dih19','hc20 ','hcc20','dih20',
c     1   'hc21 ','hcc21','dih21','hc22 ','hcc22','dih22'/
cc                                                   
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      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                                         angstroms per bohr
c      data toang  /0.52917706d0/
c
c 7: all points of chain for Mma input
c      OPEN(7,FILE='ala2gs.weg') 
c 77: all energy points for Mma input
      OPEN(77,FILE='ala2keten.weg')       
c 9: all points of chain
      OPEN(9,FILE='ala2chain.dat') 
      OPEN(11,FILE='ala2param.dat') 
c      OPEN(12,FILE='ala2proj.dat') 
c 13: actual point X of GS iteration
      OPEN(13,FILE='ala2point.dat') 
      OPEN(14,FILE='ala2ener.dat') 
      OPEN(15,FILE='ala2grad.dat')
      OPEN(18,FILE='ala2gmatr.dat')
c Pre.-Step:
      OPEN(19,FILE='ala2prest.dat')
c old tangent:      
      OPEN(25,FILE='ala2told.dat')
      OPEN(48,FILE='ala2gmold.dat')
      OPEN(38,FILE='ala2rpl.dat')
      OPEN(21,FILE='ala2iflag.dat')
      OPEN(33,FILE='ala2poig.dat')
      OPEN(44,FILE='ala2protocol.txt',status='unknown',access='append')
      OPEN(54,FILE='ala2dih917p.txt',status='unknown',access='append')
c
c GS has actual point at node LA on chain to final minimum
c Constants
c (bohr,radian)->(angstroem,degree) for internal coordinates
       torad=3.141592653589793d0/180.0d0
c      torad=dacos(-1.0d0)/180.0d0
       btoa =0.52917706d0
      WRITE(44,*)
      WRITE(44,*) ' Predictor '       
      IFLAG=0
      rewind 21
      write(21,*) IFLAG 
      istatus=0
c  itGS iteration counter in the GS loop
      rewind 11
      READ(11,*)  J, LA, EPS, ITALL, itGS
      J=L
      IF(LA.eq.L+1) then
        istatus=10
        goto 999
      Endif
c
ccccccccccccccccccccc
c        IF(LA.le.24) EPS=EPS*1.035
ccccccccccccccccccccc
c     WRITE(44,*) L, LA, EPS, ITALL, itGS
c Notbremse: emergency brake
c      IF(itGS .gt. 2*N+1) istatus=10
c
      rewind 9
      Do 11, J=1,L+1
      READ(9,*) (Y(J,I),I=1,N)
 11   READ(9,*)
c
      rewind 14
      READ(14,*) ENall(LA)
C  include here the ground state energy of Ala_2 !
      kcal=(ENall(LA)+492.637165d0 )*627.52d0
      WRITE(44,*) ' energy before predictor ', ENall(LA),
     +  '  in a.u. absolute ' 
      WRITE(44,*) ' in kcal to global Min: ',kcal  
       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
         Do 771, I=1,LA
 771     WRITE(77,*) ENall(I)  , ' , '                
       ENDIF
c
       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 make Reaction Path Length up to point LA: rpl
c
c       If(LA.eq.1) goto 1
cccccccccccccccccccccccccccccccccccccccccccccccc
cc !! rpl not adapted to fixed dih4, and dih22;
cc    thus does not work here
        Do I=1,LA-1
           rpl(I)=1.0d0
        enddo
ccccccccccccccccccccccccccccccccccccccccccccccccc                         
        If(LA.gt.1) goto 1       
       rewind 18
       Do 10, I=1,N
10      READ(18,50)(gmat(I,J),J=1,N)
        READ(18,*)
       Do 110, I=1,N
110     READ(18,50)(gmat(I,J),J=1,N)
c note gmat is coming after ginv ...       
       rewind 48
       Do 101, I=1,N
101     READ(48,50)(gmatold(I,J),J=1,N)
       rewind 48
       Do 102, I=1,N
102     Write(48,50)(gmat(I,J),J=1,N)
       Do 103, I=1,N
       Do 103, J=1,N
103    gmat(I,J)=(gmat(I,J)+gmatold(I,J))/2.0d0
       rewind 38
       do I=1,LA-1
        READ(38,50) rpl(I)
       enddo
c (bohr,radian)<<--(angstroem,degree) for internal coordinates
c        Do i=1,N3-1
c         w3(i)=(Y(LA-1,i)-Y(LA,i))/btoa
c        enddo
c        Do i=N3,N
c         w3(i)=(Y(LA-1,i)-Y(LA,i))*torad
c        Enddo
        w3(1)=(Y(LA-1,1)-Y(LA,1))/btoa
        w3(2)=(Y(LA-1,2)-Y(LA,2))/btoa
        w3(3)=(Y(LA-1,3)-Y(LA,3))*torad
        If(N.gt.3) then
         Do 123 I=1,N3-3
          w3(3*I+1)=(Y(LA-1,3*I+1)-Y(LA,3*I+1))/btoa
          w3(3*I+2)=(Y(LA-1,3*I+2)-Y(LA,3*I+2))*torad
123       w3(3*I+3)=(Y(LA-1,3*I+3)-Y(LA,3*I+3))*torad
         Endif
        ww=0.0d0
        Do 9, I=1,N
        Do 9, J=1,N
        ww=ww+gmat(I,J)*w3(I)*w3(J)
9      continue
        IF(ww.lt.0.0d0) then
         write(44,*)' ERROR: length step ',ww,'    < 0 ??    '
         ww=-ww
        endif
        ww=Dsqrt(ww)
        rpl(LA)=rpl(LA-1)+ww
       Write(38,50) rpl(LA)
       Write(44,*)'Integrated Reaction Path increases by: ',ww
       Write(44,*)'Reaction Path Length (rpl) up to here: ',rpl(LA)
       Write(*,*) 'integrated rpl ', rpl(LA)
1      continue
      LA=LA+1
      rewind 11
      WRITE(11,*) L, LA, EPS ,ITALL , itGS 
c      WRITE(44,*) 'new LA in predictor '
c      WRITE(44,*)  L, LA, EPS, ITALL, itGS
      If(LA.eq.L) then
        rewind 13
        write(13,50)(Y(L,I),I=1,N)
        rewind 33  
        Do 17, J=1,5
 17      WRITE(33,*) Zcoor(J), Y(L,J)
        Do  J=6,N
         WRITE(33,*) Zcoor(J+1), Y(L,J)
        enddo
        WRITE(33,*) '  '
        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 Extremum  --- '      
c   Mma - output of Reaction Path  
c        rewind 7  
c        WRITE(7,*) '{'
c        Do 15, J=1,L-1
c        WRITE(7,22)(Y(J,I),I=1,N)
c 15     continue        
c        WRITE(7,21)(Y(L,I),I=1,N)
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 of the GS: wq papers of 2005/2007 
      If(LA.lt.35) then 
        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 or do alternately the step along tangent like in old MH programs 
       else 
        WRITE(44,*) ' --- Pred step along TANGENT, no GS pred !!  -- '
        rewind 25
        read(25,*) (told(I),I=1,N)
        call bohrtoang(told,N)
c        told(1)=told(1)*btoa
c        told(2)=told(2)*btoa
c        told(3)=told(3)/torad
c        If(N.gt.3) then
c          Do 120 I=1,N3-3
c           told(3*I+1)=told(3*I+1)*btoa
c           told(3*I+2)=told(3*I+2)/torad
c120        told(3*I+3)=told(3*I+3)/torad
c        Endif
        Do 131, I=1,N
131     Y(LA,I) = rpl(2)*1.50d0*told(I) + Y(LA-1,I)
c and the remainder of the chain in old fashion...
        Do 122, J=2,L-LA
        Do 132, I=1,N
132    Y(LA-1+J,I)=( Y(LA-1,I)*(L-LA-J+1)+ Y(L,I)*J )/(L-LA+1)
122     continue
       Endif
c
       rewind 9
       Do 14, J=1,L+1
        WRITE(9,50) (Y(J,I),I=1,N)
 14     WRITE(9,50)
c 
       rewind 13
       write(13,50) (Y(LA,I),I=1,N)
       WRITE(44,*) ' guess point of predictor '
       rewind 33  
       Do 7, J=1,5
       write(44,51) Zcoor(J), Y(LA,J)
 7     WRITE(33,* ) Zcoor(J), Y(LA,J)
       Do  J=6,N
       write(44,51) Zcoor(J+1), Y(LA,J) 
       WRITE(33,* ) Zcoor(J+1), Y(LA,J)
       enddo
        WRITE(33,* ) '  '
       Do J=1,N
        If(Zcoor(J).eq.'dih9') WRITE(54,*) ' {  ',Y(LA,J-1),' ,'
        If(Zcoor(J).eq.'dih17') WRITE(54,*) Y(LA,J-1),' }, '
       Enddo
       Do J=1,N
        predstep(J)=Y(LA,J)-Y(LA-1,J)
       enddo
       write(19,50) (predstep(J),J=1,N) 
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              
      WRITE(*,*)' Exit-Status Predictor ', istatus
      WRITE(*,*)
      WRITE(44,*)' Exit-Status Predictor ', istatus
      WRITE(44,*) 
 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)
 51   FORMAT(1X,A5,3F20.10)
      close(44) 
      call exit(istatus)
      STOP                                                       
      END
cccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_dinit(arra1,ndim,wert)
c Set all elements of ARRA1 to the double value of WERT
       real*8   arra1,wert
       integer  ndim,i
       dimension arra1(ndim)
       if(ndim .le. 0) return
       do i=1,ndim
         arra1(i)=wert
       end do
       return
       end
ccccccccccccccccccccccccccc
      subroutine vec_dcopy(arra1,arra2,ndim)
c  Copies ARRA1 to ARRA2
      real*8 arra1,arra2
      integer i,ndim
      dimension arra1(ndim),arra2(ndim)
      if(ndim .le. 0) return
      do 10 i = 1,ndim
           arra2(i) = arra1(i)
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccc
      subroutine vec(small,ohoh,u,c,j,k)
c
      implicit real*8  (a-h,p-w),integer   (i-n),logical    (o)
      implicit character *8 (z),character *1 (x)
      implicit character *4 (y)
      dimension c(*),r(3),u(3)
      data dzero/0.0d0/
      jtemp = (j-1)*3
      ktemp = (k-1)*3
      r2 = dzero
      do 20 i = 1 , 3
         r(i) = c(i+jtemp) - c(i+ktemp)
         r2 = r2 + r(i)*r(i)
 20   continue
      r2 = dsqrt(r2)
      ohoh = r2.lt.small
      if (ohoh) return
      do 30 i = 1 , 3
         u(i) = r(i)/r2
 30   continue
      return
      end
cccccccccccccccccccccccccccccc
      subroutine vprod(vp,p,q)
c
      implicit real*8  (a-h,p-w),integer   (i-n),logical    (o)
      implicit character *8 (z),character *1 (x)
      implicit character *4 (y)
      dimension vp(*),p(*),q(*)
c
c     --- evaluate the cross product of x and y
c
      vp(1) = p(2)*q(3) - p(3)*q(2)
      vp(2) = p(3)*q(1) - p(1)*q(3)
      vp(3) = p(1)*q(2) - p(2)*q(1)
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine bohrtoang(point,NN)
        real*8 torad,btoa,point(nn)
        integer i,nn,n3
c (bohr,radian)-->>(angstroem,degree) for internal coordinates
c        torad=dacos(-1.0d0)/180.0d0
        torad=3.141592653589793d0/180.0d0
        btoa =0.52917706d0
        n3=(nn+6 +2 )/3
c  note dih6 is fixed, also dih22
        point(1)=point(1)*btoa
        point(2)=point(2)*btoa
        point(3)=point(3)/torad
        point(4)=point(4)*btoa
        point(5)=point(5)/torad
        If(nn.gt.5) then
         Do 1 I=2,N3-4
         point(3*I)=point(3*I)*btoa
         point(3*I+1)=point(3*I+1)/torad
1        point(3*I+2)=point(3*I+2)/torad
        Endif
         point(nn-1)=point(nn-1)*btoa
         point(nn)=point(nn)/torad
        return
        end

ccccccccccccccccccccccc