********************************************************************************
*                                                                              *
*  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 ### locenergy ###
c  
c     Calculates the locale energy for a given orbital.
c
      subroutine locenergy (iorb,psi,pot,excp,e,f0,f4,wgt1,wgt2,
     &			wk0,wk1,wk2,wk3)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      character*13 fn
      include 'commons8.inc'
      dimension psi(*),pot(*),excp(*),e(*),f0(*),wgt1(*),wgt2(*),
     &          wk0(*),wk1(*),wk2(*),wk3(*),f4(*)

      ioutmat=30

      open(ioutmat,file='vni',status='unknown',form='formatted')
      write(ioutmat,*) nni
      write(ioutmat,1000) (vni(in),in=1,nni)
      close(ioutmat)

      open(ioutmat,file='vmu',status='unknown',form='formatted')
      write(ioutmat,*) i1mu(1)
      write(ioutmat,1000) (vmu(im),im=1,i1mu(1))
      close(ioutmat)
01000 format(F25.15)

      wtwoel=0.d0

      i1beg=i1b(iorb)
      i2beg=i2b(iorb)	
      nmut=i1mu(iorb)
      ngorb=i1si(iorb)
      ngpot=i2si(iorb)	

      isd=iorb+(iorb-1)*norb
      engo(isd)=0.0d0
      ipe=mgx(6,iorb)
      ipex=ipe-2*(ipe/2)

c     calculate derivatives over mu and ni 

      if (ipex.eq.0) then
         isym= 1
      else
         isym=-1
      endif
      
      call fill (nni,nmut,isym,psi(i1beg),wk3)
      call difni (nmut,wk3,wk0,wk1,wk2)
      call refill(nni,nmut,wk1,wk0)
      
      call difmu (nmut,wk3,wk2)
      call refill(nni,nmut,wk0,wk2)

c     add contribution from derivatives over mu and ni
      
      call add (ngorb,wk0,wk1)
      
      if (ipe.eq.0) then
         call dcopy (ngorb,f0,1,wk0,1)
      else
         
c        e enter the expression with minus sign which is already 
c        incorporated in e
         
         w=dble(ipe*ipe)
         call dcopy (ngorb,f0,1,wk0,1)
         call daxpy (ngorb,w,e,1,wk0,1)
      endif
      
      call proda (ngorb,psi(i1beg),wk0,wk1)
      call prod (ngorb,wgt1,wk1)

c  now wk1 contains (T+V)|psi>

c	  call prod  (ngorb,psi(i1beg),wk1)
c	  w=ddot(ngorb,wgt1,1,wk1,1)
c	  write(*,*) '<T+V> ',w
c          stop

      do i=1,ngorb
         wk2(i)=0.d0
      enddo

      if (nel.gt.1) then
         if (norb.eq.1) then

c           contribution from coulomb potential of one sigma orbital

	    call dcopy (ngpot,pot(i2beg),1,wk2,1)
	    call prod  (ngorb,psi(i1beg),wk2)
         else
	    do i=1,mxsize
               wk2(i)=0.d0
	    enddo

c    add contributions from coulomb and exchange potentials

	    do iorb1=1,norb
               i1beg1=i1b(iorb1)
               i2beg1=i2b(iorb1)
               ngorb1=i1si(iorb1)		
               ngpot1=i2si(iorb1)
               oc=occ(iorb1)				
               kex=iorb+norb*(iorb1-1)

               if (iorb.le.iorb1) then
		  ihc=iorb+iorb1*(iorb1-1)/2
               else
		  ihc=iorb1+iorb*(iorb-1)/2
               endif

               i3beg=i3b(ihc)
               ngex =i3si(ihc)

               ngrid=min(ngorb,ngpot1)		
               do i=1,mxsize
	          wk0(i)=0.d0
               enddo	    

               call dcopy (ngpot1,pot(i2beg1),1,wk0,1)
               call prod  (ngrid,psi(i1beg),wk0)

               if (iorb.eq.iorb1) oc=oc-1.d0
               if (oc.ne.1.d0) then
		  call dscal (ngrid,oc,wk0,1)
               endif

               if (iorb1.ne.iorb)  then
                  ngrid2=min(ngorb1,ngex,ngorb)		
		  oc=gec(kex)
		  call prodas (ngrid2,-oc,psi(i1beg1),excp(i3beg),wk0)
		  if (ilc(ihc).gt.1) then
                     oc=gec(kex+norb*norb)
                     call prodas (ngrid2,-oc,psi(i1beg1),
     &                    excp(i3beg+ngex),wk0)
		  endif
		else
                   if ((ipe.gt.0).and.(ilc(ihc).gt.0)) then
                      ngrid2=min(ngorb1,ngex,ngorb)		
                      oc=gec(kex)
                      call prodas (ngrid2,-oc,psi(i1beg1),
     &                     excp(i3beg),wk0)
                   endif
		endif

		if (iorb1.eq.1) then
                   call dcopy (ngrid,wk0,1,wk2,1)
		else
                   call add   (ngrid,wk0,wk2)
		endif
             enddo
	  endif

