Working on fixing n_int.

src/csf/configuration_CI_sigma_helpers.irp.f

@@ -12,6 +12,7 @@ use bitmasks
 
   integer                        :: idxI ! The id of the Ith CFG
   integer(bit_kind)              :: Icfg(N_int,2)
+  integer(bit_kind)              :: Jcfg(N_int,2)
   integer                        :: NalphaIcfg
   logical,dimension(:,:),allocatable :: tableUniqueAlphas
   integer                        :: listholes(mo_num)
@@ -20,10 +21,10 @@ use bitmasks
   integer                        :: nvmos
   integer                        :: listvmos(mo_num)
   integer                        :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
-  integer*8                      :: Idomo
+  integer*8                      :: Idomo, Idomop, Idomoq 
-  integer*8                      :: Isomo
+  integer*8                      :: Isomo, Isomop, Isomoq
-  integer*8                      :: Jdomo
+  integer*8                      :: Jdomo, Jdomop, Jdomoq
-  integer*8                      :: Jsomo
+  integer*8                      :: Jsomo, Jsomop, Jsomoq
   integer*8                      :: diffSOMO
   integer*8                      :: diffDOMO
   integer*8                      :: xordiffSOMODOMO
@@ -31,20 +32,21 @@ use bitmasks
   integer                        :: ndiffDOMO
   integer                        :: nxordiffSOMODOMO
   integer                        :: ndiffAll
-  integer                        :: i,ii
+  integer                        :: i,ii,iii
-  integer                        :: j,jj
+  integer                        :: j,jj, i_s, i_d
   integer                        :: k,kk
   integer                        :: kstart
   integer                        :: kend
-  integer                        :: Nsomo_I
+  integer                        :: Nsomo_I, Nsomo_J
-  integer                        :: hole
+  integer                        :: hole, n_core_orb_64
-  integer                        :: p
+  integer                        :: p, pp, p_s
-  integer                        :: q
+  integer                        :: q, qq, q_s
   integer                        :: countalphas
   logical                        :: pqAlreadyGenQ
   logical                        :: pqExistsQ
   logical                        :: ppExistsQ
   integer*8                      :: MS
+  integer :: listall(N_int*bit_kind_size), nelall
 
   double precision               :: t0, t1
   call wall_time(t0)
@@ -57,6 +59,9 @@ use bitmasks
   do idxI = 1, N_configuration
 
     Icfg  = psi_configuration(:,:,idxI)
+    Jcfg  = psi_configuration(:,:,idxI)
+    !print *," Jcfg somo=",Jcfg(1,1), " ", Jcfg(2,1)
+    !print *," Jcfg domo=",Jcfg(1,2), " ", Jcfg(2,2)
 
     Isomo = iand(act_bitmask(1,1),Icfg(1,1))
     Idomo = iand(act_bitmask(1,1),Icfg(1,2))
@@ -64,38 +69,89 @@ use bitmasks
     ! find out all pq holes possible
     nholes = 0
     ! holes in SOMO
-    do ii = 1,n_act_orb
+    !do ii = 1,n_act_orb
-      i = list_act(ii)
+    !  i = list_act(ii)
-      if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
+    !  if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
-        nholes += 1
+    !    nholes += 1
-        listholes(nholes) = i
+    !    listholes(nholes) = i
-        holetype(nholes) = 1
+    !    holetype(nholes) = 1
-      endif
+    !  endif
-    end do
+    !end do
+        call bitstring_to_list(psi_configuration(1,1,idxI),listall,nelall,N_int)
+
+        !print *,'list somo'
+        do iii=1,nelall
+          nholes += 1
+          listholes(nholes) = listall(iii)
+            !print *,listall(iii)
+          holetype(nholes) = 1
+        end do
+
+        Nsomo_I = nelall
+
     ! holes in DOMO
-    do ii = 1,n_act_orb
+    !do ii = 1,n_act_orb
-      i = list_act(ii)
+    !  i = list_act(ii)
-      if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
+    !  if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
-        nholes += 1
+    !    nholes += 1
-        listholes(nholes) = i
+    !    listholes(nholes) = i
-        holetype(nholes) = 2
+    !    holetype(nholes) = 2
-      endif
+    !  endif
-    end do
+    !end do
+
+    !do iii=1,N_int
+    !  print *,' iii=',iii, psi_configuration(iii,2,idxI), ' idxI=',idxI
+    !end do
+        call bitstring_to_list(psi_configuration(1,2,idxI),listall,nelall,N_int)
+
+        !print *,'list domo ncore=',n_core_orb, ' nelall=',nelall
+        do iii=1,nelall
+          if(listall(iii) .gt. n_core_orb)then
+            nholes += 1
+            listholes(nholes) = listall(iii)
+            !print *,listall(iii)
+            holetype(nholes) = 2
+          endif
+        end do
 
     ! find vmos
     listvmos = -1
     vmotype = -1
     nvmos = 0
-    do ii = 1,n_act_orb
+    !do ii = 1,n_act_orb
-      i = list_act(ii)
+    !  i = list_act(ii)
+    !  if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then
+    !    if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then
+    !      nvmos += 1
+    !      listvmos(nvmos) = i
+    !      print *,'1 i=',i
+    !      vmotype(nvmos) = 1
+    !    else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then
+    !      nvmos += 1
+    !      listvmos(nvmos) = i
+    !      print *,'2 i=',i
+    !      vmotype(nvmos) = 2
+    !    end if
+    !  end if
+    !end do
+    !print *,'-----------'
+
+    ! Take into account N_int
+    do ii = 1, n_act_orb
+      iii = list_act(ii)
+      i_s = (1+((iii-1)/63))
+      i = iii - ( i_s -1 )*63 
+      Isomo = iand(act_bitmask(i_s,1),Icfg(i_s,1))
+      Idomo = iand(act_bitmask(i_s,1),Icfg(i_s,2))
+
       if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then
         if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then
           nvmos += 1
-          listvmos(nvmos) = i
+          listvmos(nvmos) = iii
           vmotype(nvmos) = 1
         else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then
           nvmos += 1
-          listvmos(nvmos) = i
+          listvmos(nvmos) = iii
           vmotype(nvmos) = 2
         end if
       end if
@@ -106,7 +162,7 @@ use bitmasks
     ! Now find the allowed (p,q) excitations
     Isomo = iand(act_bitmask(1,1),Icfg(1,1))
     Idomo = iand(act_bitmask(1,1),Icfg(1,2))
-    Nsomo_I = POPCNT(Isomo)
+    !Nsomo_I = POPCNT(Isomo)
     if(Nsomo_I .EQ. 0) then
       kstart = 1
     else
@@ -115,24 +171,54 @@ use bitmasks
     kend = idxI-1
 
     do i = 1,nholes
-      p = listholes(i)
+      pp  = listholes(i)
+      p_s = (1+((pp-1)/63))
+      p   = pp - (p_s - 1)*63
+      !print *,' pp=',pp, ' p_s=',p_s, ' p=',p
       do j = 1,nvmos
