********************************************************************************
*                                                                              *
*  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 ### lagra ###
c
c     Calculates the off-diagonal lagrange multipliers.
c     Damping factors are set to 1, i.e. dmp=1
c 
      subroutine lagra (iorb,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(1),wk0(*),wk1(*)
      dimension engl(3600)

c     there are two possibilities to calculate Lagrange multipliers.
c     if abs((cw1*occ(i)-cw2*occ(iorb))/(c12*occ(i)*occ(iorb))  is
c              > 5.e-02  take <j/h/i> - <i/h/j>  
c              < 5.e-02  take <j/h/i> + <i/h/j>

      if    (idbg(370).eq.0) then
         dmp=1.d0
      elseif(idbg(371).eq.0) then
         dmp=0.8d0
      elseif(idbg(372).eq.0) then
         dmp=0.6d0
      elseif(idbg(373).eq.0) then
         dmp=0.4d0
      elseif(idbg(374).eq.0) then
         dmp=0.2d0
      elseif(idbg(375).eq.0) then
         dmp=0.1d0
      endif

      if (norb.le.1) return
      if (nel.eq.1)  return
      if (iorb.eq.norb) return

      iorbbeg=iorb+1
      do iorb1=1,norb
         if (ifix(iorb1).ne.0) iorbbeg=1
      enddo

      do iorb1=iorbbeg,norb

c        if break is on and the two orbitals have different symmetry
c        off-diagonal lm is not calculated

c        if (ihomon.eq.2.and.ihomo(iorb1)*ihomo(iorb).lt.0) goto 10
         if (ihomo(iorb1)*ihomo(iorb).lt.0) goto 10

         if (icou(iorb1,iorb).eq.0) goto 10

         ipc1=iorb1+(iorb-1)*norb
         ipc2=iorb+(iorb1-1)*norb

         cw1=1.0d0
         cw2=1.0d0
         c12=cw1*cw2
         
         engl(ipc1)=engo(ipc1)
         engl(ipc2)=engo(ipc2)

         call rayln (iorb1,iorb,psi,pot,excp,wgt2,wk0,wk1)
         call rayln (iorb,iorb1,psi,pot,excp,wgt2,wk0,wk1)

         ent=0.d0
         divt=(cw1*occ(iorb1)-cw2*occ(iorb))/(c12*occ(iorb1)*occ(iorb))

         if (abs(divt).lt.5.d-02) then
            divt=(cw1*occ(iorb1)+cw2*occ(iorb))/
     &           (c12*occ(iorb1)*occ(iorb))
            ent=(engo(ipc1)+engo(ipc2))/divt
         else
            ent=(engo(ipc1)-engo(ipc2))/divt
         endif	   	   

         if (idbg(333).ne.0) then
            write(*,*) 'E(ab) ',iorn(iorb),iorn(iorb1),
     &           engo(ipc1),engo(ipc1)/occ(iorb1)
            write(*,*) 'E(ba) ',iorn(iorb1),iorn(iorb),
     &           engo(ipc2),engo(ipc2)/occ(iorb)
         endif
         engo(ipc1)=sflagra*ent/(cw2*occ(iorb))
         engo(ipc2)=sflagra*ent/(cw1*occ(iorb1))

         engo(ipc1)=dmp*ent/(cw2*occ(iorb ))+(1.d0-dmp)*engl(ipc1)
         engo(ipc2)=dmp*ent/(cw1*occ(iorb1))+(1.d0-dmp)*engl(ipc2)

 10      continue
      enddo

      return
      end