********************************************************************************
*                                                                              * 
*  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     ### relexch2 ###
c
      subroutine relexch2 (iorb,cw_sor,psi,pot,excp,bpot,d,e,f3,g,
     &			 lhs,rhs,wk2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      real*8 lhs
      integer*4 cw_sor
      include 'commons8.inc'
      dimension psi(*),pot(*),excp(*),cw_sor(*),
     &		rhs(*),bpot(*),d(*),e(*),f3(*),
     &		g(*),lhs(*),wk2(*)	

      if (nel.lt.2 ) return
      
      if     (exlexp.eq.0) then
	  ifirst=1
      elseif (exlexp.eq.2) then
c        recalculate only those exchange potentials which depend on orbitals 
c        modified	so far (note the reverse order of relaxation in scf)
         ifirst=iorb
      endif
      
      ngorb=i1si(iorb)	  
      do iorb1=ifirst,norb
         iswtch=0

         if (iorb.eq.iorb1.and.mgx(6,iorb).eq.0 ) goto 10
         if ((iorb.eq.iorb1).and.(ilc(iorb*(iorb+1)/2).lt.1)) goto 10

c        orbitals in increasing order

         if (iorb.lt.iorb1) then
	    in1=iorb
	    in2=iorb1	  
         else	    
	    in1=iorb1
	    in2=iorb
         endif	    

         ib1=i1b(in1)
         ng1=i1si(in1)
         ib2=i1b(in2)
         ng2=i1si(in2)	  

         ipc=in1+in2*(in2-1)/2
         iwexch(ipc)=-1

         ibexp=i3b(ipc)
         ngexp=i3si(ipc)

         idel=iabs(mgx(6,in1)-mgx(6,in2))
         if (in1.eq.in2) idel=2*mgx(6,in1)
         idel1=mgx(6,in1)+mgx(6,in2)
         
c        prepare right-hand side of the poisson equation
         
         do i=1,ngorb
            wk2(i)=0.d0
         enddo	
         
         ngrid=min(ng1,ng2)
         call prod2 (ngrid,psi(ib1),psi(ib2),wk2)
         call prod2 (ngorb,wk2,g,rhs)
         
         ipex=idel-2*(idel/2)
c        if (mod(idel,2).eq.0) then
         if (ipex.eq.0) then	  
            isym= 1
         else
            isym=-1
         endif
 1000    continue
         do itr1=1,maxsor1	  	
            do ig=1,i1ng(iorb)
               do i=1,i1si(iorb)
                  wk2(i)=0.d0
               enddo	
               
               omega=ovfexch(ig)
               omega1=1.d0-omega
               ioffst=ioffs(ig)	    
               ioffs1=ioffst+1
               ibexp1=ibexp+ioffst	  	  
               
               do i=1,4
                  dmu2t(i)=dmu2(i,ig)
                  dmu1t(i)=dmu1(i,ig)
               enddo
               
c              prepare left-hand side of the poisson equation
c              include the diagonal part of the differentiation 
c              operator in lhs
               
               if (idel.eq.0) then
                  do i=1,ngsize(ig)
                     lhs(ioffst+i)=f3(ioffst+i)+diag(ig)
                  enddo		
               else
                  xdel=dble(idel*idel)
                  do i=1,ngsize(ig)
                     lhs(ioffst+i)=f3(ioffst+i)+
     &                    xdel*e(ioffst+i)+diag(ig)
                  enddo				
               endif
               
               if (i1ng(iorb).eq.1) then 
c                 icase=1
                  ifill=1
                  ngrd1 =ingr1(1,ig)
                  ngrd6a=ingr1(2,ig)
                  ngrd6b=ingr1(3,ig)
                  ngrd7 =ingr1(4,ig)	  	  
                  call fill (nni,nmu(ig),isym,excp(ibexp),wk2)
                  do itr2=1,maxsor2/lsor	  		
                     call mcsor
     &                    (wk2,lhs(ioffs1),rhs(ioffs1),
     &                    bpot(ioffs1),d(ioffs1),
     &                    cw_sor(iadext(ig)),cw_sor(iadnor(ig)),
     &                    cw_sor(iadex1(ig)),cw_sor(iadex2(ig)),
     &                    cw_sor(iadex3(ig)))
                  enddo     
                  call refill (nni,nmu(ig),excp(ibexp),wk2)     
               else
                  if (ig.eq.1) then 
c                    icase=2
                     ifill=2
                     ngrd1 =ingr2(1,ig)
                     ngrd6a=ingr2(2,ig)
                     ngrd6b=ingr2(3,ig)
                     ngrd7 =ingr2(4,ig)	  	  	    
                  elseif (ig.ne.1.and.ig.ne.i1ng(iorb)) then 
c                    icase=2
                     ifill=3
                     ngrd1 =ingr2(1,ig)
                     ngrd6a=ingr2(2,ig)
                     ngrd6b=ingr2(3,ig)
                     ngrd7 =ingr2(4,ig)	  	  	    
                  elseif (ig.eq.i1ng(iorb)) then 
                     ifill=4
                     ngrd1 =ingr1(1,ig)
                     ngrd6a=ingr1(2,ig)
                     ngrd6b=ingr1(3,ig)
                     ngrd7 =ingr1(4,ig)	  	  	    
                     muoffs=iemu(ig-1)-1
                  endif
               endif
            enddo
         enddo
         if (iswtch.eq.1) goto 10
         if (ilc(ipc).lt.2) goto 10
         iwexch(ipc)=-2
         idel=idel1
         ibexp=ibexp+ngexp
         ipc=ipc+norb*(norb+1)/2
         iswtch=1
         goto 1000
 10      continue
      enddo
 1021 format(/,i5,10d12.3,/(5x,10d12.3))
      return
      end