qp2/src/csf/sigma_vector.irp.f

 real*8 function logabsgamma(x)
 implicit none
 real*8, intent(in) :: x
 logabsgamma = 1.d32  ! Avoid floating point exception
 if (x>0.d0) then
   logabsgamma = log(abs(gamma(x)))
 endif
 end function logabsgamma

  BEGIN_PROVIDER [ integer, NSOMOMax]
 &BEGIN_PROVIDER [ integer, NSOMOMin]
 &BEGIN_PROVIDER [ integer, NCSFMax]
 &BEGIN_PROVIDER [ integer*8, NMO]
 &BEGIN_PROVIDER [ integer, NBFMax]
 &BEGIN_PROVIDER [ integer, n_CSF]
 &BEGIN_PROVIDER [ integer, maxDetDimPerBF]
  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
  integer MS, ialpha
  MS = elec_alpha_num-elec_beta_num
  NSOMOMax = min(elec_num, cfg_nsomo_max + 2)
  NSOMOMin = max(0,cfg_nsomo_min-2)
  ! Note that here we need NSOMOMax + 2 sizes
  ialpha = (NSOMOMax + MS)/2
  NCSFMax = max(1,nint((binom(NSOMOMax,ialpha)-binom(NSOMOMax,ialpha+1)))) ! TODO: NCSFs for MS=0 (CHECK)
  NBFMax = NCSFMax
  maxDetDimPerBF = max(1,nint((binom(NSOMOMax,ialpha))))
  NMO = n_act_orb
  integer i,j,k,l
  integer startdet,enddet
  integer ncfg,ncfgprev
  integer NSOMO
  integer dimcsfpercfg
  integer detDimperBF
  real*8 :: coeff, binom1, binom2
  integer ncfgpersomo
  real*8, external :: logabsgamma
  detDimperBF = 0
  ! number of cfgs = number of dets for 0 somos
  n_CSF = cfg_seniority_index(NSOMOMin)-1
  ncfgprev = cfg_seniority_index(NSOMOMin)
  !do i = 0-iand(MS,1)+2, NSOMOMax,2
  !!print *," i=",0," dimcsf=",1," ncfg=",ncfgprev, " senor=",cfg_seniority_index(0)
  !!do i = NSOMOMin+2, NSOMOMax,2
  !!   if(cfg_seniority_index(i) .EQ. -1)then
  !!      ncfgpersomo = N_configuration + 1
  !!   else
  !!      ncfgpersomo = cfg_seniority_index(i)
  !!   endif
  !!ncfg = ncfgpersomo - ncfgprev
  !!!detDimperBF = max(1,nint((binom(i,(i+1)/2))))
  !!!dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1))))
  !!n_CSF += ncfg * dimcsfpercfg
  !!!if(cfg_seniority_index(i+2) == -1) EXIT
  !!!if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF
  !!ncfgprev = cfg_seniority_index(i)
  !!print *," i=",i," dimcsf=",dimcsfpercfg," ncfg=",ncfg, " senor=",cfg_seniority_index(i)
  !!enddo
  !!print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration
  n_CSF = 0
  !ncfgprev = cfg_seniority_index(0)
  !ncfgpersomo = ncfgprev
  !do i = iand(MS,1), NSOMOMax-2,2
  !  if(cfg_seniority_index(i) .EQ. -1) then
  !    cycle
  !  endif
  !  if(cfg_seniority_index(i+2) .EQ. -1) then
  !    ncfgpersomo = N_configuration + 1
  !  else
  !    if(cfg_seniority_index(i+2) > ncfgpersomo) then
  !        ncfgpersomo = cfg_seniority_index(i+2)
  !    else
  !      k = 0
  !      do while(cfg_seniority_index(i+2+k) < ncfgpersomo)
  !        k = k + 2
  !        ncfgpersomo = cfg_seniority_index(i+2+k)
  !      enddo
  !    endif
  !  endif
  !  ncfg = ncfgpersomo - ncfgprev
  !  if(i .EQ. 0 .OR. i .EQ. 1) then
  !    dimcsfpercfg = 1
  !  elseif( i .EQ. 3) then
  !    dimcsfpercfg = 2
  !  else
  !    if(iand(MS,1) .EQ. 0) then
  !      dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1))))
  !    else
  !      dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2))))
  !    endif
  !  endif
  !  n_CSF += ncfg * dimcsfpercfg
  !  print *," i=",i," dimcsf=",dimcsfpercfg," ncfg=",ncfg, " senor=",cfg_seniority_index(i)
  !  if(cfg_seniority_index(i+2) > ncfgprev) then
  !    ncfgprev = cfg_seniority_index(i+2)
  !  else
  !    k = 0
  !    do while(cfg_seniority_index(i+2+k) < ncfgprev)
  !      k = k + 2
  !      ncfgprev = cfg_seniority_index(i+2+k)
  !    enddo
  !  endif
  !enddo
  n_CSF = 0
  print *," -9(((((((((((((( NSOMOMin=",NSOMOMin
  ncfgprev = cfg_seniority_index(NSOMOMin) ! can be -1 
  if(ncfgprev.eq.-1)then
    ncfgprev=1
  endif
  do i=NSOMOMin,NSOMOMax+2,2
    !k=0
    !do while((cfg_seniority_index(i+2+k) .eq. -1) .and. (k.le.NSOMOMax))
    !  k=k+2
    !end do
    if(cfg_seniority_index(i).eq.-1)cycle
    if(cfg_seniority_index(i+2).eq.-1)then
      ncfg = N_configuration - ncfgprev + 1
      if(ncfg .eq. 0)then
        ncfg=1
      endif
    else
      ncfg = cfg_seniority_index(i+2) - ncfgprev
    endif
    if(i .EQ. 0 .OR. i .EQ. 1) then
      dimcsfpercfg = 1
    elseif( i .EQ. 3) then
      dimcsfpercfg = 2
    else
      if(iand(MS,1) .EQ. 0) then
        ialpha = (i + MS)/2
        dimcsfpercfg = max(1,nint((binom(i,ialpha)-binom(i,ialpha+1))))
      else
        ialpha = (i + MS)/2
        dimcsfpercfg = max(1,nint((binom(i,ialpha)-binom(i,ialpha+1))))
      endif
    endif
    n_CSF += ncfg*dimcsfpercfg
    print *," i=",i," dimcsf=",dimcsfpercfg," ncfg=",ncfg, " ncfgprev=",ncfgprev, " senor=",cfg_seniority_index(i)
    ncfgprev = cfg_seniority_index(i+2)
  end do
  print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration

END_PROVIDER


subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout)
  use bitmasks
  implicit none
  BEGIN_DOC
  ! Documentation for get_phase_qp_to_cfg
  !
  ! This function converts from (aaaa)(bbbb)
  ! notation to (ab)(ab)(ab)(ab)
  ! notation.
  ! The cfgCI code works in (ab)(ab)(ab)(ab)
  ! notation throughout.
  END_DOC
  integer(bit_kind),intent(in)   :: Ialpha(N_int)
  integer(bit_kind),intent(in)   :: Ibeta(N_int)
  real*8,intent(out)             :: phaseout
  integer(bit_kind)              :: mask, deta(N_int), detb(N_int)
  integer                        :: nbetas
  integer                        :: count, k

  ! Initliaze deta and detb
  deta = Ialpha
  detb = Ibeta

  ! Find how many alpha electrons there are in all the N_ints
  integer :: Na(N_int)
  do k=1,N_int
    Na(k) = popcnt(deta(k))
  enddo

  integer :: shift, ipos, nperm
  phaseout = 1.d0
  do k=1,N_int

    do while(detb(k) /= 0_bit_kind)
      ! Find the lowest beta electron and clear it
      ipos = trailz(detb(k))
      detb(k) = ibclr(detb(k),ipos)

      ! Create a mask will all MOs higher than the beta electron
      mask = not(shiftl(1_bit_kind,ipos + 1) - 1_bit_kind)

      ! Apply the mask to the alpha string to count how many electrons to cross
      nperm = popcnt( iand(mask, deta(k)) )

      ! Count how many alpha electrons are above the beta electron in the other integers
      nperm = nperm + sum(Na(k+1:N_int))
      if (iand(nperm,1) == 1) then
        phaseout  = -phaseout
      endif

    enddo

  enddo
