better isInCassd

plugins/mrcepa0/dressing.irp.f

@@ -55,7 +55,12 @@ use bitmasks
 END_PROVIDER
 
 
-
+BEGIN_PROVIDER [ integer, HP, (2,N_det_non_ref) ]
+  integer :: i
+  do i=1,N_det_non_ref
+    call getHP(psi_non_ref(1,1,i), HP(1,i), HP(2,i), N_int)
+  end do 
+END_PROVIDER
 
  BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ]
 &BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ]
@@ -94,6 +99,7 @@ END_PROVIDER
     end do
   end if
   
+  
   do i=1, N_det_non_ref
   do k=1, N_int
     det_noactive(k,1,i) = iand(psi_non_ref(k,1,i), not(active_sorb(k,1)))
@@ -277,6 +283,7 @@ logical function detEq(a,b,Nint)
    detEq = .true.
 end function
 
+
 logical function isInCassd(a,Nint)
   use bitmasks
   implicit none
@@ -288,7 +295,6 @@ logical function isInCassd(a,Nint)
   
   isInCassd = .false.
      
-   
   deg = 0
    do i=1,2
     do ni=1,Nint
@@ -306,10 +312,46 @@ logical function isInCassd(a,Nint)
       if(deg > 2) return
     end do
   end do
-  
   isInCassd = .true.
 end function
 
+
+subroutine getHP(a,h,p,Nint)
+  use bitmasks
+  implicit none
+  integer, intent(in) :: Nint
+  integer(bit_kind), intent(in) :: a(Nint,2)
+  integer, intent(out) :: h, p
+  integer(bit_kind) :: inac, virt
+  integer :: ni, i, deg
+
+  
+  !isInCassd = .false.
+  h = 0
+  p = 0
+   
+  deg = 0
+  lp : do i=1,2
+    do ni=1,Nint
+      virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i)))
+      deg += popcnt(iand(virt, a(ni,i)))
+      if(deg > 2) exit lp
+    end do
+  end do lp
+  p = deg
+
+  deg = 0
+  lh : do i=1,2
+    do ni=1,Nint
+      inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i)))
+      deg += popcnt(xor(iand(inac,a(ni,i)), inac))
+      if(deg > 2) exit lh
+    end do
+  end do lh
+  h = deg
+  !isInCassd = .true.
+end function
+
 integer function detCmp(a,b,Nint)
    use bitmasks
    implicit none

plugins/mrcepa0/dressing_slave.irp.f

@@ -158,12 +158,15 @@ subroutine mrsc2_dressing_slave(thread,iproc)
         
         i = komon(m)
         
+        if(HP(1,i) + HP(1,k) <= 2 .and. HP(2,i) + HP(2,k) <= 2) cycle
+
         hJi = h_(J,i)
         hIi = h_(i_I,i)
         
+
         call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
         if(.not. ok) cycle
-        if(isInCassd(det_tmp, N_int)) cycle
+        !if(isInCassd(det_tmp, N_int)) cycle
           
         do i_state = 1, N_states 
           if(lambda_mrcc(i_state, i) == 0d0) cycle