********************************************************************************
*                                                                              *  
*  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 ### raylloc ###
c
c     Solves the eigenvalue to the schroedinger equation as e = <w(i)/h/w(j)>
c
      subroutine raylloc(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,1000) (vni(in),in=1,nni)
      close(ioutmat)

      open(ioutmat,file='vmu',status='unknown',form='formatted')
      write(ioutmat,1000) (vmu(im),im=1,i1mu(1))
      close(ioutmat)
01000 format(5F15.6)

      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 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

      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)     
      
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 daxpy(ngorb,w,wk1,1,wk2,1)

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

      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
      write(*,1010) iorn(iorb),bond(iorb),sqrt(w)
01010 format(18x,i2,1x,a5,'  ',d10.2)

      return
      end