-        q = listvmos(j)
+        qq  = listvmos(j)
+        q_s = (1+((qq-1)/63))
+        q   = qq - (q_s - 1)*63
+        !print *,' qq=',qq, ' q_s=',q_s, ' q=',q
+        Isomop = iand(act_bitmask(i_s,1),Icfg(p_s,1))
+        Idomop = iand(act_bitmask(i_s,1),Icfg(p_s,2))
+        Isomop = iand(act_bitmask(i_s,1),Icfg(q_s,1))
+        Idomop = iand(act_bitmask(i_s,1),Icfg(q_s,2))
         if(p .EQ. q) cycle
         if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
           ! SOMO -> VMO
-          Jsomo = IBCLR(Isomo,p-1)
+          !print *,'SOMO -> VMO'
-          Jsomo = IBSET(Jsomo,q-1)
+          if (p_s .eq. q_s) then
-          Jdomo = Idomo
+            Jsomop = IBCLR(Isomop,p-1)
+            Jsomop = IBSET(Jsomop,q-1)
+            Jsomoq = Jsomop
+          else
+            Jsomop = IBCLR(Isomop,p-1)
+            Jsomoq = IBSET(Isomoq,q-1)
+          endif
+
+          ! Domo remains the same
+          Jdomop = Idomop
+          Jdomoq = Idomoq
+
           kstart = max(1,cfg_seniority_index(max(NSOMOMin,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)
+          !print *,'SOMO -> SOMO'
-          Jsomo = IBCLR(Jsomo,q-1)
+          if(p_s .eq. q_s) then
-          Jdomo = IBSET(Idomo,q-1)
+            Jsomop = IBCLR(Isomop,p-1)
+            Jsomop = IBCLR(Jsomop,q-1)
+            Jsomoq = Jsomop
+          else
+            Jsomop = IBCLR(Isomop,p-1)
+            Jsomoq = IBCLR(Isomoq,q-1)
+          endif
+
+          Jdomoq = IBSET(Idomoq,q-1)
+
           ! Check for Minimal alpha electrons (MS)
-          if(POPCNT(Jsomo).ge.MS)then
+          if(POPCNT(Jsomoq).ge.MS)then
             kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
             kend = idxI-1
           else
@@ -140,24 +226,60 @@ use bitmasks
           endif
         else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
           ! DOMO -> VMO
-          Jsomo = IBSET(Isomo,p-1)
+          !print *,'DOMO -> VMO', Isomop, p, q, Jsomop
-          Jsomo = IBSET(Jsomo,q-1)
+          if(p_s .eq. q_s) then
-          Jdomo = IBCLR(Idomo,p-1)
+            Jsomop = IBSET(Isomop,p-1)
+            Jsomop = IBSET(Jsomop,q-1)
+            Jsomoq = Jsomop
+          else
+            Jsomop = IBSET(Isomop,p-1)
+            Jsomoq = IBSET(Jsomoq,q-1)
+          endif
+          !print *, 'Jsomop=', Jsomop
+
+          Jdomop = IBCLR(Idomop,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)
+          !print *,'DOMO -> SOMO'
-          Jsomo = IBCLR(Jsomo,q-1)
+          if(p_s .eq. q_s) then
-          Jdomo = IBCLR(Idomo,p-1)
+            Jsomop = IBSET(Isomop,p-1)
-          Jdomo = IBSET(Jdomo,q-1)
+            Jsomop = IBCLR(Jsomop,q-1)
+            Jsomoq = Jsomop
+
+            Jdomop = IBCLR(Idomop,p-1)
+            Jdomop = IBSET(Jdomop,q-1)
+            Jdomoq = Jdomop
+          else
+            Jsomop = IBSET(Isomop,p-1)
+            Jsomoq = IBCLR(Jsomoq,q-1)
+
+            Jdomop = IBCLR(Idomop,p-1)
+            Jdomoq = IBSET(Jdomoq,q-1)
+          endif
+
           kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
           kend = idxI-1
         else
           print*,"Something went wrong in obtain_associated_alphaI"
         endif
+
+        ! Save it to Jcfg
+        !print *,i,j,"0| nalpha=",NalphaIcfg, " somo=",Jcfg(1,1),Jcfg(2,1)
+        Jcfg(p_s,1) = Jsomop
+        Jcfg(q_s,1) = Jsomoq
+        Jcfg(p_s,2) = Jdomop
+        Jcfg(q_s,2) = Jdomoq
+        !print *,'p_s=',p_s,' q_s=', q_s
+        !print *,'Jsomop=',Jsomop, ' Jsomoq=', Jsomoq, ' Jdomop=', Jdomop, ' Jdomoq=', Jdomo
+        !print *,i,j,"1| nalpha=",NalphaIcfg, " somo=",Jcfg(1,1),Jcfg(2,1)
+        call bitstring_to_list(Jcfg(1,1),listall,nelall,N_int)
+        Nsomo_J = nelall
+
         ! Check for Minimal alpha electrons (MS)
-        if(POPCNT(Jsomo).lt.MS)then
+        if(Nsomo_J.lt.MS)then
           cycle
         endif
 
@@ -169,15 +291,32 @@ use bitmasks
         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)))
+          !diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)))
-          ndiffSOMO = POPCNT(diffSOMO)
+          !ndiffSOMO = POPCNT(diffSOMO)
-          if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle
+          !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)))
+          !diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
-          xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+          !xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
-          ndiffDOMO = POPCNT(diffDOMO)
+          !ndiffDOMO = POPCNT(diffDOMO)
-          nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+          !nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
-          nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
+          !nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
-          !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
+
+          ndiffSOMO = 0
+          ndiffDOMO = 0
+          nxordiffSOMODOMO = 0
+          do ii = 1, N_int
+            Jsomo = Jcfg(ii,1)
+            Jdomo = Jcfg(ii,2)
+            diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k)))
+            ndiffSOMO += POPCNT(diffSOMO)
+            diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k)))
+            xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+            ndiffDOMO += POPCNT(diffDOMO)
+            nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
+            nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
+          end do
+
+          if((ndiffSOMO .ne. 0) .and. (ndiffSOMO .ne. 2)) cycle
+
           if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
             pqAlreadyGenQ = .TRUE.
             ppExistsQ = .TRUE.
@@ -208,22 +347,57 @@ use bitmasks
     Jdomo = Icfg(1,2)
     NalphaIcfg = 0
     do i = 1, nholes
-      p = listholes(i)
+      !p = listholes(i)
+      pp  = listholes(i)
+      p_s = (1+((pp-1)/63))
+      p   = pp - (p_s - 1)*63
       do j = 1, nvmos
