********************************************************************************
*                                                                              *
*  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 ### inihyd ###  
c
c     This routine initializes molecular orbitals as linear
c     combinations of hydrogenic functions on centres A and B.
c     In the case of HF or HFS calculations Coulomb (exchange) potentials are 
c     approximated as a linear combination of Thomas-Fermi (1/r) potentials 
c     of the two centres; if method OED is chosen the potential functions are 
c     set to zero.
c
      subroutine inihyd (psi,pot,excp,f2,f4,wgt2,wk0)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension psi(*),pot(*),excp(*),f2(*),f4(*),wgt2(*),wk0(*)

c     Initialization of molecular orbitals 

        print *,'... initializing molecular orbitals from hydrogenic',
     &          ' functions ...'

c       loop over orbitals

        if (ini.eq.4.) then
           norbt=1
        else
           norbt=norb
        endif

	do iorb=1,norbt
           
           ishift=i1b(iorb)-1
           n1=mgx(1,iorb)
           l1=mgx(2,iorb)
           m1=mgx(3,iorb)
           ez1=eza1(iorb)

c          normalization factor for Laguere polynomials

           fn1=(2.d0*ez1/dble(n1))**(3.d0/2.d0+dble(l1))
     &          *sqrt(factor(n1+l1)/(2.d0*dble(n1)*factor(n1-l1-1)))
     &          /factor(2*l1+1)

c          normalization factor for spherical harmonics

           shn1=(-1.d0)**m1/sqrt(4.d0*pii)*
     &		sqrt((2*l1+1)*factor(l1-m1)/factor(l1+m1))
           if (m1.eq.0) shn1=1.d0/sqrt(4.d0*pii)*
     &		sqrt((2*l1+1)*factor(l1-m1)/factor(l1+m1))

c          multiplication factor fot the associate Legendre polynomials
c          alp1=(-1.d0/2.d0)**m1*factor(l1+m1)/factor(l1-m1)/factor(m1)

           n2=mgx(4,iorb)
           l2=mgx(5,iorb)
           m2=mgx(6,iorb)
           ez2=eza2(iorb)	  	
           fn2=(2.d0*ez2/dble(n2))**(3.d0/2.d0+dble(l2))
     &          *sqrt(factor(n2+l2)/(2.d0*dble(n2)*factor(n2-l2-1)))
     &          /factor(2*l1+1)

           shn2=(-1.d0)**m2/sqrt(4.d0*pii)*
     &		sqrt((2*l2+1)*factor(l2-m2)/factor(l2+m2))

c	  alp2=(-1.d0/2.d0)**m2*factor(l2+m2)/factor(l2-m2)/factor(m2)     
           
c          loop over grid points

           do imu=1,mxnmu
              inioff=(imu-1)*nni
              do in=1,nni
                 igp=ishift+inioff+in

c                for each grid point, e.i. for (vmu(imu),vni(ini)) 
c                determine its distance |_r1| and |_r2| from the nuclei A 
c                and B and cosine of the polar angles costh1 and costh2 
c                between z axis and the vectors _r1 and _r2
c
c                rr=dsqrt(xi(j,6)-eta(i,3)+1.d0)
c                rr=dsqrt(vxi12(j)-veta12(i)+1.d0)	
c                rr=dsqrt(vxisq(j)+vetasq(i)-1.d0)	
c                costh=eta(i,1)*xi(j,4)/rr

                 rr=(r/2.d0)*sqrt(vxisq(imu)+vetasq(in)-1.d0)		
                 z=(r/2.d0)*vxi(imu)*veta(in)
                 r1t=(r/2.d0)*(vxi(imu)+veta(in))
                 r2t=(r/2.d0)*(vxi(imu)-veta(in))
                 if (r1t.lt.precis) then
                    costh1=0.d0
                    r1t=1.0d-10
                    psi1=0.d0
                    psi2=0.d0
                    goto 99
                 else   
                    costh1=(z+r/2.d0)/r1t
                 endif
c                
                 if (r2t.lt.precis) then
                    costh2=0.d0
                    r2t=1.0d-10
                    psi1=0.d0
                    psi2=0.d0
                    goto 99
                 else   
                    costh2=-(z-r/2.d0)/r2t	
                 endif
c                
c                calculate radial part of the hydrogenic orbital centered 
c                on both the nuclei
                 
                 psi1=0.d0
                 psi2=0.d0

                 if     (co1(iorb).ne.0.d0) then
                    psi1=fn1*plaguer(n1,l1,ez1,r1t)
     &                   *shn1*plegendg(l1,m1,costh1)
                 endif

                 if (co2(iorb).ne.0.d0) then	    
                    psi2=fn2*plaguer(n2,l2,ez2,r2t)
     &                   *shn2*plegendg(l2,m2,costh2)	  
                 endif
00099            continue
                 psi(igp)=co1(iorb)*psi1+co2(iorb)*psi2
                 if(m1.eq.1) psi(igp)=-psi(igp)
              enddo
           enddo

c          call norm94 (iorb,psi,f4,wgt2,wk0,xnorm)
c          write (*,1115) iorn(iorb),bond(iorb),xnorm,ipnzero

01114   format(/1x,'  orbital        norm       (# primitive bf)')
01115   format(1x,i2,1x,a5,d20.12,' (',i4,')')
        enddo

c  initialize Coulomb and exchange potentials

        call tfpot(psi,pot,excp,f2,f4,wk0)	

        return
	end