c         to complete the integrand wk2 has to be multiplied by psi(i1beg)

  	  call prod (ngorb,wgt2,wk2)
       endif	  	  

c      wk2 contains V(fock)|psi>
       
       i1beg=i1b(iorb)
       ngorb=i1si(iorb)

       call add   (ngorb,wk1,wk2)
       
       call dcopy (ngorb,wk2,1,wk3,1)

c      call prod  (ngorb,psi(i1beg),wk2)
c      w=ddot(ngorb,wgt1,1,wk2,1)
c      write(*,*) '<T+V> ',w

       w=-eng(iorb)

       call dcopy(ngorb,f4,1,wk1,1)
       call prod(ngorb,psi(i1beg),wk1)
       call prod(ngorb,wgt2,wk1)
       call dcopy (ngorb,wk1,1,wk0,1)
       call daxpy(ngorb,w,wk1,1,wk2,1)

       w=0.d0
       w1=0.d0
       w2=0.d0
       wk2max=0.d0
       do i=1,ngorb

c         exclude grid points where nuclei are located

          if (abs(wk2(i)).gt.wk2max) then
             wk2max=abs(wk2(i))
             imax=i
          endif
          if (.not.(i.eq.1.or.i.eq.(ngrid-mxnmu))) then 
             w=w+wk2(i)*wk2(i)
             w1=w1+wk2(i)*wk2(i)*f4(i)*wgt2(i)
c            w2=w2+wk2(i)*wk2(i)*f4(i)
          endif
       enddo

       write(*,*) 'imax, wk2max',imax,wk2max
       
       if (idbg(496).ne.0) then
          go to ( 11, 12, 13, 14, 15, 16), iorb
00011     fn='bf-1p.lenergy'
          goto 100
00012     fn='bf-5s.lenergy'
          goto 100
00013     fn='bf-4s.lenergy'
          goto 100
00014     fn='bf-3s.lenergy'
          goto 100
00015     fn='bf-2s.lenergy'
	  goto 100
00016     fn='bf-1s.lenergy'
00100     continue
          open(ioutmat,file=fn,status='unknown',form='formatted')
  	  ibeg=i1b(iorb)
          call prtmat (nni,i1mu(1),wk2,ioutmat)
          close(ioutmat)
        endif

c       write(*,1010) iorn(iorb),bond(iorb),w,sqrt(w),w1,w2,sqrt(w2)
c       write(*,1010) iorn(iorb),bond(iorb),w,sqrt(w),w1,sqrt(w1)
        write(*,1010) iorn(iorb),bond(iorb),w1,sqrt(w1),
     &       sqrt(w1)/(-eng(iorb))
c01010   format(10x,i2,1x,a5,'  ',2d15.5/16x,d15.5/16x,2d15.5)
01010   format(10x,i2,1x,a5,'  ',4d15.5)

c        call enrradius(wk2)

        if (idbg(497).ne.0) then
           call leexact1(wk3,wk0)
c          call divide(wk0,psi(i1beg))
        else
c          call leexact(wk1,wk2)
        endif

	return
	end

c Several auxiliary routines follow.

c ### prtmatle ###
c
      subroutine prtmatle (m,n,a,ioutmat)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      dimension a(m,n)

c  loop over grid points

      if (n.gt.500) then
         print *,'prtmatle: nonadequate format'
         stop
      endif

      do im=1,m
         write(ioutmat,1000) (a(im,in),in=1,n)
      enddo
01000 format(500F25.15)

      return
      end