-        q = listvmos(j)
+        !q = listvmos(j)
+        qq  = listvmos(j)
+        q_s = (1+((qq-1)/63))
+        q   = qq - (q_s - 1)*63
+        Isomop = iand(act_bitmask(i_s,1),Icfg(p_s,1))
+        Idomop = iand(act_bitmask(i_s,1),Icfg(p_s,2))
+        Isomoq = iand(act_bitmask(i_s,1),Icfg(q_s,1))
+        Idomoq = iand(act_bitmask(i_s,1),Icfg(q_s,2))
         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 = IBCLR(Isomo,p-1)
-            Jsomo = IBSET(Jsomo,q-1)
+            !Jsomo = IBSET(Jsomo,q-1)
-            Jdomo = Idomo
+            !Jdomo = Idomo
+            if (p_s .eq. q_s) then
+              Jsomop = IBCLR(Isomop,p-1)
+              Jsomop = IBSET(Jsomop,q-1)
+              Jsomoq = Jsomop
+            else
+              Jsomop = IBCLR(Isomop,p-1)
+              Jsomoq = IBSET(Isomoq,q-1)
+            endif
+
+            ! Domo remains the same
+            Jdomop = Idomop
+            Jdomoq = Idomoq
+
           else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
             ! SOMO -> SOMO
-            Jsomo = IBCLR(Isomo,p-1)
+            !Jsomo = IBCLR(Isomo,p-1)
-            Jsomo = IBCLR(Jsomo,q-1)
+            !Jsomo = IBCLR(Jsomo,q-1)
-            Jdomo = IBSET(Idomo,q-1)
+            !Jdomo = IBSET(Idomo,q-1)
-            if(POPCNT(Jsomo).ge.MS)then
+
+            if(p_s .eq. q_s) then
+              Jsomop = IBCLR(Isomop,p-1)
+              Jsomop = IBCLR(Jsomop,q-1)
+              Jsomoq = Jsomop
+            else
+              Jsomop = IBCLR(Isomop,p-1)
+              Jsomoq = IBCLR(Isomoq,q-1)
+            endif
+
+            Jdomoq = IBSET(Idomoq,q-1)
+
+            if(POPCNT(Jsomoq).ge.MS)then
               kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
               kend = idxI-1
             else
@@ -231,26 +405,74 @@ use bitmasks
             endif
           else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
             ! DOMO -> VMO
-            Jsomo = IBSET(Isomo,p-1)
+            !Jsomo = IBSET(Isomo,p-1)
-            Jsomo = IBSET(Jsomo,q-1)
+            !Jsomo = IBSET(Jsomo,q-1)
-            Jdomo = IBCLR(Idomo,p-1)
+            !Jdomo = IBCLR(Idomo,p-1)
+
+            if(p_s .eq. q_s) then
+              Jsomop = IBSET(Isomop,p-1)
+              Jsomop = IBSET(Jsomop,q-1)
+              Jsomoq = Jsomop
+            else
+              Jsomop = IBSET(Isomop,p-1)
+              Jsomoq = IBSET(Jsomoq,q-1)
+            endif
+
+            Jdomop = IBCLR(Idomop,p-1)
+
           else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
             ! DOMO -> SOMO
-            Jsomo = IBSET(Isomo,p-1)
+            !Jsomo = IBSET(Isomo,p-1)
-            Jsomo = IBCLR(Jsomo,q-1)
+            !Jsomo = IBCLR(Jsomo,q-1)
-            Jdomo = IBCLR(Idomo,p-1)
+            !Jdomo = IBCLR(Idomo,p-1)
-            Jdomo = IBSET(Jdomo,q-1)
+            !Jdomo = IBSET(Jdomo,q-1)
+            if(p_s .eq. q_s) then
+              Jsomop = IBSET(Isomop,p-1)
+              Jsomop = IBCLR(Jsomop,q-1)
+              Jsomoq = Jsomop
+
+              Jdomop = IBCLR(Idomop,p-1)
+              Jdomop = IBSET(Jdomop,q-1)
+              Jdomoq = Jdomop
+            else
+              Jsomop = IBSET(Isomop,p-1)
+              Jsomoq = IBCLR(Jsomoq,q-1)
+
+              Jdomop = IBCLR(Idomop,p-1)
+              Jdomoq = IBSET(Jdomoq,q-1)
+            endif
+
           else
             print*,"Something went wrong in obtain_associated_alphaI"
           endif
 
+          ! Save it to Jcfg
+          Jcfg(p_s,1) = Jsomop
+          Jcfg(q_s,1) = Jsomoq
+          Jcfg(p_s,2) = Jdomop
+          Jcfg(q_s,2) = Jdomoq
+
           ! SOMO
           !print *,i,j,"|",NalphaIcfg, Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
           if(POPCNT(Jsomo) .ge. NSOMOMin) then
             NalphaIcfg += 1
-            alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo
+            alphasIcfg_list(:,1,idxI,NalphaIcfg) = Jcfg(:,1)
-            alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+            !alphasIcfg_list(:,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+            if(n_core_orb .le. 63)then
+              alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jcfg(1,2),ISHFT(1_8,n_core_orb)-1)
+            else
+              n_core_orb_64 = n_core_orb
+              do ii=1,N_int
+                if(n_core_orb_64 .gt. 0)then
+                  alphasIcfg_list(ii,2,idxI,NalphaIcfg) = IOR(Jcfg(ii,2),ISHFT(1_8,n_core_orb_64)-1)
+                else
+                  alphasIcfg_list(ii,2,idxI,NalphaIcfg) = Jcfg(ii,2)
+                endif
+                n_core_orb_64 = ISHFT(n_core_orb_64,-6)
+              end do
+            endif
             NalphaIcfg_list(idxI) = NalphaIcfg
+            !print *,i,j,"2| nalpha=",NalphaIcfg, " somo=",Jcfg(1,1),Jcfg(2,1)
           endif
         endif
       end do
@@ -261,14 +483,24 @@ use bitmasks
     Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
     Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
     kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
+    ndiffDOMO = 0
     do k = kstart, idxI-1
-      diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k)))
+      do ii=1,N_int
-      ndiffSOMO = POPCNT(diffSOMO)
+        diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k)))
+        ndiffSOMO += POPCNT(diffSOMO)
+      end do
+      ! ndiffSOMO cannot be 0 (I /= k)
+      ! if ndiffSOMO /= 2 then it has to be greater than 2 and hense
+      ! this Icfg could not have been generated before.
       if (ndiffSOMO /= 2) cycle
-      diffDOMO = IEOR(Idomo,iand(act_bitmask(1,1),psi_configuration(1,2,k)))
+      ndiffDOMO = 0
-      xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+      nxordiffSOMODOMO = 0
-      ndiffDOMO = POPCNT(diffDOMO)
+      do ii=1,N_int
-      nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+        diffDOMO = IEOR(Icfg(ii,2),iand(act_bitmask(ii,1),psi_configuration(ii,2,k)))
+        xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+        ndiffDOMO += POPCNT(diffDOMO)
+        nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
+      end do
       if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then
         ppExistsQ = .TRUE.
         EXIT
