********************************************************************************
*                                                                              *
*  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 ### rayl ###  
c
c     Solves the eigenvalue to the schroedinger equation as e = <w(i)/h/w(j)>

      subroutine rayl(iorb,psi,pot,excp,e,f0,wgt1,wgt2,
     &			wk0,wk1,wk2,wk3)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'
      dimension psi(*),pot(*),excp(*),e(*),f0(*),wgt1(*),wgt2(*),
     &     wk0(*),wk1(*),wk2(*),wk3(*)

c  Slater exchange approximation is used	

      if (islat.ne.0) then
         call raylsl (iorb,psi,pot,excp,e,f0,wgt1,wgt2,wk0,wk1,wk2,wk3)
         return
      endif

      wtwoel=0.d0
      
      i1beg=i1b(iorb)
      i2beg=i2b(iorb)	
      nmut=i1mu(iorb)
      ngorb=i1si(iorb)
      ngpot=i2si(iorb)	

      if (ngorb.ne.ngpot) then
         write(*,*) 'Orbitals and corresponding Coulomb potentials ',
     &        'have to be defined on the same number of grid points'
         stop 'raylsl'
      endif

      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 (idbg(70).ne.0) then
         call dcopy (ngorb,wk1,1,wk2,1)
         call prod  (ngorb,psi(i1beg),wk2)
         w=ddot(ngorb,wgt1,1,wk2,1)
c         write(*,'(i6,2d22.14)') ngorb,wk1(ngorb),wgt1(ngorb)
         write(*,'(" rayl: kinetic energy ",d22.14)') w
      endif
      
      if (ipe.eq.0) then
         call dcopy (ngorb,f0,1,wk0,1)
      else
         
c        e enters 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

      if (idbg(70).ne.0) then
         call prod2 (ngorb,psi(i1beg),wk0,wk2)
         call prod  (ngorb,psi(i1beg),wk2)
         w=ddot(ngorb,wgt1,1,wk2,1)
c        nuclear energy for non-sigma orbitals contains contribution from e term
c        (in toten this term is correctly added to the kinetic energy contribution)
         write(*,*) ' rayl: nuclear energy ',w
      endif

      call proda (ngorb,psi(i1beg),wk0,wk1)
      call prod (ngorb,psi(i1beg),wk1)
      
      woneel=ddot(ngorb,wgt1,1,wk1,1)
      
      if(idbg(96).ne.0) then
         write(*,*) 'woneel', woneel
         write(*,*) 'rayl: wk0'
         call pmtx (nni,nmu(1),wk0,1,1,incrni,incrmu)
         write(*,*) 'rayl: wk1'
         call pmtx (nni,nmu(1),wk1,1,1,incrni,incrmu)
      endif

      if (nel.gt.1) then
         do i=1,mxsize
            wk2(i)=0.d0
         enddo
         
         if (norb.eq.1) then
            
c           contribution from coulomb potential of one sigma orbital
 
c           14/12/02
c	    call dcopy (ngpot,pot(i2beg),1,wk2,1)
            oc=occ(norb)
            call daxpy (ngpot,oc,pot(i2beg),1,wk2,1)
	    call prod  (ngorb,psi(i1beg),wk2)
         else

c           add contributions from coulomb and exchange potentials

	    do iorb1=1,norb

c              it is asumed that exchange potentials involving iorb1 
c              has already been retrieved from disk

               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,psi(i1beg),wk2)
         wtwoel=ddot(ngorb,wgt2,1,wk2,1)
      endif	  	  

      engo(isd) = woneel+wtwoel
      eng(iorb) = engo(isd)

      if(idbg(95).ne.0) then
         write(*,'(i6,3d22.14)') 'iorb,woneel,wtwoel', iorb,woneel,wtwoel
      endif

      if(idbg(96).ne.0) then
         write(*,*) 'wtwoel', wtwoel
         write(*,*) 'rayl: wk0'
         call pmtx (nni,nmu(1),wk0,1,1,incrni,incrmu)
         write(*,*) 'rayl: wk2'
         call pmtx (nni,nmu(1),wk2,1,1,incrni,incrmu)
      endif
      
      if (idbg(70).ne.0) then
         write(*,'("rayl: woneel, wtwoel, eng ",3d22.14)') 
     &        woneel,wtwoel,eng(iorb)
      endif
      
      return
      end