c ### enrradius ###
c
      subroutine enrradius (vt)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      dimension vt(nni,mxnmu)
      include 'commons8.inc'

      dev1=0.d0
      dev2=0.d0
      idev2=0
      r1max=0.2d0
      ipoints=0
      izero=0

      do ini=1,nni
         do imu=1,mxnmu
            r1=(r/2.d0)*(vxi(imu)+veta(ini))            
            t=vt(ini,imu)**2
            dev1=dev1+t
            if (abs(r1).gt.r1max) then
               dev2=dev2+t
               idev2=idev2+1
            endif
         enddo
      enddo
      write(*,*) 'Eucleadian norm ',sqrt(dev1),-0.5d0-eng(1)
      write(*,*) 'Eucleadian norm for r>rmax ',idev2,r1max,sqrt(dev2)

      return
      end

c ### leexact ###
c
      subroutine leexact (a,b)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      dimension a(nni,mxnmu),b(nni,mxnmu)
      include 'commons8.inc'
      iorb=1
      dev1=0.d0
      dev2=0.d0
      dev3=0.d0
      ipoints=0
      izero=0
      do ini=1,nni
         do imu=1,mxnmu

            r1=(r/2.d0)*(vxi(imu)+veta(ini))            

            dnum=0.d0
            ddnom=0.d0
            a(ini,imu)=0.d0
            s1=0.d0
            s2=0.d0
            s3=0.d0

            if (abs(r1).gt.precis) then
               do ipb=npbasis,1,-1
                  d1=primexp(ipb)
                  c1=primcoef(iorb,ipb)
                  e=(2.d0*d1/pii)**(3.d0/4.d0)*exp(-d1*r1*r1)
                  s1=s1+c1*(-2.d0*d1*d1*r1*r1)*e
                  s2=s2+c1*(3.d0*d1)*e
                  if (abs(r1).gt.precis) s3=s3+c1*(-z1/r1)*e
                  ddnom=ddnom+c1*e
               enddo
               dnum=s1+s2+s3
               if (idbg(497).ne.0) then
                  a(ini,imu)=dnum
               else
                  a(ini,imu)=dnum+0.5d0*ddnom
               endif
               ipoints=ipoints+1
               dev1=dev1+(a(ini,imu))**2
            else
               izero=izero+1
            endif
         enddo
      enddo
      write(*,*) 'deviation of gaussian approximation from -Z**2/2'
      write(*,*)  ipoints,sqrt(dev1),sqrt(dev1/dble(ipoints))

      dev2=0.d0
      dev3=0.d0
      do inii=1,nni
         do imu=1,mxnmu
c            dev2=dev2+(a(ini,imu)+z1*z1/2.d0-b(ini,imu))**2
            dev2=dev2+(a(ini,imu)-b(ini,imu))**2
         enddo
      enddo

      write(*,*) 'deviation of gaussian approximation from exact one'
      write(*,*) sqrt(dev2)
      write(*,'(2d20.6,2i10)') sqrt(dev1), sqrt(dev2), ipoints, izero

      stop

      dev1=0.d0
      dev2=0.d0
      dev3=0.d0
      ipoints=0
      do ini=1,nni
         do imu=1,mxnmu
            r1=(r/2.d0)*(vxi(imu)+veta(ini))            
            dnum=0.d0
            ddnom=0.d0
            a(ini,imu)=0.d0
            do ipb=1,npbasis
               ic1 =icgau(ipb)
               d1  =primexp(ipb)
               c1 =primcoef(iorb,ipb)
               if (abs(r1).gt.precis.and.(d1*r1*r1).lt.200.0d0) then
                  dnum=dnum+c1*(-2.d0*d1*d1*r1*r1+3.d0*d1-z1/r1)*
     &                 (2.d0*d1/pii)**(3.d0/4.d0)*exp(-d1*r1*r1)
                  ddnom=ddnom+c1*(2.d0*d1/pii)**(3.d0/4.d0)*
     &                 exp(-d1*r1*r1)
               endif
            enddo
            if (abs(ddnom).gt.precis) then
               a(ini,imu)=dnum+0.5d0*ddnom
            else
               a(ini,imu)=0.d0
            endif
            ipoints=ipoints+1
            dev1=dev1+(a(ini,imu))**2
         enddo
      enddo
      write(*,*) 'deviation of gaussian approximation from -Z**2/2'
      write(*,*)  ipoints,sqrt(dev1),sqrt(dev1/dble(ipoints))

      dev2=0.d0
      dev3=0.d0
      do ini=1,nni
         do imu=1,mxnmu
            dev2=dev2+(a(ini,imu)-b(ini,imu))**2
         enddo
      enddo

      write(*,*) 'deviation of gaussian approximation from exact one'
      write(*,*) sqrt(dev2)
      write(*,'(3d20.6,i10)') sqrt(dev1), sqrt(dev2), sqrt(dev3),ipoints

      return
      end

