********************************************************************************
*                                                                              *
*  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 ### gauss_ovlap ###
c
c     Evaluates the overlap matrix of gaussian basis functions.
c
      subroutine gauss_ovlap (psi,pot,excp,f2,f4,wgt2,wk0)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
c set the maximum value of an argument of the exp function 
      parameter (exp_max=700.d0)
      dimension psi(*),pot(*),excp(*),f2(*),f4(*),wgt2(*),wk0(*)

      do ipb1=1,npbasis
         do ipb2=ipb1,npbasis
            ovl(ipb1,ipb2)=0.d0
            ovl(ipb2,ipb1)=0.d0
         enddo
      enddo

      igauss=0
      if (igauss.eq.0) then
         ishift=0  
         ishift2=nni*mxnmu
         
c        outer loop over gaussian basis set functions

         do i=1,ishift2
            psi(i)=0.d0
         enddo
         do ipb1=1,npbasis
            l1  =lprim(ipb1)
            mp1  =mprim(ipb1)
            if (mp1.lt.0) goto 999
            m1  =iabs(mprim(ipb1))
            d1  =primexp(ipb1)
            if (d1.gt.99999999.d0) goto 999
            ic1 =icgau(ipb1)
            fn1 =fngau2(ipb1)
            shn1=shngau(ipb1)
            fnorm1=fn1*shn1
            
01111       format(5i4,f20.9,3d15.5)
01112       format(4i4,f20.9)
01113       format(5i4,2f20.9)
01114       format(4i4,2f20.9)

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, i.e. for (vmu(imu),vni(ini)) determine 
c                 its distance |_r1| from the nuclei Z_1 and Z_2 and cosine 
c                 of the polar angles costh1 and costh2 between z axis and 
c                 the vectors _r1 and _r2

c                 rr=(r/2.d0)*sqrt(vxisq(imu)+vetasq(in)-1.d0)		

                  z=(r/2.d0)*vxi(imu)*veta(in)
                  if (ic1.eq.1) then
c                    gaussians centred on Z1
                     r1=(r/2.d0)*(vxi(imu)+veta(in))
                     if (r1.lt.precis) then
                        costh1=0.d0
                        psi1=0.d0
                        goto 99
                     else   
                        costh1=(z+r/2.d0)/r1
                     endif
                  elseif (ic1.eq.3) then
c                    gaussians centred on Z2
                     r1=(r/2.d0)*(vxi(imu)-veta(in))
                     if (r1.lt.precis) then
                        costh1=0.d0
                        psi1=0.d0
                        goto 99
                     else   
                        costh1=(z-r/2.d0)/r1	
                     endif
                  elseif (ic1.eq.2) then
c                    gaussians centred on the bond centre
                     r1=sqrt(vxisq(imu)+vetasq(in)-1.d0)*r/2.d0
                     if (r1.lt.precis) then
                        costh1=0.d0
                        psi1=0.d0
                        goto 99
                     else   
                        costh1= (z)/r1	
                     endif
                  endif

c                 calculate the gaussian function centered on one of the nuclei
                  
                  d1r=d1*r1*r1
                  if (d1r.gt.exp_max) then
                     expt=0.d0
                  else
                     expt=1.d0/exp(d1r)
                  endif
                  psi1=fnorm1*r1**l1*expt*plegendg(l1,m1,costh1)
00099             continue
                  psi(igp)=psi1
               enddo
            enddo

c           c1.ne.0 inner loop over gaussian functions
            do i=1,ishift2
               psi(ishift2+i)=0.d0
            enddo
            do ipb2=ipb1,npbasis
               l2  =lprim(ipb2)
               mp2  =mprim(ipb2)
               if (mp2.ne.mp1) goto 1999
               m2  =iabs(mprim(ipb2))
               d2  =primexp(ipb2)
               if (d2.gt.99999999.d0) goto 1999
               ic2 =icgau(ipb2)
               fn2 =fngau2(ipb2)
               shn2=shngau(ipb2)
               fnorm2=fn2*shn2
               do imu=1,mxnmu
                  inioff=(imu-1)*nni
                  do in=1,nni
                     igp=ishift2+inioff+in
                     z=(r/2.d0)*vxi(imu)*veta(in)
                     if (ic2.eq.1) then
c                       gaussians centred on Z1
                        r1=(r/2.d0)*(vxi(imu)+veta(in))
                        if (r1.lt.precis) then
                           costh1=0.d0
                           psi1=0.d0
                           goto199
                        else   
                           costh1=(z+r/2.d0)/r1
                        endif
                     elseif (ic2.eq.3) then
c                       gaussians centred on Z2
                        r1=(r/2.d0)*(vxi(imu)-veta(in))
                        if (r1.lt.precis) then
                           costh1=0.d0
                           psi1=0.d0
                           goto199
                        else   
                           costh1=(z-r/2.d0)/r1	
                        endif
                     elseif (ic2.eq.2) then
c                       gaussians centred on the bond centre
                        r1=sqrt(vxisq(imu)+vetasq(in)-1.d0)*r/2.d0
                        if (r1.lt.precis) then
                           costh1=0.d0
                           psi1=0.d0
                           goto199
                        else   
                           costh1= (z)/r1	
                        endif
                     endif

c                    calculate the gaussian function centered on one of
c                    the nuclei

                     d2r=d2*r1*r1
                     if (d2r.gt.exp_max) then
                        expt=0.d0
                     else
                        expt=1.d0/exp(d2r)
                     endif
                     psi1=fnorm2*r1**l2*expt*plegendg(l2,m2,costh1)
00199                continue
                     psi(igp)=psi1
                  enddo
               enddo

c              evaluate overlap between ipb1 and ipb2 functions
               
               call prod2 (ishift2,psi(1),psi(ishift2+1),wk0)
               call prod  (ishift2,f4,wk0)

               xovlap=ddot (ishift2,wgt2,1,wk0,1)

               if (ipb1.eq.ipb2) then
                  write(*,1114) ipb1,ic2,l2,mprim(ipb2),d2,xovlap
                  call flush(iout6)
               endif

               ovl(ipb1,ipb2)=xovlap
               ovl(ipb2,ipb1)=xovlap

c              end inner loop over basis functions 
               
01999          continue
            enddo
00999       continue
c           end outer loop over basis functions
         enddo

c        save the matrix on disk

         open(30,file='overlap.data',status='unknown',
     &        form='unformatted')
         write(30) ovl
         close(30)
      endif

c     check normalization of finite basis molecular orbitals
      
      open(30,file='overlap.data',status='unknown',form='unformatted')
      read(30) ovl
      close(30)
      
      do iorb=1,norb
         xmon=0.d0
         do ipb1=1,npbasis
            do ipb2=1,npbasis
               xmon=xmon+primcoef(iorb,ipb1)*primcoef(iorb,ipb2)*
     &              ovl(ipb1,ipb2)
            enddo
         enddo
02123    format(2i5,2i3,4d15.4)
         write(*,*) 'norm of mol.orb ',iorb,xmon
      enddo

      return
      end