@@ -279,8 +511,8 @@ use bitmasks
       ! SOMO
       if(POPCNT(Jsomo) .ge. NSOMOMin) then
         NalphaIcfg += 1
-        alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1)
+        alphasIcfg_list(:,1,idxI,NalphaIcfg) = Icfg(:,1)
-        alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2)
+        alphasIcfg_list(:,2,idxI,NalphaIcfg) = Icfg(:,2)
         NalphaIcfg_list(idxI) = NalphaIcfg
       endif
     endif

src/csf/configurations.irp.f

@@ -352,6 +352,11 @@ end
    psi_configuration(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i))
    psi_configuration(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i))
   enddo
+ if(i.eq.1)then
+   print *,'Preparing PSI_CONFIGURATION i=',i
+            print *," Icfg somo=",psi_configuration(1,1,1), " ", psi_configuration(2,1,1)
+            print *," Icfg domo=",psi_configuration(1,2,1), " ", psi_configuration(2,2,1)
+ endif
  enddo
 
  ! Sort

src/csf/obtain_I_foralpha.irp.f

@@ -38,6 +38,7 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
   integer :: holetype(mo_num)
   integer :: end_index
   integer :: Nsomo_I
+  integer :: listall(N_int*bit_kind_size), nelall
 
   ! 
   ! 2 2 1 1 0 0 : 1 1 0 0 0 0
@@ -65,9 +66,12 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
 
   ! Since CFGs are sorted wrt to seniority
   ! we don't have to search the full CFG list
-  Isomo = givenI(1,1)
+  Nsomo_I = 0
-  Idomo = givenI(1,2)
+  do i=1,N_int
-  Nsomo_I = POPCNT(Isomo)
+    Isomo = givenI(i,1)
+    Idomo = givenI(i,2)
+    Nsomo_I += POPCNT(Isomo)
+  end do
   end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_I+6,elec_num))-1)
   if(end_index .LT. 0) end_index= N_configuration
   !end_index = N_configuration
@@ -83,17 +87,24 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
      !  idxs_connectedI(nconnectedI)=i
      !  cycle
      !endif
-     Isomo = givenI(1,1)
+
-     Idomo = givenI(1,2)
+     ndiffSOMO = 0
-     Jsomo = psi_configuration(1,1,i)
+     ndiffDOMO = 0
-     Jdomo = psi_configuration(1,2,i)
+     nxordiffSOMODOMO = 0
-     diffSOMO = IEOR(Isomo,Jsomo)
+     do ii=1,N_int
-     ndiffSOMO = POPCNT(diffSOMO)
+       Isomo = givenI(ii,1)
-     diffDOMO = IEOR(Idomo,Jdomo)
+       Idomo = givenI(ii,2)
-     xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+       Jsomo = psi_configuration(ii,1,i)
-     ndiffDOMO = POPCNT(diffDOMO)
+       Jdomo = psi_configuration(ii,2,i)
-     nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+       diffSOMO = IEOR(Isomo,Jsomo)
-     nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO 
+       ndiffSOMO += POPCNT(diffSOMO)
+       diffDOMO = IEOR(Idomo,Jdomo)
+       xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+       ndiffDOMO += POPCNT(diffDOMO)
+       nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
+       nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO 
+     end do
+
      if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
        !-------
        ! MONO |
@@ -144,25 +155,45 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
         ! find out all pq holes possible
         nholes = 0
         ! holes in SOMO
-        Isomo = psi_configuration(1,1,i)
+        !Isomo = psi_configuration(1,1,i)
-        Idomo = psi_configuration(1,2,i)
+        !Idomo = psi_configuration(1,2,i)
-        do iii = 1,n_act_orb
+        !do iii = 1,n_act_orb
-          ii = list_act(iii)
+        !  ii = list_act(iii)
-           if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
+        !   if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
-              nholes += 1
+        !      nholes += 1
-              listholes(nholes) = ii
+        !      listholes(nholes) = ii
-              holetype(nholes) = 1
+        !      holetype(nholes) = 1
-           endif
+        !   endif
+        !end do
+       
+        call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int)
+
+        do iii=1,nelall
+          nholes += 1
+          listholes(nholes) = listall(iii)
+          holetype(nholes) = 1
         end do
+
         ! holes in DOMO
-        do iii = 1,n_act_orb
+        !do iii = 1,n_act_orb
-          ii = list_act(iii)
+        !  ii = list_act(iii)
-           if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
+        !   if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
-              nholes += 1
+        !      nholes += 1
-              listholes(nholes) = ii
+        !      listholes(nholes) = ii
-              holetype(nholes) = 2
+        !      holetype(nholes) = 2
-           endif
+        !   endif
+        !end do
+        
+        call bitstring_to_list(psi_configuration(1,2,i),listall,nelall,N_int)
+
+        do iii=1,nelall
+          if(listall(iii) .gt. n_core_orb)then
+            nholes += 1
+            listholes(nholes) = listall(iii)
+            holetype(nholes) = 2
+          endif
         end do
+
         ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1)*nholes)
      endif
   end do
@@ -199,6 +230,8 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
   integer*8                                :: Isomo
   integer*8                                :: Jdomo
   integer*8                                :: Jsomo
+  integer(bit_kind)                        :: Jcfg(N_int,2)
+  integer(bit_kind)                        :: Icfg(N_int,2)
   integer*8                                :: IJsomo
   integer*8                                :: diffSOMO
   integer*8                                :: diffDOMO
@@ -209,9 +242,12 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
   integer                        :: iii,ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes
   integer                        :: listholes(mo_num)
   integer                        :: holetype(mo_num)
-  integer                        :: end_index
+  integer                        :: end_index, ishift
-  integer                        :: Nsomo_alpha
+  integer                        :: Nsomo_alpha, pp,qq, nperm
   integer*8                      :: MS
+  integer                        :: exc(0:2,2,2), tz, m, n, high, low
+  integer                        :: listall(N_int*bit_kind_size), nelall
+  integer(bit_kind)              :: hole, particle, tmp
   MS = elec_alpha_num-elec_beta_num
 
   nconnectedI = 0
@@ -219,42 +255,66 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
 
   ! Since CFGs are sorted wrt to seniority
   ! we don't have to search the full CFG list
-  Isomo = Ialpha(1,1)
+  !Isomo = Ialpha(1,1)
-  Idomo = Ialpha(1,2)
+  !Idomo = Ialpha(1,2)
-  Nsomo_alpha = POPCNT(Isomo)
+  !Nsomo_alpha = POPCNT(Isomo)
+  Icfg = Ialpha
+  Nsomo_alpha = 0
+  do i=1,N_int
+    Isomo = Ialpha(i,1)
+    Idomo = Ialpha(i,2)
+    Nsomo_alpha += POPCNT(Isomo)
+  end do
   end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_alpha+4,elec_num))-1)
   if(end_index .LT. 0) end_index= N_configuration
