      subroutine init(ishake,nmomax,nmos,nocc,naos,istart,ibasis_size,
     &                ibasis,coeff,fock,wkmat,theta,omega,gamma,ccoeff,
     &                overlap)

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

      dimension istart(*),ibasis_size(*),ibasis(naos,*)

      dimension coeff(nmomax,*),fock(*),wkmat(naos,*),theta(naos,*),
     &          omega(naos,*),gamma(nmos,*),ccoeff(naos,*),
     &          overlap(naos,*)

      data zero,hundredth,tenth,one,two,ten/0.0d0,0.01d0,0.1d0,
     &                                      1.0d0,2.0d0,10.0d0/
      data tolerance/0.000000000001d0/

      save cutoff

      if(ishake.le.10) imaxloop=1
      if(ishake.gt.10) imaxloop=20

      if(ishake.eq.1) then
                            cutoff=tenth

                            rewind 45

                            do 1001 i=1,naos
 1001                       call fastrd(45,overlap(1,i),naos)
                      endif

      do 1002 j=1,nmos
      do 1003 i=1,naos
 1003 omega(i,j)=zero
      do 1002 k=1,ibasis_size(j)
      do 1002 i=1,naos
 1002 omega(i,j)=omega(i,j)+overlap(i,ibasis(k,j))*coeff(k,j)

      do 1004 jmo=1,nmos
      do 1005 j=1,ibasis_size(jmo)
 1005 theta(j,jmo)=zero
      do 1004 i=1,naos
      do 1004 j=1,ibasis_size(jmo)
 1004 theta(j,jmo)=theta(j,jmo)+overlap(ibasis(j,jmo),i)*omega(i,jmo)

      do 1006 jmo=1,nocc
      do 1007 i=1,naos
 1007 wkmat(i,jmo)=zero
      do 1006 k=1,ibasis_size(jmo)
      do 1006 i=1,naos
 1006 wkmat(i,jmo)=wkmat(i,jmo)+overlap(i,ibasis(k,jmo))*coeff(k,jmo)

      do 1008 iloop=1,imaxloop

      sum_right_max=zero

      do 1009 imo=1,nocc
      do 1009 jmo=1,imo

      if(jmo.eq.imo) then
                           sum_right=one
                     else
                           sum_right=zero
                     endif

      do 1010 i=1,ibasis_size(imo)
 1010 sum_right=sum_right-coeff(i,imo)*wkmat(ibasis(i,imo),jmo)

      if(abs(sum_right).gt.sum_right_max) then

                                          i_save=imo
                                          j_save=jmo
                                          sum_right_max=abs(sum_right)

                                          endif

      if(abs(sum_right).le.cutoff) goto 1009

      sum_left=zero

      do 1011 i=1,ibasis_size(imo)
 1011 sum_left=sum_left+coeff(i,imo)*theta(i,imo)

      do 1012 j=1,ibasis_size(jmo)
 1012 sum_left=sum_left+coeff(j,jmo)*theta(j,jmo)

      sum_left=sum_left

      if(imo.eq.jmo) sum_left=two*sum_left

      gshake=-(sum_right/sum_left)

      if(abs(gshake).gt.tenth) gshake=tenth*abs(gshake)/gshake

      do 1013 i=1,ibasis_size(imo)
 1013 coeff(i,imo)=coeff(i,imo)-gshake*omega(ibasis(i,imo),jmo)

      do 1014 i=1,ibasis_size(jmo)
 1014 coeff(i,jmo)=coeff(i,jmo)-gshake*omega(ibasis(i,jmo),imo)

      do 1015 i=1,naos
 1015 wkmat(i,imo)=zero

      do 1016 k=1,ibasis_size(imo)
      do 1016 i=1,naos
 1016 wkmat(i,imo)=wkmat(i,imo)+overlap(i,ibasis(k,imo))*coeff(k,imo)

      do 1017 i=1,naos
 1017 wkmat(i,jmo)=zero

      do 1018 k=1,ibasis_size(jmo)
      do 1018 i=1,naos
 1018 wkmat(i,jmo)=wkmat(i,jmo)+overlap(i,ibasis(k,jmo))*coeff(k,jmo)

      sum_right=zero

      do 1019 i=1,ibasis_size(imo)
 1019 sum_right=sum_right+coeff(i,imo)*wkmat(ibasis(i,imo),imo)

      factor=sqrt(one/sum_right)

      do 1020 i=1,ibasis_size(imo)
 1020 coeff(i,imo)=coeff(i,imo)*factor

      do 1021 i=1,naos
 1021 wkmat(i,imo)=wkmat(i,imo)*factor

      sum_right=zero

      do 1022 i=1,ibasis_size(jmo)
 1022 sum_right=sum_right+coeff(i,jmo)*wkmat(ibasis(i,jmo),jmo)

      factor=sqrt(one/sum_right)

      do 1023 i=1,ibasis_size(jmo)
 1023 coeff(i,jmo)=coeff(i,jmo)*factor

      do 1024 i=1,naos
 1024 wkmat(i,jmo)=wkmat(i,jmo)*factor

 1009 continue

      write(6,2223) ishake,iloop,sum_right_max,i_save,j_save
 2223 format(' checking ',2i7,f20.12,2i5)

      if(sum_right_max.lt.tolerance) then

            if(ishake.ne.50) then
                                   return
                             else
                                   do 1025 j=1,nocc
                                   do 1025 i=1,ibasis_size(j)
 1025                              wkmat(i,j)=coeff(i,j)

                                   do 1026 j=1,nocc
                                   do 1026 i=1,naos
 1026                              coeff(i,j)=zero

                                   do 1027 j=1,nocc
                                   do 1027 i=1,ibasis_size(j)
 1027                              coeff(ibasis(i,j),j)=wkmat(i,j)

                                   return
                             endif

                                     endif

 1008 cutoff=sum_right_max/ten

      if(ishake.eq.50) stop ' failed to create initial set of MOs'

      return
      end
