Added a naive function that stores the \alphas.

src/csf/configuration_CI_sigma_helpers.irp.f

@@ -1,3 +1,268 @@
+use bitmasks
+
+ BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,N_configuration*mo_num*mo_num)]
+&BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ]
+  implicit none
+  !use bitmasks
+  BEGIN_DOC
+  ! Documentation for alphasI
+  ! Returns the associated alpha's for
+  ! the input configuration Icfg.
+  END_DOC
+
+  integer                            :: idxI ! The id of the Ith CFG
+  integer(bit_kind)                  :: Icfg(N_int,2)
+  integer                            :: NalphaIcfg
+  integer(bit_kind)                  :: alphasIcfg(N_int,2,N_configuration*mo_num*mo_num)
+  logical,dimension(:,:),allocatable :: tableUniqueAlphas
+  integer                            :: listholes(mo_num)
+  integer                            :: holetype(mo_num) ! 1-> SOMO 2->DOMO
+  integer                            :: nholes
+  integer                            :: nvmos
+  integer                            :: listvmos(mo_num)
+  integer                            :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
+  integer*8                          :: Idomo
+  integer*8                          :: Isomo
+  integer*8                          :: Jdomo
+  integer*8                          :: Jsomo
+  integer*8                          :: diffSOMO
+  integer*8                          :: diffDOMO
+  integer*8                          :: xordiffSOMODOMO
+  integer                            :: ndiffSOMO
+  integer                            :: ndiffDOMO
+  integer                            :: nxordiffSOMODOMO
+  integer                            :: ndiffAll
+  integer                            :: i
+  integer                            :: j
+  integer                            :: k
+  integer                            :: kstart
+  integer                            :: kend
+  integer                            :: Nsomo_I
+  integer                            :: hole
+  integer                            :: p
+  integer                            :: q
+  integer                            :: countalphas
+  logical                            :: pqAlreadyGenQ
+  logical                            :: pqExistsQ
+  logical                            :: ppExistsQ
+
+
+  allocate(tableUniqueAlphas(mo_num,mo_num))
+
+  do idxI = 1, N_configuration
+
+    Icfg  = psi_configuration(:,:,idxI)
+
+    Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
+    Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
+
+    ! find out all pq holes possible
+    nholes = 0
+    ! holes in SOMO
+    do i = 1,mo_num
+       if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
+          nholes += 1
+          listholes(nholes) = i
+          holetype(nholes) = 1
+       endif
+    end do
+    ! holes in DOMO
+    do i = 1,mo_num
+       if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
+          nholes += 1
+          listholes(nholes) = i
+          holetype(nholes) = 2
+       endif
+    end do
+
+    ! find vmos
+    listvmos = -1
+    vmotype = -1
+    nvmos = 0
+    do i = 1,mo_num
+       !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
+       if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then
+          nvmos += 1
+          listvmos(nvmos) = i
+          vmotype(nvmos) = 1
+       else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then
+          nvmos += 1
+          listvmos(nvmos) = i
+          vmotype(nvmos) = 2
+       end if
+    end do
+
+    tableUniqueAlphas = .FALSE.
+
+    ! Now find the allowed (p,q) excitations
+    Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
+    Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
+    Nsomo_I = POPCNT(Isomo)
+    if(Nsomo_I .EQ. 0) then
+      kstart = 1
+    else
+      kstart = cfg_seniority_index(Nsomo_I-2)
+    endif
+    kend = idxI-1
+
+    ! TODO cfg_seniority_index
+    do i = 1,nholes
+       p = listholes(i)
+       do j = 1,nvmos
+          q = listvmos(j)
+          if(p .EQ. q) cycle
+          if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
+             ! SOMO -> VMO
+             Jsomo = IBCLR(Isomo,p-1)
+             Jsomo = IBSET(Jsomo,q-1)
+             Jdomo = Idomo
+             kstart = max(0,cfg_seniority_index(Nsomo_I-2))
+             kend = idxI-1
+          else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
+             ! SOMO -> SOMO
+             Jsomo = IBCLR(Isomo,p-1)
+             Jsomo = IBCLR(Jsomo,q-1)
+             Jdomo = IBSET(Idomo,q-1)
+             kstart = max(0,cfg_seniority_index(Nsomo_I-4))
+             kend = idxI-1
+          else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
+             ! DOMO -> VMO
+             Jsomo = IBSET(Isomo,p-1)
+             Jsomo = IBSET(Jsomo,q-1)
+             Jdomo = IBCLR(Idomo,p-1)
+             kstart = cfg_seniority_index(Nsomo_I)
+             kend = idxI-1
+          else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
+             ! DOMO -> SOMO
+             Jsomo = IBSET(Isomo,p-1)
+             Jsomo = IBCLR(Jsomo,q-1)
+             Jdomo = IBCLR(Idomo,p-1)
+             Jdomo = IBSET(Jdomo,q-1)
+             kstart = max(0,cfg_seniority_index(Nsomo_I-2))
+             kend = idxI-1
+          else
+             print*,"Something went wrong in obtain_associated_alphaI"
+          endif
+
+          ! Again, we don't have to search from 1
+          ! we just use seniorty to find the
+          ! first index with NSOMO - 2 to NSOMO + 2
+          ! this is what is done in kstart, kend
+
+          pqAlreadyGenQ = .FALSE.
+          ! First check if it can be generated before
+          do k = kstart, kend
+             diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)))
+             ndiffSOMO = POPCNT(diffSOMO)
+             if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle
+             diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
+             xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+             ndiffDOMO = POPCNT(diffDOMO)
+             nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+             nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
+             !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
+             if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
+                pqAlreadyGenQ = .TRUE.
+                ppExistsQ = .TRUE.
+                EXIT
+             endif
+             if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
+                pqAlreadyGenQ = .TRUE.
+                EXIT
+             endif
+          end do
+
+          if(pqAlreadyGenQ) cycle
+
+          pqExistsQ = .FALSE.
+
+          if(.NOT. pqExistsQ) then
+             tableUniqueAlphas(p,q) = .TRUE.
+          endif
+       end do
+    end do
+
+  !print *,tableUniqueAlphas(:,:)
+
+  ! prune list of alphas
+  Isomo = Icfg(1,1)
+  Idomo = Icfg(1,2)
+  Jsomo = Icfg(1,1)
+  Jdomo = Icfg(1,2)
+  NalphaIcfg = 0
+  do i = 1, nholes
+     p = listholes(i)
+     do j = 1, nvmos
+        q = listvmos(j)
+        if(p .EQ. q) cycle
+        if(tableUniqueAlphas(p,q)) then
+           if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
+              ! SOMO -> VMO
+              Jsomo = IBCLR(Isomo,p-1)
+              Jsomo = IBSET(Jsomo,q-1)
+              Jdomo = Idomo
+           else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
+              ! SOMO -> SOMO
+              Jsomo = IBCLR(Isomo,p-1)
+              Jsomo = IBCLR(Jsomo,q-1)
+              Jdomo = IBSET(Idomo,q-1)
+           else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
+              ! DOMO -> VMO
+              Jsomo = IBSET(Isomo,p-1)
+              Jsomo = IBSET(Jsomo,q-1)
+              Jdomo = IBCLR(Idomo,p-1)
+           else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
+              ! DOMO -> SOMO
+              Jsomo = IBSET(Isomo,p-1)
+              Jsomo = IBCLR(Jsomo,q-1)
+              Jdomo = IBCLR(Idomo,p-1)
+              Jdomo = IBSET(Jdomo,q-1)
+           else
+              print*,"Something went wrong in obtain_associated_alphaI"
+           endif
+
+           ! SOMO
+           NalphaIcfg += 1
+           !print *,i,j,"|",NalphaIcfg
+           alphasIcfg(1,1,NalphaIcfg) = Jsomo
+           alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+           alphasIcfg_list(1:N_int,1:2,idxI,NalphaIcfg) = alphasIcfg(1:N_int,1:2,NalphaIcfg)
+           NalphaIcfg_list(idxI) = NalphaIcfg
+        endif
+     end do
+  end do
+
+  ! Check if this Icfg has been previously generated as a mono
+  ppExistsQ = .False.
+  Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
+  Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
+  do k = 1, idxI-1
+     diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)))
+     diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
+     xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+     ndiffSOMO = POPCNT(diffSOMO)
+     ndiffDOMO = POPCNT(diffDOMO)
+     nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+     if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
+        ppExistsQ = .TRUE.
+        EXIT
+     endif
+  end do
+  ! Diagonal part (pp,qq)
+  if(nholes > 0 .AND. (.NOT. ppExistsQ))then
+     ! SOMO
+     NalphaIcfg += 1
+     alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1)
+     alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2)
+     alphasIcfg_list(1:N_int,1:2,idxI,NalphaIcfg) = alphasIcfg(1:N_int,1:2,NalphaIcfg)
+     NalphaIcfg_list(idxI) = NalphaIcfg
+  endif
+
+
+  enddo ! end loop idxI
+
+END_PROVIDER
+
   subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg)
   implicit none
   use bitmasks