-  end_index = N_configuration
+  !end_index = N_configuration
 
 
   p = 0
   q = 0
-  if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1'
+  !if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1'
   do i=idxI,end_index
-     Isomo = Ialpha(1,1)
-     Idomo = Ialpha(1,2)
-     Jsomo = psi_configuration(1,1,i)
-     Jdomo = psi_configuration(1,2,i)
      ! Check for Minimal alpha electrons (MS)
-     if(POPCNT(Isomo).lt.MS)then
+     if(Nsomo_alpha .lt. MS)then
        cycle
      endif
-     diffSOMO = IEOR(Isomo,Jsomo)
+     !Isomo = Ialpha(1,1)
-     ndiffSOMO = POPCNT(diffSOMO)
+     !Idomo = Ialpha(1,2)
-     !if(idxI.eq.1)then
+     !Jsomo = psi_configuration(1,1,i)
-     !  print *," \t idxI=",i," diffS=",ndiffSOMO," popJs=", POPCNT(Jsomo)," popIs=",POPCNT(Isomo)
+     !Jdomo = psi_configuration(1,2,i)
-     !endif
+     !diffSOMO = IEOR(Isomo,Jsomo)
-     diffDOMO = IEOR(Idomo,Jdomo)
+     !ndiffSOMO = POPCNT(diffSOMO)
-     xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+     !diffDOMO = IEOR(Idomo,Jdomo)
-     ndiffDOMO = POPCNT(diffDOMO)
+     !xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
-     nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+     !ndiffDOMO = POPCNT(diffDOMO)
-     nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO 
+     !nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+     !nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO 
+
+     ndiffSOMO = 0
+     ndiffDOMO = 0
+     nxordiffSOMODOMO = 0
+     do ii=1,N_int
+       Isomo = Ialpha(ii,1)
+       Idomo = Ialpha(ii,2)
+       Jsomo = psi_configuration(ii,1,i)
+       Jdomo = psi_configuration(ii,2,i)
+       diffSOMO = IEOR(Isomo,Jsomo)
+       ndiffSOMO += POPCNT(diffSOMO)
+       diffDOMO = IEOR(Idomo,Jdomo)
+       xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+       ndiffDOMO += POPCNT(diffDOMO)
+       nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
+       nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO 
+     end do
+     Jcfg = psi_configuration(:,:,i)
+
      if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
         select case(ndiffDOMO)
         case (0)
            ! SOMO -> VMO
            !print *,"obt SOMO -> VMO"
            extyp = 3
+           if(N_int .eq. 1) then
            IJsomo = IEOR(Isomo, Jsomo)
 !IRP_IF WITHOUT_TRAILZ
 !           p = (popcnt(ieor( IAND(Isomo,IJsomo) , IAND(Isomo,IJsomo) -1))-1) + 1
@@ -267,6 +327,77 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
 !IRP_ELSE
            q = TRAILZ(IJsomo) + 1
 !IRP_ENDIF
+           !print *," p=",p," q=",q
+           !call get_single_excitation_cfg(Jcfg, Icfg, p, q, N_int)
+         else
+  exc = 0
+  do ii = 1,2
+    ishift = 1-bit_kind_size
+    do l=1,N_int
+      ishift = ishift + bit_kind_size
+      if (Jcfg(l,ii) == Icfg(l,ii)) then
+        cycle
+      endif
+      tmp = xor( Jcfg(l,ii), Icfg(l,ii) )
+      particle = iand(tmp, Icfg(l,ii))
+      hole     = iand(tmp, Jcfg(l,ii))
+      if (particle /= 0_bit_kind) then
+        tz = trailz(particle)
+        exc(0,2,ii) = 1
+        exc(1,2,ii) = tz+ishift
+        !print *,"part ",tz+ishift, " ii=",ii, exc(1,2,2)
+      endif
+      if (hole /= 0_bit_kind) then
+        tz = trailz(hole)
+        exc(0,1,ii) = 1
+        exc(1,1,ii) = tz+ishift
+        !print *,"hole ",tz+ishift, " ii=",ii, exc(1,1,2)
+      endif
+
+      if ( iand(exc(0,1,ii),exc(0,2,ii)) /= 1) then  ! exc(0,1,ii)/=1 and exc(0,2,ii) /= 1
+        cycle
+      endif
+
+      high = max(exc(1,1,ii), exc(1,2,ii))-1
+      low  = min(exc(1,1,ii), exc(1,2,ii))
+
+      ASSERT (low >= 0)
+      ASSERT (high > 0)
+
+      k = shiftr(high,bit_kind_shift)+1
+      j = shiftr(low,bit_kind_shift)+1
+      m = iand(high,bit_kind_size-1)
+      n = iand(low,bit_kind_size-1)
+
+      if (j==k) then
+        nperm = nperm + popcnt(iand(Jcfg(j,ii),           &
+            iand( shiftl(1_bit_kind,m)-1_bit_kind,            &
+                  not(shiftl(1_bit_kind,n))+1_bit_kind)) )
+      else
+        nperm = nperm + popcnt(                                    &
+             iand(Jcfg(j,ii),                                   &
+                  iand(not(0_bit_kind),                            &
+                       (not(shiftl(1_bit_kind,n)) + 1_bit_kind) ))) &
+             + popcnt(iand(Jcfg(k,ii),                          &
+                           (shiftl(1_bit_kind,m) - 1_bit_kind ) ))
+
+        do iii=j+1,k-1
+          nperm = nperm + popcnt(Jcfg(iii,ii))
+        end do
+
+      endif
+
+      ! Set p and q
+      q = max(exc(1,1,1),exc(1,1,2))
+      p = max(exc(1,2,1),exc(1,2,2))
+      exit
+
+    enddo
+  enddo
+endif
+           !assert ( p == pp)
+           !assert ( q == qq)
+           !print *," --- p=",p," q=",q
         case (1)
            ! DOMO -> VMO
            ! or
@@ -277,6 +408,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
               ! DOMO -> VMO
               !print *,"obt DOMO -> VMO"
               extyp = 2
+              if(N_int.eq.1)then
 !IRP_IF WITHOUT_TRAILZ
 !              p = (popcnt(ieor( IEOR(Idomo,Jdomo),IEOR(Idomo,Jdomo) -1))-1) + 1
 !IRP_ELSE
@@ -289,10 +421,83 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
 !IRP_ELSE
               q = TRAILZ(Isomo) + 1
 !IRP_ENDIF
