BEGIN_PROVIDER [ integer, NSOMOMax]
 &BEGIN_PROVIDER [ integer, NCSFMax]
  implicit none
  BEGIN_DOC
  ! Documentation for NSOMOMax
  ! The maximum number of SOMOs for the current calculation.
  ! required for the calculation of prototype arrays.
  END_DOC
  NSOMOMax = 8
  NCSFMax = 14 ! TODO: NCSFs for MS=0
  END_PROVIDER

  BEGIN_PROVIDER [ double precision, AIJpqMatrixDimsList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)]
  use cfunctions
  implicit none
  BEGIN_DOC
  ! Documentation for AIJpqMatrixList
  ! The prototype matrix containing the <I|E_{pq}|J>
  ! matrices for each I,J somo pair and orb ids.
  END_DOC
  integer i,j,k,l
  integer*8 Isomo, Jsomo
  Isomo = 0
  Jsomo = 0
  integer*8 rows, cols
  rows = -1
  cols = -1
  integer*8 MS
  MS = 0
  print *,"NSOMOMax = ",NSOMOMax
  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
  do i = 2, NSOMOMax, 2
     Isomo = ISHFT(1,i)-1
     do j = i-2,i+2, 2
        Jsomo = ISHFT(1,j)-1
        if(j .GT. NSOMOMax .OR. j .LE. 0) then
           cycle
        end if
        do k = 1,NSOMOMax
           do l = k,NSOMOMax
              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)
              print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
              ! i -> j
              AIJpqMatrixDimsList(i,j,1,k,l,1) = rows
              AIJpqMatrixDimsList(i,j,1,k,l,2) = cols
              AIJpqMatrixDimsList(i,j,1,l,k,1) = rows
              AIJpqMatrixDimsList(i,j,1,l,k,2) = cols
              ! j -> i
              AIJpqMatrixDimsList(j,i,1,k,l,1) = rows
              AIJpqMatrixDimsList(j,i,1,k,l,2) = cols
              AIJpqMatrixDimsList(j,i,1,l,k,1) = rows
              AIJpqMatrixDimsList(j,i,1,l,k,2) = cols
           end do
        end do
     end do
  end do
  END_PROVIDER