@@ -155,10 +420,10 @@
         endif
 
         ! Again, we don't have to search from 1
-        ! we just use seniorty to find the 
+        ! we just use seniorty to find the
         ! first index with NSOMO - 2 to NSOMO + 2
         ! this is what is done in kstart, kend
-        
+
         pqAlreadyGenQ = .FALSE.
         ! First check if it can be generated before
         do k = kstart, kend
@@ -259,6 +524,7 @@
            !print *,i,j,"|",NalphaIcfg
            alphasIcfg(1,1,NalphaIcfg) = Jsomo
            alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+           !print *,"I = ",idxI, " Na=",NalphaIcfg," - ",Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
         endif
      end do
   end do

src/csf/sigma_vector.irp.f

@@ -425,6 +425,7 @@ end subroutine get_phase_qp_to_cfg
   END_PROVIDER
 
   BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,NBFMax,NBFMax)]
+  !BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,-1:1,4,NSOMOMax+1,-1:1,NBFMax,NBFMax)]
   use cfunctions
   implicit none
   BEGIN_DOC
@@ -849,7 +850,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
 end subroutine calculate_preconditioner_cfg
 
 
-subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep)
+subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep)
   implicit none
   use bitmasks
   BEGIN_DOC
@@ -910,7 +911,6 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie
   real*8          :: diag_energies(n_CSF)
   !PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
 