+  else
+    exc=0
+  exc(0,1,1) = 0
+  exc(0,2,1) = 0
+  exc(0,1,2) = 0
+  exc(0,2,2) = 0
+  do ii = 1,2
+    ishift = 1-bit_kind_size
+    do l=1,N_int
+      ishift = ishift + bit_kind_size
+      if (Jcfg(l,ii) == Icfg(l,ii)) then
+        cycle
+      endif
+      tmp = xor( Jcfg(l,ii), Icfg(l,ii) )
+      particle = iand(tmp, Icfg(l,ii))
+      hole     = iand(tmp, Jcfg(l,ii))
+      if (particle /= 0_bit_kind) then
+        tz = trailz(particle)
+        exc(0,2,ii) = 1
+        exc(1,2,ii) = tz+ishift
+        !print *,"part ",tz+ishift, " ii=",ii
+      endif
+      if (hole /= 0_bit_kind) then
+        tz = trailz(hole)
+        exc(0,1,ii) = 1
+        exc(1,1,ii) = tz+ishift
+        !print *,"hole ",tz+ishift, " ii=",ii
+      endif
+
+      if ( iand(exc(0,1,ii),exc(0,2,ii)) /= 1) then  ! exc(0,1,ii)/=1 and exc(0,2,ii) /= 1
+        cycle
+      endif
+
+      high = max(exc(1,1,ii), exc(1,2,ii))-1
+      low  = min(exc(1,1,ii), exc(1,2,ii))
+
+      ASSERT (low >= 0)
+      ASSERT (high > 0)
+
+      k = shiftr(high,bit_kind_shift)+1
+      j = shiftr(low,bit_kind_shift)+1
+      m = iand(high,bit_kind_size-1)
+      n = iand(low,bit_kind_size-1)
+
+      if (j==k) then
+        nperm = nperm + popcnt(iand(Jcfg(j,ii),           &
+            iand( shiftl(1_bit_kind,m)-1_bit_kind,            &
+                  not(shiftl(1_bit_kind,n))+1_bit_kind)) )
+      else
+        nperm = nperm + popcnt(                                    &
+             iand(Jcfg(j,ii),                                   &
+                  iand(not(0_bit_kind),                            &
+                       (not(shiftl(1_bit_kind,n)) + 1_bit_kind) ))) &
+             + popcnt(iand(Jcfg(k,ii),                          &
+                           (shiftl(1_bit_kind,m) - 1_bit_kind ) ))
+
+        do iii=j+1,k-1
+          nperm = nperm + popcnt(Jcfg(iii,ii))
+        end do
+
+      endif
+
+      ! Set p and q
+      q = max(exc(1,1,1),exc(1,1,2))
+      p = max(exc(1,2,1),exc(1,2,2))
+      exit
+
+    enddo
+  enddo
+endif
+           !assert ( p == pp)
+           !assert ( q == qq)
            else
               ! SOMO -> SOMO
               !print *,"obt SOMO -> SOMO"
               extyp = 1
+              if(N_int.eq.1)then
 !IRP_IF WITHOUT_TRAILZ
 !              q = (popcnt(ieor( IEOR(Idomo,Jdomo), IEOR(Idomo,Jdomo)-1))-1) + 1
 !IRP_ELSE
@@ -309,11 +514,84 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
               !if(POPCNT(Isomo).lt.MS)then
               !  cycle
               !endif
+  else
+  exc=0
+  exc(0,1,1) = 0
+  exc(0,2,1) = 0
+  exc(0,1,2) = 0
+  exc(0,2,2) = 0
+  do ii = 1,2
+    ishift = 1-bit_kind_size
+    do l=1,N_int
+      ishift = ishift + bit_kind_size
+      if (Jcfg(l,ii) == Icfg(l,ii)) then
+        cycle
+      endif
+      tmp = xor( Jcfg(l,ii), Icfg(l,ii) )
+      particle = iand(tmp, Icfg(l,ii))
+      hole     = iand(tmp, Jcfg(l,ii))
+      if (particle /= 0_bit_kind) then
+        tz = trailz(particle)
+        exc(0,2,ii) = 1
+        exc(1,2,ii) = tz+ishift
+        !print *,"part ",tz+ishift, " ii=",ii
+      endif
+      if (hole /= 0_bit_kind) then
+        tz = trailz(hole)
+        exc(0,1,ii) = 1
+        exc(1,1,ii) = tz+ishift
+        !print *,"hole ",tz+ishift, " ii=",ii
+      endif
+
+      if ( iand(exc(0,1,ii),exc(0,2,ii)) /= 1) then  ! exc(0,1,ii)/=1 and exc(0,2,ii) /= 1
+        cycle
+      endif
+
+      high = max(exc(1,1,ii), exc(1,2,ii))-1
+      low  = min(exc(1,1,ii), exc(1,2,ii))
+
+      ASSERT (low >= 0)
+      ASSERT (high > 0)
+
+      k = shiftr(high,bit_kind_shift)+1
+      j = shiftr(low,bit_kind_shift)+1
+      m = iand(high,bit_kind_size-1)
+      n = iand(low,bit_kind_size-1)
+
+      if (j==k) then
+        nperm = nperm + popcnt(iand(Jcfg(j,ii),           &
+            iand( shiftl(1_bit_kind,m)-1_bit_kind,            &
+                  not(shiftl(1_bit_kind,n))+1_bit_kind)) )
+      else
+        nperm = nperm + popcnt(                                    &
+             iand(Jcfg(j,ii),                                   &
+                  iand(not(0_bit_kind),                            &
+                       (not(shiftl(1_bit_kind,n)) + 1_bit_kind) ))) &
+             + popcnt(iand(Jcfg(k,ii),                          &
+                           (shiftl(1_bit_kind,m) - 1_bit_kind ) ))
+
+        do iii=j+1,k-1
+          nperm = nperm + popcnt(Jcfg(iii,ii))
+        end do
+
+      endif
+
+      ! Set p and q
+      q = max(exc(1,1,1),exc(1,1,2))
+      p = max(exc(1,2,1),exc(1,2,2))
+      exit
+
+    enddo
+  enddo
+endif
+           !assert ( p == pp)
+           !assert ( q == qq)
            end if
         case (2)
            ! DOMO -> SOMO
            !print *,"obt DOMO -> SOMO"
            extyp = 4
+           if(N_int.eq.1)then
            IJsomo = IEOR(Isomo, Jsomo)
 !IRP_IF WITHOUT_TRAILZ
 !           p = (popcnt(ieor( IAND(Jsomo,IJsomo), IAND(Jsomo,IJsomo)-1))-1) + 1
@@ -326,6 +604,79 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
 !IRP_ELSE
            q = TRAILZ(IJsomo) + 1
 !IRP_ENDIF