end subroutine get_phase_qp_to_cfg



  BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)]
 &BEGIN_PROVIDER [ real*8, psi_coef_config,  (n_CSF,1)]
 &BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)]
 &BEGIN_PROVIDER [ integer, psi_csf_to_config_data, (n_CSF)]
  use cfunctions
  implicit none
  BEGIN_DOC
  ! Documentation for DetToCSFTransformationMatrix
  ! Provides the matrix of transformatons for the
  ! conversion between determinant to CSF basis (in BFs)
  END_DOC
  integer*8 :: Isomo, Idomo
  integer(bit_kind) :: Ialpha(N_int),Ibeta(N_int)
  integer   :: rows, cols, i, j, k
  integer   :: startdet, enddet, idx
  integer*8 MS, salpha
  integer ndetI
  integer :: getNSOMO
  real*8,dimension(:,:),allocatable    :: tempBuffer
  real*8,dimension(:),allocatable    :: tempCoeff
  real*8  :: norm_det1, phasedet

  integer :: nt


  norm_det1 = 0.d0
  MS = elec_alpha_num - elec_beta_num
  ! initialization
  psi_coef_config = 0.d0
  DetToCSFTransformationMatrix(0,:,:) = 1.d0
  do i = 2-iand(MS,1), NSOMOMax,2
    Isomo = IBSET(0_8, i) - 1_8
    ! rows = Ncsfs
    ! cols = Ndets
    salpha = (i+MS)/2
    bfIcfg = max(1,nint((binom(i,salpha)-binom(i,salpha+1))))
    ndetI =  max(1,nint((binom(i,salpha))))
    !bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1))))
    !ndetI =  max(1,nint((binom(i,(i+1)/2))))

    allocate(tempBuffer(bfIcfg,ndetI))
    call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer)
    DetToCSFTransformationMatrix(i,1:bfIcfg,1:ndetI) =  tempBuffer(1:bfIcfg,1:ndetI)
    deallocate(tempBuffer)
  enddo

  integer s, bfIcfg
  integer countcsf
  countcsf = 0
  integer countdet
  countdet = 0
  integer istate
  istate = 1
  psi_csf_to_config_data(1) = 1
  phasedet = 1.0d0
  call omp_set_max_active_levels(1)
  !$OMP PARALLEL
  !$OMP MASTER
  do i = 1,N_configuration
      startdet = psi_configuration_to_psi_det(1,i)
      enddet = psi_configuration_to_psi_det(2,i)
      ndetI = enddet-startdet+1

      allocate(tempCoeff(ndetI))
      countdet = 1
      do j = startdet, enddet
         idx = psi_configuration_to_psi_det_data(j)
         Ialpha(:) = psi_det(:,1,idx)
         Ibeta(:)  = psi_det(:,2,idx)
         call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet)
         tempCoeff(countdet) = psi_coef(idx, istate)*phasedet
         norm_det1 += tempCoeff(countdet)*tempCoeff(countdet)
         countdet += 1
      enddo

       !print *,"dimcoef=",bfIcfg,norm_det1
       !call printMatrix(tempCoeff,ndetI,1)

      s = 0
      do k=1,N_int
        if (psi_configuration(k,1,i) == 0_bit_kind) cycle
        s = s + popcnt(psi_configuration(k,1,i))
      enddo
      salpha = (s+MS)/2
      bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
      !bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))

      ! perhaps blocking with CFGs of same seniority
      ! can be more efficient
      allocate(tempBuffer(bfIcfg,ndetI))
      tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI)

       call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1,1), size(psi_coef_config,1))
       !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1)

      deallocate(tempCoeff)
      deallocate(tempBuffer)
      psi_config_data(i,1) = countcsf + 1
      do k=1,bfIcfg
        psi_csf_to_config_data(countcsf+k) = i
      enddo
      countcsf += bfIcfg
      psi_config_data(i,2) = countcsf
  enddo
  !$OMP END MASTER
  !$OMP END PARALLEL
  call omp_set_max_active_levels(4)

  END_PROVIDER

  BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (NSOMOMin:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)]
 &BEGIN_PROVIDER [ integer, rowsmax]
 &BEGIN_PROVIDER [ integer, colsmax]
  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, tmpsomo
  Isomo = 0
  Jsomo = 0
  integer rows, cols, nsomoi, nsomoj
  rows = -1
  cols = -1
  integer*8 MS
  MS = elec_alpha_num-elec_beta_num
  nsomomin = elec_alpha_num-elec_beta_num
  rowsmax = 0
  colsmax = 0
  print *,"NSOMOMax = ",NSOMOMax
  print *,"NSOMOMin = ",NSOMOMin
  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
  ! Type
  ! 1. SOMO -> SOMO
  !print *,"Doing SOMO->SOMO"
  AIJpqMatrixDimsList(NSOMOMin,1,1,1,1) = 1
  AIJpqMatrixDimsList(NSOMOMin,1,1,1,2) = 1
  do i = NSOMOMin, NSOMOMax, 2
     Isomo = ISHFT(1_8,i)-1
     do j = i-2,i-2, 2
        Jsomo = ISHFT(1_8,j)-1
        if(j .GT. NSOMOMax .OR. j .LT. 0) then
           cycle
        end if
        do k = 1,i
           do l = 1,i
              ! Define Jsomo
              if(k.NE.l)then
                 Jsomo = IBCLR(Isomo, k-1)
                 Jsomo = IBCLR(Jsomo, l-1)
                 nsomoi = i
                 nsomoj = j
              else
                 Isomo = ISHFT(1_8,i)-1
                 Jsomo = ISHFT(1_8,i)-1
                 nsomoi = i
                 nsomoj = i
              endif

              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)
              !print *, "SOMO->SOMO \t",i,j,k,l,">",Isomo,Jsomo,">",rows, cols
              if(rowsmax .LT. rows) then
                 rowsmax = rows
              end if
              if(colsmax .LT. cols) then
                 colsmax = cols
              end if
              ! i -> j
              AIJpqMatrixDimsList(nsomoi,1,k,l,1) = rows
              AIJpqMatrixDimsList(nsomoi,1,k,l,2) = cols
           end do
        end do
     end do
  end do
  ! Type
  ! 2. DOMO -> VMO
  !print *,"Doing DOMO->VMO"
  AIJpqMatrixDimsList(NSOMOMin,2,1,1,1) = 1
  AIJpqMatrixDimsList(NSOMOMin,2,1,1,2) = 1
  do i = NSOMOMin, NSOMOMax, 2
     Isomo = ISHFT(1_8,i)-1
     tmpsomo = ISHFT(1_8,i+2)-1
     do j = i+2,i+2, 2
        Jsomo = ISHFT(1_8,j)-1
        if(j .GT. NSOMOMax .OR. j .LT. 0) then
           cycle
        end if
        do k = 1,j
           do l = 1,j
              if(k .NE. l) then
                 Isomo = IBCLR(tmpsomo,k-1)
                 Isomo = IBCLR(Isomo,l-1)

                 ! Define Jsomo
                 Jsomo = ISHFT(1_8,j)-1;
                 nsomoi = i
                 nsomoj = j
              else
                 Isomo = ISHFT(1_8,j)-1
                 Isomo = IBCLR(Isomo,1-1)
                 Isomo = IBCLR(Isomo,j-1)
                 Jsomo = ISHFT(1_8,j)-1
                 Isomo = ISHFT(1_8,j)-1
                 nsomoi = j
                 nsomoj = j
              endif

              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)
              !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
              if(rowsmax .LT. rows) then
                 rowsmax = rows
              end if
              if(colsmax .LT. cols) then
                 colsmax = cols
              end if
              ! i -> j
              AIJpqMatrixDimsList(nsomoi,2,k,l,1) = rows
              AIJpqMatrixDimsList(nsomoi,2,k,l,2) = cols
           end do
        end do
     end do
  end do
  ! Type
  ! 3. SOMO -> VMO
  !print *,"Doing SOMO->VMO"
  AIJpqMatrixDimsList(NSOMOMin,3,1,1,1) = 1
  AIJpqMatrixDimsList(NSOMOMin,3,1,1,2) = 1
  do i = NSOMOMin, NSOMOMax, 2
     Isomo = ISHFT(1_8,i)-1
     do j = i,i, 2
        Jsomo = ISHFT(1_8,j)-1
        if(j .GT. NSOMOMax .OR. j .LE. 0) then
           cycle
        end if
        do k = 1,i+1
           do l = 1,i+1
              if(k .NE. l) then
                 Isomo = ISHFT(1_8,i+1)-1
                 Isomo = IBCLR(Isomo,l-1)
                 Jsomo = ISHFT(1_8,j+1)-1
                 Jsomo = IBCLR(Jsomo,k-1)
              else
                 Isomo = ISHFT(1_8,i)-1
                 Jsomo = ISHFT(1_8,j)-1
              endif
              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)
              !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
              if(rowsmax .LT. rows) then
                 rowsmax = rows
              end if
              if(colsmax .LT. cols) then
                 colsmax = cols
              end if
              ! i -> j
              AIJpqMatrixDimsList(i,3,k,l,1) = rows
              AIJpqMatrixDimsList(i,3,k,l,2) = cols
           end do
        end do
     end do
  end do
  ! Type
  ! 4. DOMO -> SOMO
  !print *,"Doing DOMO->SOMO"
  AIJpqMatrixDimsList(NSOMOMin,4,1,1,1) = 1
  AIJpqMatrixDimsList(NSOMOMin,4,1,1,2) = 1
  do i = NSOMOMin, NSOMOMax, 2
     do j = i,i, 2
        if(j .GT. NSOMOMax .OR. j .LE. 0) then
           cycle
        end if
        do k = 1,i+1
           do l = 1,i+1
              if(k .NE. l) then
                 Isomo = ISHFT(1_8,i+1)-1
                 Isomo = IBCLR(Isomo,k-1)
                 Jsomo = ISHFT(1_8,j+1)-1
                 Jsomo = IBCLR(Jsomo,l-1)
              else
                 Isomo = ISHFT(1_8,i)-1
                 Jsomo = ISHFT(1_8,j)-1
              endif
              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)
              !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
              if(rowsmax .LT. rows) then
                 rowsmax = rows
              end if
              if(colsmax .LT. cols) then
                 colsmax = cols
              end if
              ! i -> j
              AIJpqMatrixDimsList(i,4,k,l,1) = rows
              AIJpqMatrixDimsList(i,4,k,l,2) = cols
           end do
        end do
     end do
  end do
  print *,"Rowsmax=",rowsmax," Colsmax=",colsmax
  END_PROVIDER

  BEGIN_PROVIDER [ real*8, AIJpqContainer, (NBFMax,NBFmax,NSOMOMax+1,NSOMOMax+1,4,NSOMOMin:NSOMOMax)]
  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.
  !
  ! Due to the different types of excitations which
  ! include DOMOs and VMOs two prototype DOMOs and two
  ! prototype VMOs are needed. Therefore
  ! the 4th and 5th dimensions are NSOMOMax+4 and NSOMOMax+4
  ! respectively.
  !
  ! The type of excitations are ordered as follows:
  ! Type 1 - SOMO -> SOMO
  ! Type 2 - DOMO -> VMO
  ! Type 3 - SOMO -> VMO
  ! Type 4 - DOMO -> SOMO
  END_DOC
  integer i,j,k,l, orbp, orbq, ri, ci
  orbp = 0
  orbq = 0
  integer*8 Isomo, Jsomo, tmpsomo
  Isomo = 0
  Jsomo = 0
  integer rows, cols, nsomoi, nsomoj
  rows = -1
  cols = -1
  integer*8 MS
  MS = 0
  real*8,dimension(:,:),allocatable :: meMatrix
  integer maxdim

  ! Type
  ! 1. SOMO -> SOMO
  AIJpqContainer = 0.d0
  AIJpqContainer(1,1,1,1,1,NSOMOMin) = 1.0d0
  integer :: rows_old, cols_old
  rows_old = -1
  cols_old = -1
  allocate(meMatrix(1,1))
  do i = NSOMOMin+2, NSOMOMax, 2
     Isomo = ISHFT(1_8,i)-1
     j=i-2
     if(j .GT. NSOMOMax .OR. j .LT. 0) cycle
     nsomoi = i
     do k = 1,i
        orbp = k
        do l = 1,i

           ! Define Jsomo
           if(k .NE. l) then
              Jsomo = IBCLR(Isomo, k-1)
              Jsomo = IBCLR(Jsomo, l-1)
              nsomoj = j
           else
              Isomo = ISHFT(1_8,i)-1
              Jsomo = ISHFT(1_8,i)-1
              nsomoj = i
           endif

           call getApqIJMatrixDims(Isomo,           &
                Jsomo, &
                MS,                       &
                rows,                     &
                cols)

           orbq = l
           if ((rows /= rows_old).or.(cols /= cols_old)) then
             deallocate(meMatrix)
             allocate(meMatrix(rows,cols))
             rows_old = rows
             cols_old = cols
           endif
           meMatrix = 0.0d0
           ! fill matrix
           call getApqIJMatrixDriver(Isomo,           &
                Jsomo, &
                orbp,                     &
                orbq,                     &
                MS,                       &
                NMO,                      &
                meMatrix,                 &
                rows,                     &
                cols)
           ! i -> j
          do ri = 1,rows
              do ci = 1,cols
                 AIJpqContainer(ri,ci,k,l,1,nsomoi) = meMatrix(ri, ci)
              end do
           end do
        end do
     end do
  end do
  deallocate(meMatrix)

  ! Type
  ! 2. DOMO -> VMO
  !print *,"Doing DOMO -> VMO"
  !AIJpqContainer(NSOMOMin,2,1,1,1,1) = 1.0d0
  AIJpqContainer(1,1,1,1,2,NSOMOMin) = 1.0d0
  do i = NSOMOMin, NSOMOMax, 2
     Isomo = ISHFT(1_8,i)-1
     tmpsomo = ISHFT(1_8,i+2)-1
     do j = i+2,i+2, 2
        if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
        Jsomo = ISHFT(1_8,j)-1
        do k = 1,j
           do l = 1,j
              if(k .NE. l) then
                 Isomo = IBCLR(tmpsomo,k-1)
                 Isomo = IBCLR(Isomo,l-1)
                 ! Define Jsomo
                 Jsomo = ISHFT(1_8,j)-1;
                 nsomoi = i
                 nsomoj = j
              else
                 Isomo = ISHFT(1_8,j)-1
                 Isomo = IBCLR(Isomo,1-1)
                 Isomo = IBCLR(Isomo,j-1)
                 Jsomo = ISHFT(1_8,j)-1
                 Isomo = ISHFT(1_8,j)-1
                 nsomoi = j
                 nsomoj = j
              endif

              !print *,"k,l=",k,l
              !call debug_spindet(Jsomo,1)
              !call debug_spindet(Isomo,1)

              !AIJpqContainer(nsomoi,2,k,l,:,:) = 0.0d0
              AIJpqContainer(:,:,k,l,2,nsomoi) = 0.0d0
              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)

              orbp = k
              orbq = l
              allocate(meMatrix(rows,cols))
              meMatrix = 0.0d0
              ! fill matrix
              call getApqIJMatrixDriver(Isomo,           &
                   Jsomo, &
                   orbp,                     &
                   orbq,                     &
                   MS,                       &
                   NMO,                      &
                   meMatrix,                 &
                   rows,                     &
                   cols)
             !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax
             !call printMatrix(meMatrix,rows,cols)
              ! i -> j
             do ri = 1,rows
                 do ci = 1,cols
                    !AIJpqContainer(nsomoi,2,k,l,ri,ci) = meMatrix(ri, ci)
                    AIJpqContainer(ri,ci,k,l,2,nsomoi) = meMatrix(ri, ci)
                 end do
              end do
              deallocate(meMatrix)
           end do
        end do
     end do
  end do
  ! Type
  ! 3. SOMO -> VMO
  !print *,"Doing SOMO -> VMO"
  !AIJpqContainer(NSOMOMin,3,1,1,1,1) = 1.0d0
  AIJpqContainer(1,1,1,1,3,NSOMOMin) = 1.0d0
  do i = NSOMOMin, NSOMOMax, 2
     Isomo = ISHFT(1_8,i)-1
     do j = i,i, 2
        Jsomo = ISHFT(1_8,j)-1
        if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
        do k = 1,i+1
           do l = 1,i+1
              if(k .NE. l) then
                 Isomo = ISHFT(1_8,i+1)-1
                 Isomo = IBCLR(Isomo,l-1)
                 Jsomo = ISHFT(1_8,j+1)-1
                 Jsomo = IBCLR(Jsomo,k-1)
              else
                 Isomo = ISHFT(1_8,i)-1
                 Jsomo = ISHFT(1_8,j)-1
              endif

              !print *,"k,l=",k,l
              !call debug_spindet(Jsomo,1)
              !call debug_spindet(Isomo,1)

              !AIJpqContainer(i,3,k,l,:,:) = 0.0d0
              AIJpqContainer(:,:,k,l,3,i) = 0.0d0
              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)

              orbp = k
              orbq = l
              allocate(meMatrix(rows,cols))
              meMatrix = 0.0d0
              ! fill matrix
              call getApqIJMatrixDriver(Isomo,           &
                   Jsomo, &
                   orbp,                     &
                   orbq,                     &
                   MS,                       &
                   NMO,                      &
                   meMatrix,                 &
                   rows,                     &
                   cols)
             !call printMatrix(meMatrix,rows,cols)
             !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax
              ! i -> j
             do ri = 1,rows
                 do ci = 1,cols
                    !AIJpqContainer(i,3,k,l,ri,ci) = meMatrix(ri, ci)
                    AIJpqContainer(ri,ci,k,l,3,i) = meMatrix(ri, ci)
                 end do
              end do
              deallocate(meMatrix)
           end do
        end do
     end do
  end do
  ! Type
  ! 4. DOMO -> SOMO
  !print *,"Doing DOMO -> SOMO"
  !AIJpqContainer(NSOMOMin,4,1,1,1,1) = 1.0d0
  AIJpqContainer(1,1,1,1,4,NSOMOMin) = 1.0d0
  do i = NSOMOMin+2, NSOMOMax, 2
     Isomo = ISHFT(1_8,i)-1
     do j = i,i, 2
        Jsomo = ISHFT(1_8,i)-1
        if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
        do k = 1,i+1
           do l = 1,i+1
              if(k .NE. l) then
                 Isomo = ISHFT(1_8,i+1)-1
                 Isomo = IBCLR(Isomo,k-1)
                 Jsomo = ISHFT(1_8,j+1)-1
                 Jsomo = IBCLR(Jsomo,l-1)
              else
                 Isomo = ISHFT(1_8,i)-1
                 Jsomo = ISHFT(1_8,j)-1
              endif

              !AIJpqContainer(i,4,k,l,:,:) = 0.0d0
              AIJpqContainer(:,:,k,l,4,i) = 0.0d0
              call getApqIJMatrixDims(Isomo,           &
                   Jsomo, &
                   MS,                       &
                   rows,                     &
                   cols)

              orbp = k
              orbq = l

              allocate(meMatrix(rows,cols))
              meMatrix = 0.0d0
              ! fill matrix
              call getApqIJMatrixDriver(Isomo,           &
                   Jsomo, &
                   orbp,                     &
                   orbq,                     &
                   MS,                       &
                   NMO,                      &
                   meMatrix,                 &
                   rows,                     &
                   cols)
             !call printMatrix(meMatrix,rows,cols)
             !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax
              ! i -> j
             do ri = 1,rows
                 do ci = 1,cols
                    !AIJpqContainer(i,4,k,l,ri,ci) = meMatrix(ri, ci)
                    AIJpqContainer(ri,ci,k,l,4,i) = meMatrix(ri, ci)
                 end do
              end do
              deallocate(meMatrix)
           end do
        end do
     end do
  end do
  END_PROVIDER

