********************************************************************************
*                                                                              *
*  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 ### rayln ###
c 
c     Evaluates the off-diagonal Lagrange multiplier e(iorb2,iorb1)
c
      subroutine rayln(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(*)

c     the Slater exchange approximation is used

      if (islat.ne.0) then
         call raynsl(iorb1,iorb2,psi,pot,excp,wgt2,wk0,wk1)
         return
      endif

      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
      i1beg1=i1b(iorb1)
      i1beg2=i1b(iorb2)
      ngorb1=i1si(iorb1)
      ngpot1=i2si(iorb1)	

      ngorb2=i1si(iorb2)     
      ngpot2=i2si(iorb2)	

      norb2=norb*norb

c     add contributions from coulomb and exchange potentials

      do iorb=1,norb
         i1beg=i1b(iorb)
         ngorb=i1si(iorb)		
         ngpot=i2si(iorb)
         
         kex=iorb2+norb*(iorb-1)
         hh2=1.d0

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

         i2beg=i2b(iorb)
         i3beg=i3b(ihc)
         ngex =i3si(ihc)
         
         ngrid=min(ngorb2,ngpot)		
         do i=1,mxsize
            wk0(i)=0.d0
         enddo	    

         coo=occ(iorb)
         if (iorb2.eq.iorb) coo=coo-1.d0
         
         call dcopy (ngrid,pot(i2beg),1,wk0,1)
         call prod  (ngrid,psi(i1beg2),wk0)

         if (coo.ne.1.d0) then
            call dscal (ngrid,coo,wk0,1)
         endif

         ngrid2=min(ngorb2,ngex,ngorb)		
         if (iorb.ne.iorb2)  then
            coo=hh2*gec(kex)
            call prodas (ngrid2,-coo,psi(i1beg),excp(i3beg),wk0)
            if (ilc(ihc).gt.1) then
               coo=hh2*gec(kex+norb*norb)
               call prodas (ngrid2,-coo,psi(i1beg),
     &              excp(i3beg+i3si(ihc)),wk0)
            endif
         else
            if ((ipe.gt.0).and.(ilc(ihc).gt.0)) then
               coo=hh2*gec(kex)
               call prodas (ngrid2,-coo,psi(i1beg),excp(i3beg),wk0)
            endif
         endif
       
         if (iorb.eq.1) then
            call dcopy (ngrid,wk0,1,wk1,1)
         else
            call add   (ngrid,wk0,wk1)
         endif
      enddo
      
c     to complete the integrand wk1 has to be multiplied by psi(i1beg1)

      call prod (ngrid,psi(i1beg1),wk1)
      wtwoel=ddot(ngrid,wgt2,1,wk1,1)
      
      engo(isd)= wtwoel
      
      return
      end