+
+  else
+  exc=0
+  exc(0,1,1) = 0
+  exc(0,2,1) = 0
+  exc(0,1,2) = 0
+  exc(0,2,2) = 0
+  do ii = 1,2
+    ishift = 1-bit_kind_size
+    do l=1,N_int
+      ishift = ishift + bit_kind_size
+      if (Jcfg(l,ii) == Icfg(l,ii)) then
+        cycle
+      endif
+      tmp = xor( Jcfg(l,ii), Icfg(l,ii) )
+      particle = iand(tmp, Icfg(l,ii))
+      hole     = iand(tmp, Jcfg(l,ii))
+      if (particle /= 0_bit_kind) then
+        tz = trailz(particle)
+        exc(0,2,ii) = 1
+        exc(1,2,ii) = tz+ishift
+        !print *,"part ",tz+ishift, " ii=",ii
+      endif
+      if (hole /= 0_bit_kind) then
+        tz = trailz(hole)
+        exc(0,1,ii) = 1
+        exc(1,1,ii) = tz+ishift
+        !print *,"hole ",tz+ishift, " ii=",ii
+      endif
+
+      if ( iand(exc(0,1,ii),exc(0,2,ii)) /= 1) then  ! exc(0,1,ii)/=1 and exc(0,2,ii) /= 1
+        cycle
+      endif
+
+      high = max(exc(1,1,ii), exc(1,2,ii))-1
+      low  = min(exc(1,1,ii), exc(1,2,ii))
+
+      ASSERT (low >= 0)
+      ASSERT (high > 0)
+
+      k = shiftr(high,bit_kind_shift)+1
+      j = shiftr(low,bit_kind_shift)+1
+      m = iand(high,bit_kind_size-1)
+      n = iand(low,bit_kind_size-1)
+
+      if (j==k) then
+        nperm = nperm + popcnt(iand(Jcfg(j,ii),           &
+            iand( shiftl(1_bit_kind,m)-1_bit_kind,            &
+                  not(shiftl(1_bit_kind,n))+1_bit_kind)) )
+      else
+        nperm = nperm + popcnt(                                    &
+             iand(Jcfg(j,ii),                                   &
+                  iand(not(0_bit_kind),                            &
+                       (not(shiftl(1_bit_kind,n)) + 1_bit_kind) ))) &
+             + popcnt(iand(Jcfg(k,ii),                          &
+                           (shiftl(1_bit_kind,m) - 1_bit_kind ) ))
+
+        do iii=j+1,k-1
+          nperm = nperm + popcnt(Jcfg(iii,ii))
+        end do
+
+      endif
+
+      ! Set p and q
+      q = max(exc(1,1,1),exc(1,1,2))
+      p = max(exc(1,2,1),exc(1,2,2))
+      exit
+
+    enddo
+  enddo
+endif
+           !assert ( p == pp)
+           !assert ( q == qq)
         case default
            print *,"something went wront in get connectedI"
         end select
@@ -345,26 +696,46 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
         ! find out all pq holes possible
         nholes = 0
         ! holes in SOMO
-        Isomo = psi_configuration(1,1,i)
+        !Isomo = psi_configuration(1,1,i)
-        Idomo = psi_configuration(1,2,i)
+        !Idomo = psi_configuration(1,2,i)
-        do iii = 1,n_act_orb
+        !do iii = 1,n_act_orb
-          ii = list_act(iii)
+        !  ii = list_act(iii)
-           if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
+        !   if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
-              nholes += 1
+        !      nholes += 1
-              listholes(nholes) = ii
+        !      listholes(nholes) = ii
-              holetype(nholes) = 1
+        !      holetype(nholes) = 1
-           endif
+        !   endif
+        !end do
+        call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int)
+
+        do iii=1,nelall
+          nholes += 1
+          listholes(nholes) = listall(iii)
+          holetype(nholes) = 1
         end do
+
         ! holes in DOMO
-        do iii = 1,n_act_orb
+        !do iii = 1,n_act_orb
-          ii = list_act(iii)
+        !  ii = list_act(iii)
-           if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
+        !   if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
-              nholes += 1
+        !      nholes += 1
-              listholes(nholes) = ii
+        !      listholes(nholes) = ii
-              holetype(nholes) = 2
+        !      holetype(nholes) = 2
-           endif
+        !   endif
+        !end do
+        nelall=0
+        listall=0
+        call bitstring_to_list(psi_configuration(1,2,i),listall,nelall,N_int)
+
+        do iii=1,nelall
+          if(listall(iii) .gt. n_core_orb)then
+            nholes += 1
+            listholes(nholes) = listall(iii)
+            holetype(nholes) = 2
+          endif
         end do
 
+
         do k=1,nholes
            p = listholes(k)
            q = p

src/csf/sigma_vector.irp.f

@@ -835,7 +835,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
   ! the configurations in psi_configuration
   ! returns : diag_energies :
   END_DOC
-  integer :: i,j,k,kk,l,p,q,noccp,noccq, ii, jj
+  integer :: i,j,k,kk,l,p,q,noccp,noccq, ii, jj, iii
   real*8,intent(out) :: diag_energies(n_CSF)
   integer                            :: nholes
   integer                            :: nvmos
@@ -863,6 +863,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
   real*8            :: meCC
   real*8            :: ecore
   real*8            :: core_act_contrib
+  integer                        :: listall(N_int*bit_kind_size), nelall
 
   !PROVIDE h_core_ri
   PROVIDE core_fock_operator
@@ -894,47 +895,61 @@ subroutine calculate_preconditioner_cfg(diag_energies)
      ! find out all pq holes possible
      nholes = 0
      ! holes in SOMO
-     !do k = 1,mo_num
+     !do kk = 1,n_act_orb
-     do kk = 1,n_act_orb
+     !  k = list_act(kk)
-       k = list_act(kk)
+     !   if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
-        if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
+     !      nholes += 1
-           nholes += 1
+     !      listholes(nholes) = k
-           listholes(nholes) = k
+     !      holetype(nholes) = 1
-           holetype(nholes) = 1
+     !   endif
-        endif
+     !enddo
-     enddo
+     call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int)
-     ! holes in DOMO
-     !do k = n_core_orb+1,n_core_orb + n_act_orb
-     !do k = 1+n_core_inact_orb,n_core_orb+n_core_inact_act_orb
-     !do k = 1,mo_num
-     do kk = 1,n_act_orb
-       k = list_act(kk)
-        if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
-           nholes += 1
-           listholes(nholes) = k
-           holetype(nholes) = 2
-        endif
-     enddo
 
-     ! find vmos
+     do iii=1,nelall
-     listvmos = -1
+       nholes += 1
-     vmotype = -1
+       listholes(nholes) = listall(iii)
-     nvmos = 0
+       holetype(nholes) = 1
-     !do k = n_core_orb+1,n_core_orb + n_act_orb
+     end do
-     !do k = 1,mo_num
+
-     do kk = 1,n_act_orb
+     ! holes in DOMO
-       k = list_act(kk)
+     !do kk = 1,n_act_orb
-        !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
+     !  k = list_act(kk)
-        if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0) then
+     !   if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
-           nvmos += 1
+     !      nholes += 1
-           listvmos(nvmos) = k
+     !      listholes(nholes) = k
-           vmotype(nvmos) = 0
+     !      holetype(nholes) = 2
-        else if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0 ) then
+     !   endif
-           nvmos += 1
+     !enddo
-           listvmos(nvmos) = k
+     call bitstring_to_list(psi_configuration(1,2,i),listall,nelall,N_int)
-           vmotype(nvmos) = 1
+
-        end if
+     do iii=1,nelall
-     enddo
+       if(listall(iii) .gt. n_core_orb)then
+         nholes += 1
+         listholes(nholes) = listall(iii)
+         holetype(nholes) = 2
+       endif
+     end do
+
+
+     !!! 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
+     !!do kk = 1,n_act_orb
+     !!  k = list_act(kk)
+     !!   !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
      !print *,"I=",i
      !call debug_spindet(psi_configuration(1,1,i),N_int)
      !call debug_spindet(psi_configuration(1,2,i),N_int)
