********************************************************************************
*                                                                              * 
*  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 ### relorb1 ###  
c
      subroutine relorb1 (iorb,cw_sor,psi,pot,excp,
     &     b,d,e,f0,f1,f2,f4,wgt2,lhs,rhs,wk0,wk1,wk2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      integer*4 cw_sor
      real*8 lhs
      include 'commons8.inc'

      dimension psi(*),pot(*),excp(*),cw_sor(*),
     &     rhs(*),b(*),d(*),e(*),
     &     f0(*),f1(*),f2(*),f4(*),wgt2(*),
     &     lhs(*),wk0(*),wk1(*),wk2(*)
     
      ddmax=abs(demax(1))
      if (alphaf.ne.0.d0) then
         if (sltthr.gt.ddmax.and.ddmax.ne.0.d0) then
	    if(islat.eq.1) then
               write(*,*) '*** the full exchange has been turned on ***'
               islat=0
               exlexp=0.d0
	    endif
         endif
      endif

c     prepare right and left-hand side of the Poisson equation
c     include the diagonal part of the differentiation operator in lhs

      if (islat.ne.0) then
c        the slater exchange approximation is used
         call focksl(iorb,psi,pot,excp,e,f0,f1,f2,f4,wk0,wk1,wk2)
      else
         call fock (iorb,psi,pot,excp,e,f0,f1,f2,f4,wk0,wk1,wk2)
      endif

      call dcopy (mxsize,foc1,1,lhs,1)
      call dcopy (mxsize,foc2,1,rhs,1)

      ipex=mgx(6,iorb)
      ipex=ipex-2*(ipex/2)
      if (ipex.eq.0) then
         isym= 1 
      else
         isym=-1
      endif

c     start the sor iterations
c     fill wk2 array

      iborb=i1b(iorb)

      do itr1=1,maxsor1	  
         do ig=1,i1ng(iorb)
            muoffs=0
            omega=ovforb(ig)
            omega1=1.d0-omega
            ioffst=ioffs(ig)	  	  
            ioffs1=ioffst+1	  
            iborb1=iborb+ioffst
           
            do i=1,4
               dmu2t(i)=dmu2(i,ig)
               dmu1t(i)=dmu1(i,ig)
            enddo

            do i=1,ngsize(ig)
               lhs(ioffst+i)=foc1(ioffst+i)+diag(ig)
            enddo	    

            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)	  	  

c              initialize aorb array needed to evaluate orbital values in
c              the asymptotic region

               nmut=i1mu(iorb)
               do i=nmut-4,nmut
                  ii=(i-1)*nni
                  do j=1,nni
                     ij=ii+j
                     if (f1(ij).eq.0.d0) then
                        wk1(ij)=0.d0
                     else
                        wk1(ij)=lhs(ij)/f1(ij)
                     endif
                  enddo
               enddo
               call iniaorb (nmut,wk0,wk1)
               
               call fill (nni,nmu(ig),isym,psi(iborb),wk2)
               do itr2=1,maxsor2	  		
		  call sor
     &                 (wk2,lhs(ioffs1),rhs(ioffs1),b(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),psi(iborb),wk2)     
               call asymorb (nmut,psi(iborb),wk0)				
            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)	  	  	    
                  call fill2 (nni,nmu(ig),psi(iborb),wk2)
                  do itr2=1,maxsor2	  		
                     call sor 
     &                    (wk2,lhs(ioffs1),rhs(ioffs1),b(ioffs1),
     &                    d(ioffs1),
     &                    cw_sor(iadext(ig+ngrids)),
     &                    cw_sor(iadnor(ig+ngrids)),
     &                    cw_sor(iadex1(ig+ngrids)),
     &                    cw_sor(iadex2(ig+ngrids)),
     &                    cw_sor(iadex3(ig+ngrids)))
                  enddo     
                  call refill (nni,nmu(ig),psi(iborb),wk2)     
               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)	  	  	    
                  
                  muoffs=iemu(ig-1)-1
                  call fill3 (nni,nmu(ig),psi(iborb),wk2)
                  do itr2=1,maxsor2	  		
                     call sor 
     &                    (wk2,lhs(ioffs1),rhs(ioffs1),b(ioffs1),
     &                    d(ioffs1),
     &                    cw_sor(iadext(ig+ngrids)),
     &                    cw_sor(iadnor(ig+ngrids)),
     &                    cw_sor(iadex1(ig+ngrids)),
     &                    cw_sor(iadex2(ig+ngrids)),
     &                    cw_sor(iadex3(ig+ngrids)))
                  enddo     
                  call refill34 (nni,nmu(ig),psi(iborb),wk2)     
                  muoffs=0
               elseif (ig.eq.i1ng(iorb)) then 
                  
c                 initialize aorb array needed to evaluate orbital values in
c                 the asymptotic region
                  
                  nmut=i1mu(iorb)
                  do i=nmut-4,nmut
                     ii=(i-1)*nni
                     do j=1,nni
                        ij=ii+j
                        if (f1(ij).eq.0.d0) then
                           wk1(ij)=0.d0
                        else
                           wk1(ij)=lhs(ij)/f1(ij)
                        endif
                     enddo
                  enddo
                  call iniaorb (nmut,wk0,wk1)
                  
                  ifill=4
                  ngrd1 =ingr1(1,ig)
                  ngrd6a=ingr1(2,ig)
                  ngrd6b=ingr1(3,ig)
                  ngrd7 =ingr1(4,ig)	  	  	    
                  
                  muoffs=iemu(ig-1)-1
                  
                  call fill4 (nni,nmu(ig),psi(iborb),wk2)
                  do itr2=1,maxsor2	  		
                     call sor 
     &                    (wk2,lhs(ioffs1),rhs(ioffs1),b(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 refill34 (nni,nmu(ig),psi(iborb),wk2)     
                  call asymorb (nmut,psi(iborb),wk0)
                  muoffs=0	    
               endif
            endif
         enddo
         
c        if iortho is nonzero the strong orthogonality conditions are imposed
c        i.e. the orthogonalization is performed every maxsor2 iterations
         
         if (iortho.eq.1) then
            iprt=0
            if (iorb.gt.1) call ortho (iorb,psi,f4,wgt2,wk2,iprt)
         endif
      enddo
      
      return
      end