subroutine calculate_preconditioner_cfg(diag_energies)
  implicit none
  use bitmasks
  BEGIN_DOC
  ! Documentation for calculate_preconditioner
  !
  ! Calculates the diagonal energies of
  ! the configurations in psi_configuration
  ! returns : diag_energies :
  END_DOC
  integer :: i,j,k,kk,l,p,q,noccp,noccq, ii, jj
  real*8,intent(out) :: diag_energies(n_CSF)
  integer                            :: nholes
  integer                            :: nvmos
  integer                            :: listvmos(mo_num)
  integer                            :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
  integer                            :: listholes(mo_num)
  integer                            :: holetype(mo_num) ! 1-> SOMO 2->DOMO
  integer*8                          :: Idomo
  integer*8                          :: Isomo
  integer*8                          :: Jdomo
  integer*8                          :: Jsomo
  integer*8                          :: diffSOMO
  integer*8                          :: diffDOMO
  integer                            :: NSOMOI
  integer                            :: NSOMOJ
  integer                            :: ndiffSOMO
  integer                            :: ndiffDOMO
  integer                            :: starti, endi, cnti, cntj, rows,cols
  integer                            :: extype,pmodel,qmodel
  integer(bit_kind) :: Icfg(N_INT,2)
  integer(bit_kind) :: Jcfg(N_INT,2)
  integer,external  :: getNSOMO
  real*8, external  :: mo_two_e_integral
  real*8            :: hpp
  real*8            :: meCC
  real*8            :: ecore
  real*8            :: core_act_contrib

  !PROVIDE h_core_ri
  PROVIDE core_fock_operator
  PROVIDE h_act_ri
  ! initialize energies
  diag_energies = 0.d0
  !print *,"Core energy=",core_energy," nucler rep=",nuclear_repulsion, " n_core_orb=",n_core_orb," n_act_orb=",n_act_orb," mo_num=",mo_num

  ! calculate core energy
  !call get_core_energy(ecore)
  diag_energies = core_energy - nuclear_repulsion

  ! calculate the core energy
  !print *,"Core 2energy=",ref_bitmask_energy

  do i=1,N_configuration

     Isomo = psi_configuration(1,1,i)
     Idomo = psi_configuration(1,2,i)
     Icfg(1,1) = psi_configuration(1,1,i)
     Icfg(1,2) = psi_configuration(1,2,i)
     NSOMOI = getNSOMO(psi_configuration(:,:,i))

     starti = psi_config_data(i,1)
     endi   = psi_config_data(i,2)

     core_act_contrib = 0.0d0

     ! find out all pq holes possible
     nholes = 0
     ! holes in SOMO
     !do k = 1,mo_num
     do kk = 1,n_act_orb
       k = list_act(kk)
        if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
           nholes += 1
           listholes(nholes) = k
           holetype(nholes) = 1
        endif
     enddo
     ! holes in DOMO
     !do k = n_core_orb+1,n_core_orb + n_act_orb