********************************************************************************
*                                                                              *
*  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 ### excont ###
c
c     Calculates the number of exchange integrals between two given 
c     (open) shells.

      subroutine excont (iorb1,iorb2,ox1,ox2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      character*8 alpha,beta

      include 'commons8.inc'
      data alpha,beta /'+','-'/

      ipe1 =mgx(6,iorb1)
      ipe2 =mgx(6,iorb2)

      io1=iorb1
      io2=iorb2

      ox1=0.d0
      ox2=0.d0

      ip1=4*(io1-1)
      ip2=4*(io2-1)

      if (ipe1.eq.0) then
         icase=1
         i1end=2
         if(ipe2.eq.0) then
            i2end=2
         else
            i2end=4
         endif
      else
         i1end=4
         if(ipe2.eq.0) then
            icase=1
            i2end=2
         else
            icase=2
         endif
      endif

      if (icase.eq.1) then

c        interaction between sigma-sigma or sigma-nonsigma shells

         do i1=1,i1end
            if(spin(ip1+i1).ne.alpha.and.spin(ip1+i1).ne.beta) goto 10
            do i2=1,i2end
               if(spin(ip2+i2).ne.alpha.and.spin(ip2+i2).ne.beta) 
     &              goto 12
               if(spin(ip1+i1).eq.spin(ip2+i2)) ox1=ox1+1.d0
00012          continue
            enddo
00010       continue
         enddo
      elseif(icase.eq.2) then

c        interaction between nonsigma-nonsigma shells
c        lambda positive for the both orbitals

	   do i1=1,2
              if(spin(ip1+i1).ne.alpha.and.spin(ip1+i1).ne.beta) goto 20
              do i2=1,2
                 if(spin(ip2+i2).ne.alpha.and.spin(ip2+i2).ne.beta) 
     &                goto 22
                 if(spin(ip1+i1).eq.spin(ip2+i2)) ox1=ox1+1.d0
00022            continue
              enddo
00020         continue
           enddo

c          lambda negative for the both orbitals

	   do i1=3,4
              if(spin(ip1+i1).ne.alpha.and.spin(ip1+i1).ne.beta) goto 30
              do i2=3,4
                 if(spin(ip2+i2).ne.alpha.and.spin(ip2+i2).ne.beta) 
     &                goto 32
                 if(spin(ip1+i1).eq.spin(ip2+i2)) ox1=ox1+1.d0
00032            continue
              enddo
00030         continue
           enddo

c          lambda positive and negative

	   do i1=1,2
              if(spin(ip1+i1).ne.alpha.and.spin(ip1+i1).ne.beta) goto 40
              do i2=3,4
                 if(spin(ip2+i2).ne.alpha.and.spin(ip2+i2).ne.beta) 
     &                goto 42
                 if(spin(ip1+i1).eq.spin(ip2+i2)) ox2=ox2+1.d0
00042            continue
              enddo
00040         continue
           enddo

c          lambda negative positive

	   do i1=3,4
              if(spin(ip1+i1).ne.alpha.and.spin(ip1+i1).ne.beta) goto 50
              do i2=1,2
                 if(spin(ip2+i2).ne.alpha.and.spin(ip2+i2).ne.beta) 
     &                goto 52
                 if(spin(ip1+i1).eq.spin(ip2+i2)) ox2=ox2+1.d0
00052            continue
              enddo
00050         continue
           enddo
      endif

      return
      end