Added a function to perform a single excitation on cfg and identify the type.

src/determinants/create_excitations.irp.f

@@ -183,6 +183,118 @@ subroutine do_single_excitation_cfg(key_in,key_out,i_hole,i_particle,ok)
 
 end
 
+subroutine do_single_excitation_cfg_with_type(key_in,key_out,i_hole,i_particle,ex_type,ok)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Applies the single excitation operator to a configuration
+  ! Returns the type of excitation in ex_type
+  ! where the following convention is used
+  ! 1 = (SOMO -> SOMO) 1 change in Nsomo
+  ! 2 = (DOMO -> VMO)  1 change in Nsomo
+  ! 3 = (SOMO -> VMO)  0 change in Nsomo
+  ! 4 = (DOMO -> SOMO) 0 change in Nsomo
+  ! If the excitation is possible, ok is True
+  END_DOC
+  integer, intent(in)            :: i_hole,i_particle
+  integer(bit_kind), intent(in)  :: key_in(N_int,2)
+  integer , intent(out)          :: ex_type
+  logical , intent(out)          :: ok
+  integer                        :: k,j,i
+  integer(bit_kind)              :: mask
+  integer(bit_kind)              :: key_out(N_int,2)
+  logical                        :: isholeSOMO
+  logical                        :: isparticleSOMO
+  logical                        :: isholeDOMO
+  logical                        :: isparticleVMO
+
+  ASSERT (i_hole > 0)
+  ASSERT (i_particle <= mo_num)
+
+  ok = .True.
+  key_out(:,:) = key_in(:,:)
+
+  ! hole
+  k = shiftr(i_hole-1,bit_kind_shift)+1
+  j = i_hole-shiftl(k-1,bit_kind_shift)-1
+  mask = ibset(0_bit_kind,j)
+
+  ! Check if the position j is singly occupied
+  ! 1  ->  0  (SOMO)
+  ! 0      0  (DOMO)
+  if (iand(key_out(k,1),mask) /= 0_bit_kind) then
+     key_out(k,1) = ibclr(key_out(k,1),j)
+     isholeSOMO = .True.
+
+  ! Check if the position j is doubly occupied
+  ! 0  ->  1  (SOMO)
+  ! 1      0  (DOMO)
+  else if (iand(key_out(k,2),mask) /= 0_bit_kind) then
+     key_out(k,1) = ibset(key_out(k,1),j)
+     key_out(k,2) = ibclr(key_out(k,2),j)
+     isholeDOMO = .True.
+
+  ! The position j is unoccupied: Not OK
+  ! 0  ->  0  (SOMO)
+  ! 0      0  (DOMO)
+  else
+     ok =.False.
+     return
+  endif
+
+
+  ! particle
+  k = shiftr(i_particle-1,bit_kind_shift)+1
+  j = i_particle-shiftl(k-1,bit_kind_shift)-1
+  mask = ibset(0_bit_kind,j)
+
+  ! Check if the position j is singly occupied
+  ! 1  ->  0  (SOMO)
+  ! 0      1  (DOMO)
+  if (iand(key_out(k,1),mask) /= 0_bit_kind) then
+     key_out(k,1) = ibclr(key_out(k,1),j)
+     key_out(k,2) = ibset(key_out(k,2),j)
+     isparticleSOMO = .True.
+
+  ! Check if the position j is doubly occupied : Not OK
+  ! 0  ->  1  (SOMO)
+  ! 1      0  (DOMO)
+  else if (iand(key_out(k,2),mask) /= 0_bit_kind) then
+     ok = .False.
+     return
+
+  ! Position at j is unoccupied
+  ! 0  ->  0  (SOMO)
+  ! 0      0  (DOMO)
+  else
+     key_out(k,1) = ibset(key_out(k,1),j)
+     isparticleVMO = .True.
+  endif
+
+  if(isholeSOMO) then
+      ! two possibilities
+      ! particle is SOMO or VMO
+      if(isparticleSOMO) then
+          ! SOMO -> SOMO
+          ex_type = 1
+      else
+          ! SOMO -> VMO
+          ex_type = 3
+      endif 
+   else
+       ! two possibilities
+       ! particle is SOMO or VMO
+       if(isparticleSOMO) then
+           ! DOMO -> SOMO
+           ex_type = 4
+       else
+           ! DOMO -> VMO
+           ex_type = 2
+       endif
+    endif
+
+end
+
 subroutine generate_all_singles_cfg(cfg,singles,n_singles,Nint)
   implicit none
   use bitmasks