********************************************************************************
*                                                                              * 
*  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 ### raylsl ###
c 
c     Solves the eigenvalue  to the Schroedinger equation assuming 
c     the local form of the exchange potential
c
      subroutine raylsl(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(*)

      do i=1,mxsize
         wk0(i)=0.d0	  
         wk1(i)=0.d0
      enddo

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

      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
      endif

      if (nel.gt.1) then
         do i=1,mxsize
            wk2(i)=0.d0
         enddo
        
         if (norb.eq.1) then

c           calculate the Coulomb potential contribution 
c           from one sigma orbital

            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 local exchange potential.
c           In the local exchange approximation the coulomb potential 
c           includes also the contribution from the orbital

            do iorb1=1,norb
               i1beg1=i1b(iorb1)
               i2beg1=i2b(iorb1)
               ngorb1=i1si(iorb1)		
               ngpot1=i2si(iorb1)
               oc=occ(iorb1)
               call daxpy (ngpot1,oc,pot(i2beg1),1,wk2,1)
            enddo
            call prod  (ngorb,psi(i1beg),wk2)
            
c           multiply the local exchange potential by psi(i1beg) and 
c           add the result to the coulomb potential (multiplied by
c           psi(i1beg))

	    call prod2 (ngorb,psi(i1beg),excp(1),wk1)
	    call add (ngorb,wk1,wk2)
         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(*,*) 'iorb,woneel,wtwoel', iorb,woneel,wtwoel
      endif

      if(idbg(96).ne.0) then
         write(*,*) 'raylsl: wk1'
         call pmtx (nni,nmu(1),wk1,1,1,incrni,incrmu)
         write(*,*) 'raylsl: wk2'
         call pmtx (nni,nmu(1),wk2,1,1,incrni,incrmu)
      endif

      if (idbg(70).ne.0) then
         write(*,*) 'raylsl: woneel, wtwoel, eng ',woneel,wtwoel,
     &        eng(iorb)
      endif
      
      return
      end