      subroutine baker(istep,nint,dmt,dmt1,dmt2,
     &                 f,vec1,vec2,vec3,vec4,vec5,grad,h)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      dimension simvec(10000)

      dimension dmt(nint,*),dmt1(nint,*),dmt2(nint,*),f(*),vec1(*),
     &          vec2(*),vec3(*),vec4(*),vec5(*),grad(*),h(nint,*)

      data zero,half,four,eigmin/0.0d0,0.5d0,4.0d0,0.02d0/

      save simvec
    
      do 1001 j=1,nint
      do 1001 i=1,nint
 1001 dmt1(i,j)=h(i,j)

      call tred9(nint,nint,dmt1,vec1,vec2,vec3,dmt)
      call tql9(nint,nint,vec1,vec2,dmt,vec3,vec4,ierr)   

      inegative=0

      do 1002 i=1,nint
 1002 if(vec1(i).lt.zero) inegative=inegative+1

      do 1003 i=1,nint
      f(i)=zero
      do 1003 j=1,nint
 1003 f(i)=f(i)+dmt(j,i)*grad(j)

      imode=1

      if(istep.gt.1) then
                           dotmax=zero

                           imode=0   

                           do 1004 i=1,nint

                           dotprod=zero

                           do 1005 j=1,nint
 1005                      dotprod=dotprod+dmt(j,i)*simvec(j)

                           dotprod=abs(dotprod)

                           if(dotprod.gt.dotmax) then    
                                                       imode=i
                                                       dotmax=dotprod
                                                 endif

 1004                      continue
                     endif

      do 1006 i=1,nint
 1006 simvec(i)=dmt(i,imode)

      dunp=half*(vec1(imode)+sqrt((vec1(imode)**2)+four*(f(imode)**2)))

      do 1007 j=1,nint
      do 1007 i=1,nint
 1007 dmt1(i,j)=zero
 
      l=0

      do 1008 i=1,nint
      if(i.ne.imode) then
                           l=l+1
                           dmt1(l,l)=vec1(i)
                           dmt1(l,nint)=f(i)
                           dmt1(nint,l)=f(i)
                     endif
 1008 continue

      call tred9(nint,nint,dmt1,vec5,vec2,vec3,dmt2)
      call tql9(nint,nint,vec5,vec2,dmt2,vec3,vec4,ierr)

      dunn=vec5(1)

      if(imode.eq.1) then
                           dunnmax=vec1(2)
                     else
                           dunnmax=vec1(1)
                     endif

      if(dunn.ge.dunnmax) dunn=dunnmax-eigmin

      do 1009 i=1,nint
 1009 vec2(i)=zero  

      do 1010 i=1,nint
      if(i.eq.imode) eigval=dunp
      if(i.ne.imode) eigval=dunn
      do 1010 j=1,nint
 1010 vec2(j)=vec2(j)-(f(i)*dmt(j,i))/(vec1(i)-eigval)

      return
      end
