********************************************************************************
*                                                                              *
*  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 ### oshell ###
c
c     Calculates the weights of the two-electron contributions to the total 
c     energy expression for an open shell scf case.
c
      subroutine oshell
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      character*8 alpha,beta
      include 'commons8.inc'
      data alpha/'+'/, beta/'-'/

c !!! isum=norb

      do i=1,2*maxorb*maxorb
         gec(i)=0.d0
      enddo

      do i=1,4*norb
         spn(i)=spin(i)
      enddo

      do i=1,isum
         isumi=4*(i-1)
*        gec(i)=occ(i)-1.0
         iloop(i)=2
         iqm(1+isumi)=1
         iqm(2+isumi)=1
         if (mgx(3,i).gt.0) then
            iloop(i)=4
            iqm(3+isumi)=-1
            iqm(4+isumi)=-1
         endif
      enddo
      
c     initialize spn and spin arrays
c     throw out the electrons + = alpha and - = beta
c     the orbital is locked = 1 (open shell case)
c                           = 0 (close shell case) 
      
      do i=1,isum
         if (lock(i).eq.0) then
            jj=4*(i-1)
            spn(jj+1)=alpha
            spn(jj+2)=beta
            if (mgx(3,i).gt.0) then
               spn(jj+3)=alpha
               spn(jj+4)=beta
            endif
         endif
      enddo
      
      do i=1,4*isum
         spin(i)=spn(i)
      enddo
      
c     loop now through all electrons
      
      isu2=isum*isum
      do i=1,isum
         ixx=4*(i-1)
         div(i)=0.d0
         do ii=1,iloop(i)
            ix=ii+ixx
            if ((spn(ix).ne.alpha).and.(spn(ix).ne.beta)) goto 12
            div(i)=div(i)+1.d0
            
            do j=1,isum
               jxx=4*(j-1)
               kk=i+isum*(j-1)

               mpl=iabs(mgx(3,i)-mgx(3,j))
               mmi=iabs(mgx(3,i)+mgx(3,j))

               do jj=1,iloop(j)
                  jx=jxx+jj
                  if ((spn(jx).ne.alpha).and.(spn(jx).ne.beta)) goto 15
c                 div(i)=div(i)+1.0d0
                  
                  if ((mpl.eq.mmi).and.(spn(ix).eq.spn(jx)).
     &                 and.(i.ne.j)) then
                     gec(kk)=gec(kk)+1.d0
                  endif
                  
                  if ((mpl.ne.mmi).and.(spn(ix).eq.spn(jx))) then
                     if ((i.eq.j).and.(iqm(ix).ne.iqm(jx))) then
                        gec(kk)=gec(kk)+1.d0
                     endif

                     if ((i.ne.j).and.(iqm(ix).eq.iqm(jx))) then
                        gec(kk)=gec(kk)+1.d0
                     endif
                     
                     if ((i.ne.j).and.(iqm(ix).ne.iqm(jx))) then
                        gec(kk+isu2)=gec(kk+isu2)+1.d0
                     endif
                  endif
 15               continue
               enddo
            enddo
 12         continue
         enddo
         
         if (abs(div(i)).lt.1.d-06) goto 1000
         
         if (idbg(150).ne.0) then
            write(6,1101)
01101       format(2x,' i,idiv,gec(kk),gec(kk+isu2),div(i)')
         endif
         do idiv=1,isum
            kk=i+isum*(idiv-1)
            gec(kk)=gec(kk)/div(i)
            gec(kk+isu2)=gec(kk+isu2)/div(i)
            if (idbg(151).ne.0) then
               write(6,1102)  iorn(i),bond(i),gut(i),
     &              iorn(idiv),bond(idiv),gut(idiv),
     &              gec(kk),gec(kk+isu2),div(i)
01102          format(i4,1x,a8,a1,i4,1x,a8,a1,3f5.1)
            endif
         enddo
 10      continue
      enddo
            
      ial=0
      ibe=0
      imag=0
      do i=1,isum
         isumi=4*(i-1)
         if (lock(i).eq.0) goto 100

         do ii=1,iloop(i)
            imu=ii+isumi
            if (spn(imu).eq.alpha) ial=ial+1
            if (spn(imu).eq.beta) ibe=ibe+1
            if (spn(imu).eq.alpha.or.spn(imu).eq.beta) then
               imag=imag+iqm(imu)
            endif
 110        continue
         enddo
 100     continue
      enddo

c     loop over occupied shells
      
c     no distinction is made between close and open shell cases
c     as far as storage reservation is concerned.
c     count has been modified to cover the both cases.

      do j=1,(isum-1)
         do i=(j+1),isum
            icou(j,i)=0
            icou(i,j)=0

            if (mgx(6,j).ne.mgx(6,i)) goto 300
            if (gut(j).ne.gut(i)) goto 300
            
            if ((lock(j).gt.0).or.(lock(i).gt.0)) then
               icou(j,i)=1
               icou(i,j)=1
            endif
 300        continue
         enddo
      enddo

      return
 1000 continue
      write(iout6,1001)
 1001 format(/1x,'... error in spin orbitals ...'//)
      stop 'oshell'
      end