********************************************************************************
*                                                                              *
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### inifun ###
c
c     Initializes orbitals and potentials.
c
      subroutine inifun (cw_orb,cw_coul,cw_exch,cw_suppl,cw_sctch)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension cw_orb(*),cw_coul(*),cw_exch(*),
     & 		cw_suppl(*),cw_sctch(*)

c     zeroizing arrays explicitly

      do i=1,length1
         cw_orb(i)=0.d0
      end do

      do i=1,length2
         cw_coul(i)=0.d0
      end do

      do i=1,length3
         cw_exch(i)=0.d0
      end do

c      return

      if (ini.eq.1) then

c        to initialize properly functions and orbital energies in the case of
c        the HFS calculations islat (which is nonzero) parameter has to be kept 
c        zero before the program reaches the scf routine.

         if (islat.ne.0) islat=0
         call inihyd (cw_orb,cw_coul,cw_exch,
     &        cw_suppl(i4b(7)),cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &        cw_sctch(i5b(1)))
         idump=1
      elseif (ini.eq.2) then

c  initial values of orbitals are provided by the GAUSSIAN92 program

         write(*,*) 'Initializing orbitals using GAUSSIAN94 output'
         call prepg94
         call inigauss(cw_orb,cw_coul,cw_exch,
     &        cw_suppl(i4b(7)),cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &        cw_sctch(i5b(1)))
         idump=1
      elseif (ini.eq.3) then

         write(*,*) 'Initializing orbitals using GAUSSIAN94 long output'
         call prepg94l
         call inigauss(cw_orb,cw_coul,cw_exch,
     &        cw_suppl(i4b(7)),cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &        cw_sctch(i5b(1)))
         idump=1
         
      elseif (ini.eq.4) then
         idump=0
         call rfdisk (cw_orb,cw_coul,cw_exch,cw_sctch(i5b(1)))

c        To generate a higher lying one-electron diatomic state all the lower
c        ones have to be retrieved from disk file and kept for the sake of 
c        the orthogonalization.
c        The current highest orbital is being initialized as a linear combination 
c        of hyfrogenic atomic orbitals.
c        When the incomplete orbital disk file is encountered the values of ini 
c        is changed from 5 to 4.
         
         call inihyd (cw_orb,cw_coul,cw_exch,
     &        cw_suppl(i4b(7)),cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &        cw_sctch(i5b(1)))
         idump=1

      elseif (ini.eq.5) then
c        see ini=4
         idump=0
         call rfdisk (cw_orb,cw_coul,cw_exch,cw_sctch(i5b(1)))

         if (ini.eq.4) then
            call inihyd (cw_orb,cw_coul,cw_exch,
     &           cw_suppl(i4b(7)),cw_suppl(i4b(9)),cw_suppl(i4b(14)),
     &           cw_sctch(i5b(1)))
            idump=1
         endif
      elseif (ini.eq.6) then
         idump=0
         call rfdisk (cw_orb,cw_coul,cw_exch,cw_sctch(i5b(1)))
         print *,'... initializing exchange potentials ...'
	 call tfpot(cw_orb,cw_coul,cw_exch,
     &        cw_suppl(i4b(7)),cw_suppl(i4b(9)),
     &        cw_sctch(i5b(1)))
      endif
      return
      end