********************************************************************************
*                                                                              * 
*  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 ### raynsl ###
c
c     Evaluates the off-diagonal Lagrange multiplier e(iorb2,iorb1) 
c     assuming the local exchange approximation

      subroutine raynsl(iorb1,iorb2,psi,pot,excp,wgt2,wk0,wk1)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)

      include 'commons8.inc'
      dimension psi(*),pot(*),excp(*),wgt2(*),wk0(*),wk1(*)

      if (mgx(6,iorb1).ne.mgx(6,iorb2)) return
      if (ige(iorb1).ne.ige(iorb2)) return

      isd=iorb1+(iorb2-1)*norb
      ipe=mgx(6,iorb2)
      engo(isd)=0.d0
      iborb1=i1b(iorb1)
      isiorb1=i1si(iorb1)	
      iborb2=i1b(iorb2)
      isiorb2=i1si(iorb2)	
      norb2=norb*norb

c     Add contributions from coulomb and local exchange potential.
c     In the local exchange approximtion the coulomb potential includes
c     also the contribution from the orbital

c     calculate the coulomb potential contribution from all the orbitals

      do i=1,mxsize
         wk0(i)=0.d0
      enddo
	
      do iorb=1,norb
         iborb=i1b(iorb)
         ibpot=i2b(iorb)	
         nmut=i1mu(iorb)
         isipot=i2si(iorb)
         oc=occ(iorb)
         call daxpy (isipot,oc,pot(ibpot),1,wk0,1)
      enddo

      call prod  (isiorb2,psi(iborb2),wk0)

c     multiply the local exchange potENTIAL by psi(iborb2) 
c     and add the result to the coulomb potential

      call prod2 (isiorb2,psi(iborb2),excp(1),wk1)
      call add (isiorb2,wk0,wk1)

c     to complete the integrand wk1 has to be multiplied by psi(iborb1)

      ngrid=min(isiorb1,isiorb2)
      call prod (ngrid,psi(iborb1),wk1)
      wtwoel=ddot(ngrid,wgt2,1,wk1,1)

      engo(isd)= wtwoel

      return
      end