-
   !print *," sze = ",sze
   call calculate_preconditioner_cfg(diag_energies)
 
@@ -927,6 +927,11 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie
   !!! Single Excitations !!!
   do i=istart_cfg,iend_cfg
 
+    ! if Seniority_range > 8 then
+    ! continue
+    ! else
+    ! cycle
+
      Icfg(1,1) = psi_configuration(1,1,i)
      Icfg(1,2) = psi_configuration(1,2,i)
      Isomo = Icfg(1,1)
@@ -1041,6 +1046,11 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie
   ! Loop over all selected configurations
   do i = istart_cfg,iend_cfg
 
+     ! if Seniority_range > 8 then
+     ! continue
+     ! else
+     ! cycle
+
      Icfg(1,1) = psi_configuration(1,1,i)
      Icfg(1,2) = psi_configuration(1,2,i)
      starti = psi_config_data(i,1)
@@ -1048,6 +1058,350 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie
 
      ! Returns all unique (checking the past) singly excited cfgs connected to I
      Nalphas_Icfg = 0
+     ! TODO:
+     ! test if size(alphas_Icfg,1) < Nmo**2) then deallocate + allocate
+     !call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
+     Nalphas_Icfg = NalphaIcfg_list(i)
+     alphas_Icfg(1:n_int,1:2,1:Nalphas_Icfg) = alphasIcfg_list(1:n_int,1:2,i,1:Nalphas_Icfg)
+
+     ! TODO : remove doubly excited for return
+     ! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
+     do k = 1,Nalphas_Icfg
+        ! Now generate all singly excited with respect to a given alpha CFG
+        call obtain_connected_I_foralpha(i,alphas_Icfg(1,1,k),connectedI_alpha,idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors)
+
+        if(nconnectedI .EQ. 0) then
+           cycle
+        endif
+        totcolsTKI = 0
+        rowsTKI = -1
+        do j = 1,nconnectedI
+           NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
+           NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
+           p = excitationIds(1,j)
+           q = excitationIds(2,j)
+           extype = excitationTypes(j)
+           call convertOrbIdsToModelSpaceIds(alphas_Icfg(1,1,k), connectedI_alpha(1,1,j), p, q, extype, pmodel, qmodel)
+           ! for E_pp E_rs and E_ppE_rr case
+           if(p.EQ.q) then
+              NSOMOalpha = NSOMOI
+           endif
+           rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
+           colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
+           totcolsTKI += colsikpq
+           if(rowsTKI .LT. rowsikpq .AND. rowsTKI .NE. -1) then
+              print *,">",j,"Something is wrong in sigma-vector", rowsTKI, rowsikpq, "(p,q)=",pmodel,qmodel,"ex=",extype,"na=",NSOMOalpha," nI=",NSOMOI
+              !rowsTKI = rowsikpq
+           else
+              rowsTKI = rowsikpq
+           endif
+        enddo
+
+        allocate(TKI(rowsTKI,n_st,totcolsTKI)) ! coefficients of CSF
+        ! Initialize the inegral container
+        ! dims : (totcolsTKI, nconnectedI)
+        allocate(GIJpqrs(totcolsTKI,nconnectedI))  ! gpqrs
+        allocate(TKIGIJ(rowsTKI,n_st,nconnectedI))  ! gpqrs
+
+        totcolsTKI = 0
+        do j = 1,nconnectedI
+           NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
+           NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
+           p = excitationIds(1,j)
+           q = excitationIds(2,j)
+           extype = excitationTypes(j)
+           call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
+           rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
+           colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
+           do kk = 1,n_st
+           do l = 1,rowsTKI
+              do m = 1,colsikpq
+                 TKI(l,kk,totcolsTKI+m) = AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * psi_in(idxs_connectedI_alpha(j)+m-1,kk)
+              enddo
+           enddo
+           enddo
+           do m = 1,colsikpq
+              do l = 1,nconnectedI
+                 ! <ij|kl> = (ik|jl)
+                 moi = excitationIds(1,j) ! p
+                 mok = excitationIds(2,j) ! q
+                 moj = excitationIds(2,l) ! s
+                 mol = excitationIds(1,l) ! r
+                 if(moi.EQ.mok .AND. moj.EQ.mol)then
+                    diagfac = diagfactors(j)
+                    diagfac *= diagfactors(l)
+                    !print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
+                    GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
+                 else
+                       diagfac = diagfactors(j)*diagfactors(l)
+                       !print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
+                       GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
+                    !endif
+                 endif
+              enddo
+           enddo
+           totcolsTKI += colsikpq
+        enddo
+
+
+
+        ! Do big BLAS
+        ! TODO TKI, size(TKI,1)*size(TKI,2)
+        call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0,  &
+             TKI, size(TKI,1)*n_st, GIJpqrs, size(GIJpqrs,1), 0.d0, &
+             TKIGIJ , size(TKIGIJ,1)*n_st )
+
+
+        ! Collect the result
+        totcolsTKI = 0
+        do j = 1,nconnectedI
+           NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
+           NSOMOI     = getNSOMO(connectedI_alpha(:,:,j))
+           p = excitationIds(1,j)
+           q = excitationIds(2,j)
+           extype = excitationTypes(j)
+           call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
+           rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
+           colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
+           !print *,">j=",j,rowsikpq,colsikpq, ">>",totcolsTKI,",",idxs_connectedI_alpha(j)
+           do kk = 1,n_st
+           do m = 1,colsikpq
+              do l = 1,rowsTKI
+                 psi_out(idxs_connectedI_alpha(j)+m-1,kk) += AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * TKIGIJ(l,kk,j)
+                 psi_coef_out_init(idxs_connectedI_alpha(j)+m-1) = .True.
+              enddo
+           enddo
+           enddo
+           totcolsTKI += colsikpq
+        enddo
+
+        deallocate(TKI) ! coefficients of CSF
+        ! Initialize the inegral container
+        ! dims : (totcolsTKI, nconnectedI)
+        deallocate(GIJpqrs)  ! gpqrs
+        deallocate(TKIGIJ)  ! gpqrs
+
+     enddo ! loop over alphas
+  enddo ! loop over I
+
+
+  ! Add the diagonal contribution
+  do i = 1,n_CSF
+     psi_out(i,1) += 1.0d0*diag_energies(i)*psi_in(i,1)
+  enddo
+
+
+end subroutine calculate_sigma_vector_cfg_nst
+subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep)
+  implicit none
+  use bitmasks
+  BEGIN_DOC
+  ! Documentation for sigma-vector calculation
+  !
+  ! Calculates the result of the
+  ! application of the hamiltonian to the
+  ! wavefunction in CFG basis once
+  ! TODO : Things prepare outside this routine
+  !  1. Touch the providers for
+  !     a. ApqIJ containers
+  !     b. DET to CSF transformation matrices
+  !  2. DET to CSF transcormation
+  !  2. CSF to DET back transcormation
+  ! returns : psi_coef_out_det :
+  END_DOC
+  integer,intent(in) :: sze, istart,iend, istep, ishift, n_st
+  real*8,intent(in):: psi_in(sze,n_st)
+  real*8,intent(out):: psi_out(sze,n_st)
+  integer(bit_kind) :: Icfg(N_INT,2)
+  integer :: i,j,k,l,p,q,noccp,noccq, ii, jj, m, n, idxI, kk, nocck,orbk
+  integer(bit_kind) :: alphas_Icfg(N_INT,2,max(sze,100))
+  integer(bit_kind) :: singlesI(N_INT,2,max(sze,100))
+  integer(bit_kind) :: connectedI_alpha(N_INT,2,max(sze,100))
+  integer           :: idxs_singlesI(max(sze,100))
+  integer           :: idxs_connectedI_alpha(max(sze,100))
+  integer(bit_kind) :: psi_configuration_out(N_INT,2,max(sze,100))
+  real*8            :: psi_coef_out(n_CSF)
+  logical           :: psi_coef_out_init(n_CSF)
+  integer           :: excitationIds_single(2,max(sze,100))
+  integer           :: excitationTypes_single(max(sze,100))
+  integer           :: excitationIds(2,max(sze,100))
+  integer           :: excitationTypes(max(sze,100))
+  real*8            :: diagfactors(max(sze,100))
+  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  :: Nalphas_Icfg, nconnectedI, rowsikpq, colsikpq, nsinglesI
+  integer  :: extype,NSOMOalpha,NSOMOI,NSOMOJ,pmodel,qmodel
+  integer :: getNSOMO
+  integer :: totcolsTKI
+  integer :: rowsTKI
+  integer :: noccpp
+  integer :: istart_cfg, iend_cfg
+  integer*8 :: MS, Isomo, Idomo, Jsomo, Jdomo, Ialpha, Ibeta
+  integer :: moi, moj, mok, mol, starti, endi, startj, endj, cnti, cntj, cntk
+  real*8  :: norm_coef_cfg, fac2eints
+  real*8  :: norm_coef_det
+  real*8  :: meCC1, meCC2, diagfac
+  real*8,dimension(:,:,:),allocatable :: TKI
+  real*8,dimension(:,:),allocatable  :: GIJpqrs
+  real*8,dimension(:,:,:),allocatable  :: TKIGIJ
+  real*8, external :: mo_two_e_integral
+  real*8, external :: get_two_e_integral
+  real*8          :: diag_energies(n_CSF)
+  !PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
+
+  !print *," sze = ",sze
+  call calculate_preconditioner_cfg(diag_energies)
+
+  MS = 0
+  norm_coef_cfg=0.d0
+
+  psi_out=0.d0
+  psi_coef_out_init = .False.
+
+  istart_cfg = psi_csf_to_config_data(istart)
+  iend_cfg   = psi_csf_to_config_data(iend)
+
+
+  !!! Single Excitations !!!
+  do i=istart_cfg,iend_cfg
+
+    ! if Seniority_range > 8 then
+    ! continue
+    ! else
+    ! cycle
+
+     Icfg(1,1) = psi_configuration(1,1,i)
+     Icfg(1,2) = psi_configuration(1,2,i)
+     Isomo = Icfg(1,1)
+     Idomo = Icfg(1,2)
+     NSOMOI = getNSOMO(Icfg)
+
+     ! find out all pq holes possible
+     nholes = 0
+     ! holes in SOMO
+     ! list_act
+     ! list_core
+     ! list_core_inact
+     ! bitmasks
+     !do k = n_core_orb+1,n_core_orb + n_act_orb
+     do k = 1,mo_num
+        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
+     do k = 1,mo_num
+        if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
+           nholes += 1
+           listholes(nholes) = k
+           holetype(nholes) = 2
+        endif
+     enddo
+
+     ! find vmos
+     listvmos = -1
+     vmotype = -1
+     nvmos = 0
+     !do k = n_core_orb+1,n_core_orb + n_act_orb
+     do k = 1,mo_num
+        !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
+        if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0) then
+           nvmos += 1
+           listvmos(nvmos) = k
+           vmotype(nvmos) = 0
+        else if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0 ) then
+           nvmos += 1
+           listvmos(nvmos) = k
+           vmotype(nvmos) = 1
+        end if
+     enddo
+
+
+     ! Icsf ids
+     starti = psi_config_data(i,1)
+     endi   = psi_config_data(i,2)
+     NSOMOI = getNSOMO(Icfg)
+
+     call generate_all_singles_cfg_with_type(Icfg,singlesI,idxs_singlesI,excitationIds_single, &
+          excitationTypes_single,nsinglesI,N_int)
+
+     do j = 1,nsinglesI
+        idxI = idxs_singlesI(j)
+        NSOMOJ = getNSOMO(singlesI(:,:,j))
+        p = excitationIds_single(1,j)
+        q = excitationIds_single(2,j)
+        extype = excitationTypes_single(j)
+        ! Off diagonal terms
+        call convertOrbIdsToModelSpaceIds(Icfg, singlesI(:,:,j), p, q, extype, pmodel, qmodel)
+        Jsomo = singlesI(1,1,j)
+        Jdomo = singlesI(1,2,j)
+
+        ! Add the hole on J
+        if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1  .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
+           nholes += 1
+           listholes(nholes) = q
+           holetype(nholes) = 1
+        endif
+        if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
+           nholes += 1
+           listholes(nholes) = q
+           holetype(nholes) = 2
+        endif
+
+        startj = psi_config_data(idxI,1)
+        endj   = psi_config_data(idxI,2)
+
+        !!! One-electron contribution !!!
+        do kk = 1,n_st
+        cnti = 0
+        do ii = starti, endi
+           cnti += 1
+           cntj = 0
+           do jj = startj, endj
+              cntj += 1
+              meCC1 = AIJpqContainer(NSOMOI,NSOMOJ,extype,pmodel,qmodel,cnti,cntj)
+              psi_out(jj,kk) += meCC1 * psi_in(ii,kk) * h_core_ri(p,q)
+              psi_coef_out_init(jj) = .True.
+           enddo
+        enddo
+        enddo
+
+        ! Undo setting in listholes
+        if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1  .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
+           nholes -= 1
+        endif
+        if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
+           nholes -= 1
+        endif
+     enddo
+  enddo
+
+  !!! Double Excitations !!!
+
+  ! Loop over all selected configurations
+  do i = istart_cfg,iend_cfg
+
+     ! if Seniority_range > 8 then
+     ! continue
+     ! else
+     ! cycle
+
+     Icfg(1,1) = psi_configuration(1,1,i)
+     Icfg(1,2) = psi_configuration(1,2,i)
+     starti = psi_config_data(i,1)
+     endi   = psi_config_data(i,2)
+
+     ! Returns all unique (checking the past) singly excited cfgs connected to I
+     Nalphas_Icfg = 0
+     ! TODO:
+     ! test if size(alphas_Icfg,1) < Nmo**2) then deallocate + allocate
      call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
      ! TODO : remove doubly excited for return
      ! Here we do 2x the loop. One to count for the size of the matrix, then we compute.