@@ -1413,8 +1428,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
   !nconnectedtotalmax = 1000
   !nconnectedmaxJ = 1000
   maxnalphas = elec_num*mo_num
-  Icfg(1,1) = psi_configuration(1,1,1)
+  Icfg(:,1) = psi_configuration(:,1,1)
-  Icfg(1,2) = psi_configuration(1,2,1)
+  Icfg(:,2) = psi_configuration(:,2,1)
   allocate(listconnectedJ(N_INT,2,max(sze,10000)))
   allocate(idslistconnectedJ(max(sze,10000)))
   call obtain_connected_J_givenI(1, Icfg, listconnectedJ, idslistconnectedJ, nconnectedmaxJ, nconnectedtotalmax)
@@ -1632,9 +1647,9 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
 
      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)
-     if(Nalphas_Icfg .GT. maxnalphas) then
+     !if(Nalphas_Icfg .GT. maxnalphas) then
-       print *,"Nalpha > maxnalpha"
+     !  print *,"Nalpha > maxnalpha"
-     endif
+     !endif
 
      call obtain_connected_J_givenI(i, Icfg, listconnectedJ, idslistconnectedJ, nconnectedJ, ntotJ)
 
@@ -1650,15 +1665,15 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
         call obtain_connected_I_foralpha(i, alphas_Icfg(1,1,k), connectedI_alpha, idxs_connectedI_alpha, &
                                          nconnectedI, excitationIds, excitationTypes, diagfactors)
 
+        !if(i .EQ. 1) then
+        !   print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' kcfgDOMO=',alphas_Icfg(1,2,k),' ',POPCNT(alphas_Icfg(1,2,k))
+        !endif
+
         
         if(nconnectedI .EQ. 0) then
            cycle
         endif
 
-        !if(i .EQ. 1) then
-        !   print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' kcfgDOMO=',alphas_Icfg(1,2,k),' ',POPCNT(alphas_Icfg(1,2,k))
-        !endif
-
         ! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
         totcolsTKI = 0
         rowsTKI = -1

src/determinants/slater_rules.irp.f

@@ -83,7 +83,7 @@ subroutine get_excitation(det1,det2,exc,degree,phase,Nint)
   !               exc(1,1,1) = q
   !               exc(1,2,1) = p
 
-  ! T^alpha_pq  : exc(0,1,2) = 1
+  ! T^beta_pq   : exc(0,1,2) = 1
   !               exc(0,2,2) = 1
   !               exc(1,1,2) = q
   !               exc(1,2,2) = p
@@ -434,6 +434,98 @@ subroutine get_single_excitation(det1,det2,exc,phase,Nint)
 
 end
 
+subroutine get_single_excitation_cfg(cfg1,cfg2,p,q,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns the excitation operator between two singly excited configurations.
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: cfg1(Nint,2)
+  integer(bit_kind), intent(in)  :: cfg2(Nint,2)
+  integer, intent(out)           :: p, q
+  integer                        :: tz
+  integer                        :: l, ispin, idx_hole, idx_particle, ishift
+  integer                        :: nperm
+  integer                        :: i,j,k,m,n
+  integer                        :: high, low
+  integer                        :: a,b,c,d
+  integer(bit_kind)              :: hole, particle, tmp
+  integer                        :: exc(0:2,2,2)
+
+  ASSERT (Nint > 0)
+  nperm = 0
+  p = 0
+  q = 0
+  exc(0,1,1) = 0
+  exc(0,2,1) = 0
+  exc(0,1,2) = 0
+  exc(0,2,2) = 0
+  do ispin = 1,2
+    ishift = 1-bit_kind_size
+    do l=1,Nint
+      ishift = ishift + bit_kind_size
+      if (cfg1(l,ispin) == cfg2(l,ispin)) then
+        cycle
+      endif
+      tmp = xor( cfg1(l,ispin), cfg2(l,ispin) )
+      particle = iand(tmp, cfg2(l,ispin))
+      hole     = iand(tmp, cfg1(l,ispin))
+      if (particle /= 0_bit_kind) then
+        tz = trailz(particle)
+        exc(0,2,ispin) = 1
+        exc(1,2,ispin) = tz+ishift
+        !print *,"part ",tz+ishift, " ispin=",ispin
+      endif
+      if (hole /= 0_bit_kind) then
+        tz = trailz(hole)
+        exc(0,1,ispin) = 1
+        exc(1,1,ispin) = tz+ishift
+        !print *,"hole ",tz+ishift, " ispin=",ispin
+      endif
+
+      if ( iand(exc(0,1,ispin),exc(0,2,ispin)) /= 1) then  ! exc(0,1,ispin)/=1 and exc(0,2,ispin) /= 1
+        cycle
+      endif
+
+      high = max(exc(1,1,ispin), exc(1,2,ispin))-1
+      low  = min(exc(1,1,ispin), exc(1,2,ispin))
+
+      ASSERT (low >= 0)
+      ASSERT (high > 0)
+
+      k = shiftr(high,bit_kind_shift)+1
+      j = shiftr(low,bit_kind_shift)+1
+      m = iand(high,bit_kind_size-1)
+      n = iand(low,bit_kind_size-1)
+
+      if (j==k) then
+        nperm = nperm + popcnt(iand(cfg1(j,ispin),           &
+            iand( shiftl(1_bit_kind,m)-1_bit_kind,            &
+                  not(shiftl(1_bit_kind,n))+1_bit_kind)) )
+      else
+        nperm = nperm + popcnt(                                    &
+             iand(cfg1(j,ispin),                                   &
+                  iand(not(0_bit_kind),                            &
+                       (not(shiftl(1_bit_kind,n)) + 1_bit_kind) ))) &
+             + popcnt(iand(cfg1(k,ispin),                          &
+                           (shiftl(1_bit_kind,m) - 1_bit_kind ) ))
+
+        do i=j+1,k-1
+          nperm = nperm + popcnt(cfg1(i,ispin))
+        end do
+
+      endif
+
+      ! Set p and q
+      q = max(exc(1,1,1),exc(1,1,2))
+      p = max(exc(1,2,1),exc(1,2,2))
+      return
+
+    enddo
+  enddo
+end
+
 subroutine bitstring_to_list_ab( string, list, n_elements, Nint)
   use bitmasks
   implicit none