c ### leexact1 ###
c
      subroutine leexact1 (vt,orb)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      dimension vt(nni,mxnmu),orb(nni,mxnmu)
      include 'commons8.inc'
      iorb=1
      w=eng(iorb)

      dev1=0.d0
      dev2=0.d0
      idev2=0
      r1max=0.2d0
      ipoints=0
      izero=0

      do ini=1,nni
         do imu=1,mxnmu
            r1=(r/2.d0)*(vxi(imu)+veta(ini))            
            dnum=0.d0
            ddnom=0.d0
            s1=0.d0
            s2=0.d0
            s3=0.d0
            if (abs(r1).gt.precis) then
               do ipb=npbasis,1,-1
                  d1=primexp(ipb)
                  c1=primcoef(iorb,ipb)
                  e=(2.d0*d1/pii)**(3.d0/4.d0)*exp(-d1*r1*r1)
                  s1=s1+c1*(-2.d0*d1*d1*r1*r1)*e
                  s2=s2+c1*(3.d0*d1)*e
                  if (abs(r1).gt.precis) s3=s3+c1*(-z1/r1)*e
                  ddnom=ddnom+c1*e
               enddo
               dnum=s1+s2+s3
               if (idbg(497).ne.0.and.abs(ddnom).gt.precis) then
                  t=(vt(ini,imu)-dnum/ddnom*orb(ini,imu))**2
                  dev1=dev1+(vt(ini,imu)-dnum/ddnom*orb(ini,imu))**2
               endif
               if (abs(r1).gt.r1max) then
                  dev2=dev2+t
                  idev2=idev2+1
               endif
            endif
         enddo
      enddo
      write(*,*) 'deviation from exact local energy',sqrt(dev1)
      write(*,*) 'deviation from exact local energy outside radius'
      write(*,*) idev2,r1max,sqrt(dev2)

      return
      end

c ### divide ###
c
      subroutine divide (a,b)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      dimension a(nni,mxnmu),b(nni,mxnmu)
      include 'commons8.inc'
      iorb=1

      dev1=0.d0
      dev2=0.d0
      dev3=0.d0
      ipoints=0
      izero=0
      do ini=1,nni
         do imu=1,mxnmu
            if (abs(b(ini,imu)).gt.precis) then
               a(ini,imu)=a(ini,imu)/b(ini,imu)
            else
               a(ini,imu)=0.d0
               izero=izero+1
            endif
         enddo
      enddo
      write(*,*) 'deviation of gaussian approximation from -Z**2/2'
      write(*,*)  sqrt(dev1),izero

      return
      end

c ### exorbitals ###
c     
C     This subroutine prints out orbitals
c
      subroutine exorbitals (ioutmat,psi)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      character*13 fn
      include 'commons8.inc'

      dimension psi(*)

      open(ioutmat,file='vni',status='unknown',form='formatted')
      write(ioutmat,*) nni
      write(ioutmat,1000) (vni(in),in=1,nni)
      close(ioutmat)

      open(ioutmat,file='vmu',status='unknown',form='formatted')
      write(ioutmat,*) i1mu(1)
      write(ioutmat,1000) (vmu(im),im=1,i1mu(1))
      close(ioutmat)
01000 format(F25.15)

      do iorb=1,norb
         iout=ioutmat+norb-iorb+1
         open(iout,status='unknown',form='formatted')
         ibeg=i1b(iorb)
         call prtmat (nni,i1mu(1),psi(ibeg),iout)
         close(iout)
         write(*,*) 'exorbitals: ',iout
      enddo
01010 format(10x,i2,1x,a5,'  ',4d15.5)
        
      return
      end