src/davidson/diagonalization_hcsf_dressed.irp.f

@@ -306,22 +306,22 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
         PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
         if ((sze > 100000).and.distributed_davidson) then
 
-!             call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
-!             call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W)
-!             call H_u_0_nstates_zmq   (W,U,N_st_diag,sze)
-!             call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1))
-!             call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
+!            call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
+!            call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W)
+!            call H_u_0_nstates_zmq   (W,U,N_st_diag,sze)
+!            call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1))
+!            call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
 !             call calculate_sigma_vector_cfg_nst(W_csf(1,shift+1),U_csf(1,shift+1),N_st_diag,sze_csf,1,sze_csf,0,1)
 !            ! TODO : psi_det_size ? for psi_det
             do kk=1,N_st_diag
                 call calculate_sigma_vector_cfg_nst(W_csf(1,shift+kk),U_csf(1,shift+kk),1,sze_csf,1,sze_csf,0,1)
             enddo
         else
-!             call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
-!             call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W)
-!             call H_u_0_nstates_openmp(W,U,N_st_diag,sze)
-!             call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1))
-!             call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
+            !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
+            !call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W)
+            !call H_u_0_nstates_openmp(W,U,N_st_diag,sze)
+            !call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1))
+            !call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
 !             call calculate_sigma_vector_cfg_nst(W_csf(1,shift+1),U_csf(1,shift+1),N_st_diag,sze_csf,1,sze_csf,0,1)
             do kk=1,N_st_diag
                 call calculate_sigma_vector_cfg_nst(W_csf(1,shift+kk),U_csf(1,shift+kk),1,sze_csf,1,sze_csf,0,1)