src/bitmask/core_inact_act_virt.irp.f

@@ -322,7 +322,6 @@ END_PROVIDER
    enddo
    print *,  'Active MOs:'
    print *,  list_act(1:n_act_orb)
-   print*,   list_act_reverse(1:n_act_orb)
    
 END_PROVIDER
  

src/two_body_rdm/all_2rdm_routines.irp.f

@@ -392,7 +392,7 @@ subroutine all_two_rdm_dm_nstates_work_$N_int(big_array_aa,big_array_bb,big_arra
          c_1(l) = u_t(l,l_a)
          c_2(l) = u_t(l,k_a)
        enddo
-       call off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
+       call off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
        ASSERT (l_a <= N_det)
        
      enddo

src/two_body_rdm/all_states_routines.irp.f

@@ -442,7 +442,7 @@ subroutine orb_range_all_states_two_rdm_work_$N_int(big_array,dim1,norb,list_orb
           c_2(l) = u_t(l,k_a)
           c_contrib(l) = c_1(l) * c_2(l) 
         enddo
-        call orb_range_off_diagonal_double_to_two_rdm_bb_dm_all_states(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_contrib,N_st,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
+        call orb_range_off_diagonal_double_to_two_rdm_bb_dm_all_states(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_contrib,N_st,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
         ASSERT (l_a <= N_det)
         
       enddo

src/two_body_rdm/orb_range_routines.irp.f

@@ -31,7 +31,6 @@ subroutine orb_range_two_rdm_state_av(big_array,dim1,norb,list_orb,list_orb_reve
        size(u_t, 1),                                                 &
        N_det, N_st)
    
-   
    call orb_range_two_rdm_state_av_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1)
    deallocate(u_t)
    
@@ -136,7 +135,6 @@ subroutine orb_range_two_rdm_state_av_work_$N_int(big_array,dim1,norb,list_orb,l
     stop
    endif
    
-     
    PROVIDE N_int
 
    call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
@@ -445,7 +443,7 @@ subroutine orb_range_two_rdm_state_av_work_$N_int(big_array,dim1,norb,list_orb,l
           c_2(l) = u_t(l,k_a)
           c_average += c_1(l) * c_2(l) * state_weights(l)
         enddo
-        call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
+        call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
         ASSERT (l_a <= N_det)
         
       enddo

src/two_body_rdm/orb_range_routines_openmp.irp.f

@@ -93,9 +93,11 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
    double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
    
    integer                        :: i,j,k,l
-   integer                        :: k_a, k_b, l_a, l_b
+   integer                        :: k_a, k_b, l_a, l_b, m_a, m_b
-   integer                        :: krow, kcol
+   integer                        :: istate
+   integer                        :: krow, kcol, krow_b, kcol_b
    integer                        :: lrow, lcol
+   integer                        :: mrow, mcol
    integer(bit_kind)              :: spindet($N_int)
    integer(bit_kind)              :: tmp_det($N_int,2)
    integer(bit_kind)              :: tmp_det2($N_int,2)
@@ -107,6 +109,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
    integer, allocatable           :: singles_b(:)
    integer, allocatable           :: idx(:), idx0(:)
    integer                        :: maxab, n_singles_a, n_singles_b, kcol_prev
+   integer*8                      :: k8
    double precision               :: c_average
 
    logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@@ -133,11 +136,16 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
     stop
    endif
 
+  !do i = 1, N_int
+  ! det_1_act(i,1) =   iand(det_1(i,1),orb_bitmask(i))
+  ! det_1_act(i,2) =   iand(det_1(i,2),orb_bitmask(i))
+  !enddo
+
    
    PROVIDE N_int
 
    call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
-   sze_buff = norb ** 3 + 6 * norb
+   sze_buff = norb ** 3
    list_orb_reverse = -1000 
    do i = 1, norb
     list_orb_reverse(list_orb(i)) = i 
@@ -165,13 +173,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        !    !$OMP          psi_bilinear_matrix_columns_loc,          &
        !    !$OMP          psi_bilinear_matrix_transp_rows_loc,      &
        !    !$OMP          istart, iend, istep, irp_here, v_t, s_t,  &
-       !    !$OMP          ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin)     &
+       !    !$OMP          ishift, idx0, u_t, maxab)                 &
        !    !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,&
        !    !$OMP          lcol, lrow, l_a, l_b,                     &
        !    !$OMP          buffer, doubles, n_doubles,               &
        !    !$OMP          tmp_det2, idx, l, kcol_prev,              &
        !    !$OMP          singles_a, n_singles_a, singles_b,        &
-       !    !$OMP          n_singles_b, nkeys, keys, valus, c_average)
+       !    !$OMP          n_singles_b, k8)
    
    ! Alpha/Beta double excitations
    ! =============================
@@ -345,17 +353,9 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        enddo
        if(alpha_beta.or.spin_trace.or.alpha_alpha)then
        ! increment the alpha/beta part for single excitations
-       if (nkeys+ 2 * norb .ge. size(values)) then
+!!!!   call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
-         call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
-         nkeys = 0
-       endif
-       call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        ! increment the alpha/alpha part for single excitations
-       if (nkeys+4 * norb .ge. size(values)) then
+!!!!   call orb_range_off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
-         call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
-         nkeys = 0
-       endif
-       call orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        endif
        
      enddo
@@ -378,11 +378,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
           c_2(l) = u_t(l,k_a)
           c_average += c_1(l) * c_2(l) * state_weights(l)
         enddo
-        if (nkeys+4 .ge. size(values)) then
+!!!!    call orb_range_off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
-          nkeys = 0
-        endif
-        call orb_range_off_diag_double_to_two_rdm_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
       enddo
      endif
      
@@ -449,17 +445,9 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        enddo
        if(alpha_beta.or.spin_trace.or.beta_beta)then
         ! increment the alpha/beta  part for single excitations
-        if (nkeys+2 * norb .ge. size(values)) then
+!!!!    call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
-          nkeys = 0
-        endif
-        call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
         ! increment the beta /beta  part for single excitations
-        if (nkeys+4 * norb .ge. size(values)) then
+!!!!    call orb_range_off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
-          nkeys = 0
-        endif
-        call orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        endif
      enddo
      
@@ -481,11 +469,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
           c_2(l) = u_t(l,k_a)
           c_average += c_1(l) * c_2(l) * state_weights(l)
         enddo
-        if (nkeys+4 .ge. size(values)) then
+!!!!    call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
-          nkeys = 0
-        endif
-        call orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
         ASSERT (l_a <= N_det)
         
       enddo
@@ -525,7 +509,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
      
    end do
    !!$OMP END DO
-   deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
+   deallocate(buffer, singles_a, singles_b, doubles, idx)
    !!$OMP END PARALLEL
    
 end

src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f

@@ -13,7 +13,7 @@
  double precision, intent(in)   :: c_1
  integer                        :: occ(N_int*bit_kind_size,2)
  integer                        :: n_occ_ab(2)
- integer :: i,j,h1,h2
+ integer :: i,j,h1,h2,istate
  call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
  do i = 1, n_occ_ab(1)
   h1 = occ(i,1)
@@ -53,7 +53,7 @@
 
  integer                        :: occ(N_int*bit_kind_size,2)
  integer                        :: n_occ_ab(2)
- integer :: i,j,h1,h2
+ integer :: i,j,h1,h2,istate
  integer(bit_kind) :: det_1_act(N_int,2)
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
  do i = 1, N_int
@@ -193,7 +193,7 @@
  integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
  integer, intent(in) :: list_orb_reverse(mo_num)
  double precision, intent(in)   :: c_1
- integer :: i,j,h1,h2,p1,p2
+ integer :: i,j,h1,h2,p1,p2,istate
  integer                        :: exc(0:2,2,2)
  double precision               :: phase
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@@ -278,7 +278,7 @@
 
  integer                        :: occ(N_int*bit_kind_size,2)
  integer                        :: n_occ_ab(2)
- integer :: i,j,h1,h2,p1
+ integer :: i,j,h1,h2,istate,p1
  integer                        :: exc(0:2,2,2)
  double precision               :: phase
 
@@ -397,7 +397,7 @@
 
  integer                        :: occ(N_int*bit_kind_size,2)
  integer                        :: n_occ_ab(2)
- integer :: i,j,h1,h2,p1
+ integer :: i,j,h1,h2,istate,p1
  integer                        :: exc(0:2,2,2)
  double precision               :: phase
 
@@ -477,7 +477,7 @@
 
  integer                        :: occ(N_int*bit_kind_size,2)
  integer                        :: n_occ_ab(2)
- integer :: i,j,h1,h2,p1
+ integer :: i,j,h1,h2,istate,p1
  integer                        :: exc(0:2,2,2)
  double precision               :: phase
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@@ -510,15 +510,17 @@
    p1 =  exc(1,2,2)
    if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
    p1 = list_orb_reverse(p1)
-   do i = 1, n_occ_ab(2)
+   do istate = 1, N_states
-    h2 = occ(i,2)
+    do i = 1, n_occ_ab(2)
-    if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
+     h2 = occ(i,2)
-    h2 = list_orb_reverse(h2)
+     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
-    big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase
+     h2 = list_orb_reverse(h2)
-    big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase
+     big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase
+     big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase
  
-    big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase
+     big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase
-    big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase
+     big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase
+     enddo 
    enddo
   endif
  endif
@@ -555,7 +557,7 @@
  integer, intent(in) :: list_orb_reverse(mo_num)
  double precision, intent(in)   :: c_1
 
- integer :: i,j,h1,h2,p1,p2
+ integer :: i,j,h1,h2,p1,p2,istate
  integer                        :: exc(0:2,2)
  double precision               :: phase
 
@@ -588,11 +590,13 @@
  if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return
  p2 = list_orb_reverse(p2)
  if(alpha_alpha.or.spin_trace)then
+  do istate = 1, N_states
    big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase
    big_array(h1,h2,p2,p1) -= 0.5d0 * c_1 * phase
                                          
    big_array(h2,h1,p2,p1) += 0.5d0 * c_1 * phase
    big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase
+  enddo
  endif
  end
 
@@ -627,7 +631,7 @@
  integer, intent(in) :: list_orb_reverse(mo_num)
  double precision, intent(in)   :: c_1
 
- integer :: i,j,h1,h2,p1,p2
+ integer :: i,j,h1,h2,p1,p2,istate
  integer                        :: exc(0:2,2)
  double precision               :: phase
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace

src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f

@@ -26,7 +26,7 @@
 
  integer                        :: occ(N_int*bit_kind_size,2)
  integer                        :: n_occ_ab(2)
- integer :: i,j,h1,h2
+ integer :: i,j,h1,h2,istate
  integer(bit_kind) :: det_1_act(N_int,2)
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
  do i = 1, N_int
@@ -201,7 +201,7 @@
  double precision, intent(out)  :: values(sze_buff)
  integer         , intent(out)  :: keys(4,sze_buff)
  integer         , intent(inout):: nkeys
- integer :: i,j,h1,h2,p1,p2
+ integer :: i,j,h1,h2,p1,p2,istate
  integer                        :: exc(0:2,2,2)
  double precision               :: phase
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@@ -255,553 +255,429 @@
  endif
  end
 
-  subroutine orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
+! subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
-  use bitmasks
+! use bitmasks
-  BEGIN_DOC
+! BEGIN_DOC
- ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
- !
+!!
- ! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another
+!! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another
- ! 
+!! 
- ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
- ! 
+!! 
- ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
- !
+!!
- ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
- !
+!!
- ! ispin determines which spin-spin component of the two-rdm you will update 
+!! ispin determines which spin-spin component of the two-rdm you will update 
- !
+!!
- ! ispin   == 1 :: alpha/ alpha
+!! ispin   == 1 :: alpha/ alpha
- ! ispin   == 2 :: beta / beta
+!! ispin   == 2 :: beta / beta
- ! ispin   == 3 :: alpha/ beta
+!! ispin   == 3 :: alpha/ beta
- ! ispin   == 4 :: spin traced <=> total two-rdm 
+!! ispin   == 4 :: spin traced <=> total two-rdm 
- !
+!!
- ! here, only ispin == 3 or 4 will do something
+!! here, only ispin == 3 or 4 will do something
-  END_DOC
+! END_DOC
-  implicit none
+! implicit none
- integer, intent(in) :: ispin,sze_buff
+! integer, intent(in) :: dim1,ispin
- integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
+! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
- integer, intent(in) :: list_orb_reverse(mo_num)
+! integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
- double precision, intent(in)   :: c_1
+! integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
- double precision, intent(out)  :: values(sze_buff)
+! integer, intent(in) :: list_orb_reverse(mo_num)
- integer         , intent(out)  :: keys(4,sze_buff)
+! double precision, intent(in)   :: c_1
- integer         , intent(inout):: nkeys
+!
+! integer                        :: occ(N_int*bit_kind_size,2)
+! integer                        :: n_occ_ab(2)
+! integer :: i,j,h1,h2,istate,p1
+! integer                        :: exc(0:2,2,2)
+! double precision               :: phase
+!
+! logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
+! logical :: is_integer_in_string
+! alpha_alpha = .False.   
+! beta_beta   = .False.
+! alpha_beta  = .False.
+! spin_trace  = .False.
+! if(     ispin == 1)then
+!  alpha_alpha = .True.
+! else if(ispin == 2)then
+!  beta_beta   = .True.
+! else if(ispin == 3)then
+!  alpha_beta  = .True.
+! else if(ispin == 4)then
+!  spin_trace  = .True.
+! endif
+!
+! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
+! call get_single_excitation(det_1,det_2,exc,phase,N_int)
+! if(alpha_beta)then
+!  if (exc(0,1,1) == 1) then
+!   ! Mono alpha
+!   h1 = exc(1,1,1)
+!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+!   h1 = list_orb_reverse(h1)
+!   p1 = exc(1,2,1)
+!   if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+!   p1 = list_orb_reverse(p1)
+!    do i = 1, n_occ_ab(2)
+!     h2 = occ(i,2)
+!     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
+!     h2 = list_orb_reverse(h2)
+!     big_array(h1,h2,p1,h2) += c_1 * phase
+!    enddo 
+!  else 
+!   ! Mono beta
+!   h1 = exc(1,1,2)
+!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+!   h1 = list_orb_reverse(h1)
+!   p1 = exc(1,2,2)
+!   if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+!   p1 = list_orb_reverse(p1)
+!    do i = 1, n_occ_ab(1)
+!     h2 = occ(i,1)
+!     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
+!     h2 = list_orb_reverse(h2)
+!     big_array(h2,h1,h2,p1) += c_1 * phase
+!    enddo 
+!  endif
+! else if(spin_trace)then
+!  if (exc(0,1,1) == 1) then
+!   ! Mono alpha
+!   h1 = exc(1,1,1)
+!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+!   h1 = list_orb_reverse(h1)
+!   p1 = exc(1,2,1)
+!   if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+!   p1 = list_orb_reverse(p1)
+!    do i = 1, n_occ_ab(2)
+!     h2 = occ(i,2)
+!     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
+!     h2 = list_orb_reverse(h2)
+!     big_array(h1,h2,p1,h2) +=  0.5d0 * c_1 * phase
+!     big_array(h2,h1,h2,p1) +=  0.5d0 * c_1 * phase
+!    enddo 
+!  else 
+!   ! Mono beta
+!   h1 = exc(1,1,2)
+!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+!   h1 = list_orb_reverse(h1)
+!   p1 = exc(1,2,2)
+!   if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+!   p1 = list_orb_reverse(p1)
+!    do i = 1, n_occ_ab(1)
+!     h2 = occ(i,1)
+!     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
+!     h2 = list_orb_reverse(h2)
+!     big_array(h1,h2,p1,h2) +=  0.5d0 * c_1 * phase
+!     big_array(h2,h1,h2,p1) +=  0.5d0 * c_1 * phase
+!    enddo 
+!  endif
+! endif
+! end
 
-  integer                        :: occ(N_int*bit_kind_size,2)
+! subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
-  integer                        :: n_occ_ab(2)
+! BEGIN_DOC
-  integer :: i,j,h1,h2,p1
+!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
-  integer                        :: exc(0:2,2,2)
+!!
-  double precision               :: phase
+!! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another
+!! 
+!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+!! 
+!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!!
+!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!!
+!! ispin determines which spin-spin component of the two-rdm you will update 
+!!
+!! ispin   == 1 :: alpha/ alpha
+!! ispin   == 2 :: beta / beta
+!! ispin   == 3 :: alpha/ beta
+!! ispin   == 4 :: spin traced <=> total two-rdm 
+!!
+!! here, only ispin == 1 or 4 will do something
+! END_DOC
+! use bitmasks
+! implicit none
+! integer, intent(in) :: dim1,ispin
+! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
+! integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
+! integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
+! integer, intent(in) :: list_orb_reverse(mo_num)
+! double precision, intent(in)   :: c_1
+!
+! integer                        :: occ(N_int*bit_kind_size,2)
+! integer                        :: n_occ_ab(2)
+! integer :: i,j,h1,h2,istate,p1
+! integer                        :: exc(0:2,2,2)
+! double precision               :: phase
+!
+! logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
+! logical :: is_integer_in_string
+! alpha_alpha = .False.   
+! beta_beta   = .False.
+! alpha_beta  = .False.
+! spin_trace  = .False.
+! if(     ispin == 1)then
+!  alpha_alpha = .True.
+! else if(ispin == 2)then
+!  beta_beta   = .True.
+! else if(ispin == 3)then
+!  alpha_beta  = .True.
+! else if(ispin == 4)then
+!  spin_trace  = .True.
+! endif
+!
+! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
+! call get_single_excitation(det_1,det_2,exc,phase,N_int)
+! if(alpha_alpha.or.spin_trace)then
+!  if (exc(0,1,1) == 1) then
+!   ! Mono alpha
+!   h1 = exc(1,1,1)
+!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+!   h1 = list_orb_reverse(h1)
+!   p1 = exc(1,2,1)
+!   if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+!   p1 = list_orb_reverse(p1)
+!    do i = 1, n_occ_ab(1)
+!     h2 = occ(i,1)
+!     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
+!     h2 = list_orb_reverse(h2)
+!     big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase
+!     big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase
+! 
+!     big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase
+!     big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase
+!    enddo 
+!  else 
+!   return
+!  endif
+! endif
+! end
 
-  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
+! subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
-  logical :: is_integer_in_string
+! use bitmasks
-  alpha_alpha = .False.   
+! BEGIN_DOC
-  beta_beta   = .False.
+!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
-  alpha_beta  = .False.
+!!
-  spin_trace  = .False.
+!! a given couple of determinant det_1, det_2 being a BETA  SINGLE excitation with respect to one another
-  if(     ispin == 1)then
+!! 
-   alpha_alpha = .True.
+!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
-  else if(ispin == 2)then
+!! 
-   beta_beta   = .True.
+!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
-  else if(ispin == 3)then
+!!
-   alpha_beta  = .True.
+!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
-  else if(ispin == 4)then
+!!
-   spin_trace  = .True.
+!! ispin determines which spin-spin component of the two-rdm you will update 
-  endif
+!!
- 
+!! ispin   == 1 :: alpha/ alpha
-  call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
+!! ispin   == 2 :: beta / beta
-  call get_single_excitation(det_1,det_2,exc,phase,N_int)
+!! ispin   == 3 :: alpha/ beta
-  if(alpha_beta)then
+!! ispin   == 4 :: spin traced <=> total two-rdm 
-   if (exc(0,1,1) == 1) then
+!!
-    ! Mono alpha
+!! here, only ispin == 2 or 4 will do something
-    h1 = exc(1,1,1)
+! END_DOC
-    if(list_orb_reverse(h1).lt.0)return
+! implicit none
-    h1 = list_orb_reverse(h1)
+! integer, intent(in) :: dim1,ispin
-    p1 = exc(1,2,1)
+! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
-    if(list_orb_reverse(p1).lt.0)return
+! integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
-    p1 = list_orb_reverse(p1)
+! integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
-     do i = 1, n_occ_ab(2)
+! integer, intent(in) :: list_orb_reverse(mo_num)
-      h2 = occ(i,2)
+! double precision, intent(in)   :: c_1
-      if(list_orb_reverse(h2).lt.0)return
+!
-      h2 = list_orb_reverse(h2)
+!
-      nkeys += 1
+! integer                        :: occ(N_int*bit_kind_size,2)
-      values(nkeys) = c_1 * phase
+! integer                        :: n_occ_ab(2)
-      keys(1,nkeys) = h1
+! integer :: i,j,h1,h2,istate,p1
-      keys(2,nkeys) = h2
+! integer                        :: exc(0:2,2,2)
-      keys(3,nkeys) = p1
+! double precision               :: phase
-      keys(4,nkeys) = h2
+! logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
-     enddo 
+! logical :: is_integer_in_string
-   else 
+! alpha_alpha = .False.   
-    ! Mono beta
+! beta_beta   = .False.
-    h1 = exc(1,1,2)
+! alpha_beta  = .False.
-    if(list_orb_reverse(h1).lt.0)return
+! spin_trace  = .False.
-    h1 = list_orb_reverse(h1)
+! if(     ispin == 1)then
-    p1 = exc(1,2,2)
+!  alpha_alpha = .True.
-    if(list_orb_reverse(p1).lt.0)return
+! else if(ispin == 2)then
-    p1 = list_orb_reverse(p1)
+!  beta_beta   = .True.
-     do i = 1, n_occ_ab(1)
+! else if(ispin == 3)then
-      h2 = occ(i,1)
+!  alpha_beta  = .True.
-      if(list_orb_reverse(h2).lt.0)return
+! else if(ispin == 4)then
-      h2 = list_orb_reverse(h2)
+!  spin_trace  = .True.
-      nkeys += 1
+! endif
-      values(nkeys) = c_1 * phase
+!
-      keys(1,nkeys) = h1
+!
-      keys(2,nkeys) = h2
+! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
-      keys(3,nkeys) = p1
+! call get_single_excitation(det_1,det_2,exc,phase,N_int)
-      keys(4,nkeys) = h2
+! if(beta_beta.or.spin_trace)then
-     enddo 
+!  if (exc(0,1,1) == 1) then
-   endif
+!   return
-  else if(spin_trace)then
+!  else
-   if (exc(0,1,1) == 1) then
+!   ! Mono beta
-    ! Mono alpha
+!   h1 =  exc(1,1,2)
-    h1 = exc(1,1,1)
+!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
-    if(list_orb_reverse(h1).lt.0)return
+!   h1 = list_orb_reverse(h1)
-    h1 = list_orb_reverse(h1)
+!   p1 =  exc(1,2,2)
-    p1 = exc(1,2,1)
+!   if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
-    if(list_orb_reverse(p1).lt.0)return
+!   p1 = list_orb_reverse(p1)
-    p1 = list_orb_reverse(p1)
+!   do istate = 1, N_states
-     do i = 1, n_occ_ab(2)
+!    do i = 1, n_occ_ab(2)
-      h2 = occ(i,2)
+!     h2 = occ(i,2)
-      if(list_orb_reverse(h2).lt.0)return
+!     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
-      h2 = list_orb_reverse(h2)
+!     h2 = list_orb_reverse(h2)
-      nkeys += 1
+!     big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase
-      values(nkeys) = 0.5d0 * c_1 * phase
+!     big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase
-      keys(1,nkeys) = h1
+! 
-      keys(2,nkeys) = h2
+!     big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase
-      keys(3,nkeys) = p1
+!     big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase
-      keys(4,nkeys) = h2
+!     enddo 
-      nkeys += 1
+!   enddo
-      values(nkeys) = 0.5d0 * c_1 * phase
+!  endif
-      keys(1,nkeys) = h2
+! endif
-      keys(2,nkeys) = h1
+! end
-      keys(3,nkeys) = h2
-      keys(4,nkeys) = p1
-     enddo 
-   else 
-    ! Mono beta
-    h1 = exc(1,1,2)
-    if(list_orb_reverse(h1).lt.0)return
-    h1 = list_orb_reverse(h1)
-    p1 = exc(1,2,2)
-    if(list_orb_reverse(p1).lt.0)return
-    p1 = list_orb_reverse(p1)
-   !print*,'****************'
-   !print*,'****************'
-   !print*,'h1,p1',h1,p1
-     do i = 1, n_occ_ab(1)
-      h2 = occ(i,1)
-      if(list_orb_reverse(h2).lt.0)return
-      h2 = list_orb_reverse(h2)
-   !  print*,'h2 = ',h2
-      nkeys += 1
-      values(nkeys) = 0.5d0 * c_1 * phase
-      keys(1,nkeys) = h1
-      keys(2,nkeys) = h2
-      keys(3,nkeys) = p1
-      keys(4,nkeys) = h2
-      nkeys += 1
-      values(nkeys) = 0.5d0 * c_1 * phase
-      keys(1,nkeys) = h2
-      keys(2,nkeys) = h1
-      keys(3,nkeys) = h2
-      keys(4,nkeys) = p1
-     enddo 
-   endif
-  endif
-  end
-
-  subroutine orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-  BEGIN_DOC
- ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
- !
- ! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another
- ! 
- ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
- ! 
- ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
- !
- ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
- !
- ! ispin determines which spin-spin component of the two-rdm you will update 
- !
- ! ispin   == 1 :: alpha/ alpha
- ! ispin   == 2 :: beta / beta
- ! ispin   == 3 :: alpha/ beta
- ! ispin   == 4 :: spin traced <=> total two-rdm 
- !
- ! here, only ispin == 1 or 4 will do something
-  END_DOC
-  use bitmasks
-  implicit none
- integer, intent(in) :: ispin,sze_buff
- integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
- integer, intent(in) :: list_orb_reverse(mo_num)
- double precision, intent(in)   :: c_1
- double precision, intent(out)  :: values(sze_buff)
- integer         , intent(out)  :: keys(4,sze_buff)
- integer         , intent(inout):: nkeys
- 
-  integer                        :: occ(N_int*bit_kind_size,2)
-  integer                        :: n_occ_ab(2)
-  integer :: i,j,h1,h2,p1
-  integer                        :: exc(0:2,2,2)
-  double precision               :: phase
- 
-  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
-  logical :: is_integer_in_string
-  alpha_alpha = .False.   
-  beta_beta   = .False.
-  alpha_beta  = .False.
-  spin_trace  = .False.
-  if(     ispin == 1)then
-   alpha_alpha = .True.
-  else if(ispin == 2)then
-   beta_beta   = .True.
-  else if(ispin == 3)then
-   alpha_beta  = .True.
-  else if(ispin == 4)then
-   spin_trace  = .True.
-  endif
- 
-  call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
-  call get_single_excitation(det_1,det_2,exc,phase,N_int)
-  if(alpha_alpha.or.spin_trace)then
-   if (exc(0,1,1) == 1) then
-    ! Mono alpha
-    h1 = exc(1,1,1)
-    if(list_orb_reverse(h1).lt.0)return
-    h1 = list_orb_reverse(h1)
-    p1 = exc(1,2,1)
-    if(list_orb_reverse(p1).lt.0)return
-    p1 = list_orb_reverse(p1)
-     do i = 1, n_occ_ab(1)
-      h2 = occ(i,1)
-      if(list_orb_reverse(h2).lt.0)return
-      h2 = list_orb_reverse(h2)
-
-      nkeys += 1
-      values(nkeys) = 0.5d0 * c_1 * phase
-      keys(1,nkeys) = h1
-      keys(2,nkeys) = h2
-      keys(3,nkeys) = p1
-      keys(4,nkeys) = h2
-
-      nkeys += 1
-      values(nkeys) = - 0.5d0 * c_1 * phase
-      keys(1,nkeys) = h1
-      keys(2,nkeys) = h2
-      keys(3,nkeys) = h2
-      keys(4,nkeys) = p1
-  
-      nkeys += 1
-      values(nkeys) = 0.5d0 * c_1 * phase
-      keys(1,nkeys) = h2
-      keys(2,nkeys) = h1
-      keys(3,nkeys) = h2
-      keys(4,nkeys) = p1
-
-      nkeys += 1
-      values(nkeys) = - 0.5d0 * c_1 * phase
-      keys(1,nkeys) = h2
-      keys(2,nkeys) = h1
-      keys(3,nkeys) = p1
-      keys(4,nkeys) = h2
-     enddo 
-   else 
-    return
-   endif
-  endif
-  end
-
-  subroutine orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-  use bitmasks
-  BEGIN_DOC
- ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
- !
- ! a given couple of determinant det_1, det_2 being a BETA  SINGLE excitation with respect to one another
- ! 
- ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
- ! 
- ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
- !
- ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
- !
- ! ispin determines which spin-spin component of the two-rdm you will update 
- !
- ! ispin   == 1 :: alpha/ alpha
- ! ispin   == 2 :: beta / beta
- ! ispin   == 3 :: alpha/ beta
- ! ispin   == 4 :: spin traced <=> total two-rdm 
- !
- ! here, only ispin == 2 or 4 will do something
-  END_DOC
-  implicit none
- integer, intent(in) :: ispin,sze_buff
- integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
- integer, intent(in) :: list_orb_reverse(mo_num)
- double precision, intent(in)   :: c_1
- double precision, intent(out)  :: values(sze_buff)
- integer         , intent(out)  :: keys(4,sze_buff)
- integer         , intent(inout):: nkeys
- 
-  integer                        :: occ(N_int*bit_kind_size,2)
-  integer                        :: n_occ_ab(2)
-  integer :: i,j,h1,h2,p1
-  integer                        :: exc(0:2,2,2)
-  double precision               :: phase
-  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
-  logical :: is_integer_in_string
-  alpha_alpha = .False.   
-  beta_beta   = .False.
-  alpha_beta  = .False.
-  spin_trace  = .False.
-  if(     ispin == 1)then
-   alpha_alpha = .True.
-  else if(ispin == 2)then
-   beta_beta   = .True.
-  else if(ispin == 3)then
-   alpha_beta  = .True.
-  else if(ispin == 4)then
-   spin_trace  = .True.
-  endif
 
 
-  call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
+! subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
-  call get_single_excitation(det_1,det_2,exc,phase,N_int)
+! use bitmasks
-  if(beta_beta.or.spin_trace)then
+! BEGIN_DOC
-   if (exc(0,1,1) == 1) then
+!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
-    return
+!!
-   else
+!! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another
-    ! Mono beta
+!! 
-    h1 =  exc(1,1,2)
+!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
-    if(list_orb_reverse(h1).lt.0)return
+!! 
-    h1 = list_orb_reverse(h1)
+!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
-    p1 =  exc(1,2,2)
+!!
-    if(list_orb_reverse(p1).lt.0)return
+!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
-    p1 = list_orb_reverse(p1)
+!!
-    do i = 1, n_occ_ab(2)
+!! ispin determines which spin-spin component of the two-rdm you will update 
-     h2 = occ(i,2)
+!!
-     if(list_orb_reverse(h2).lt.0)return
+!! ispin   == 1 :: alpha/ alpha
-     h2 = list_orb_reverse(h2)
+!! ispin   == 2 :: beta / beta
-     nkeys += 1
+!! ispin   == 3 :: alpha/ beta
-     values(nkeys) = 0.5d0 * c_1 * phase
+!! ispin   == 4 :: spin traced <=> total two-rdm 
-     keys(1,nkeys) = h1
+!!
-     keys(2,nkeys) = h2
+!! here, only ispin == 1 or 4 will do something
-     keys(3,nkeys) = p1
+! END_DOC
-     keys(4,nkeys) = h2
+! implicit none
+! integer, intent(in) :: dim1,ispin
+! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
+! integer(bit_kind), intent(in)  :: det_1(N_int),det_2(N_int)
+! integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
+! integer, intent(in) :: list_orb_reverse(mo_num)
+! double precision, intent(in)   :: c_1
+!
+! integer :: i,j,h1,h2,p1,p2,istate
+! integer                        :: exc(0:2,2)
+! double precision               :: phase
+!
+! logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
+! logical :: is_integer_in_string
+! alpha_alpha = .False.   
+! beta_beta   = .False.
+! alpha_beta  = .False.
+! spin_trace  = .False.
+! if(     ispin == 1)then
+!  alpha_alpha = .True.
+! else if(ispin == 2)then
+!  beta_beta   = .True.
+! else if(ispin == 3)then
+!  alpha_beta  = .True.
+! else if(ispin == 4)then
+!  spin_trace  = .True.
+! endif
+! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
+! h1 =exc(1,1)
+! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+! h1 = list_orb_reverse(h1)
+! h2 =exc(2,1)
+! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return
+! h2 = list_orb_reverse(h2)
+! p1 =exc(1,2)
+! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+! p1 = list_orb_reverse(p1)
+! p2 =exc(2,2)
+! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return
+! p2 = list_orb_reverse(p2)
+! if(alpha_alpha.or.spin_trace)then
+!  do istate = 1, N_states
+!   big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase
+!   big_array(h1,h2,p2,p1) -= 0.5d0 * c_1 * phase
+!                                         
+!   big_array(h2,h1,p2,p1) += 0.5d0 * c_1 * phase
+!   big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase
+!  enddo
+! endif
+! end
 
-     nkeys += 1
+! subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
-     values(nkeys) = - 0.5d0 * c_1 * phase
+! use bitmasks
-     keys(1,nkeys) = h1
+! BEGIN_DOC
-     keys(2,nkeys) = h2
+!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
-     keys(3,nkeys) = h2
+!!
-     keys(4,nkeys) = p1
+!! a given couple of determinant det_1, det_2 being a BETA /BETA  DOUBLE excitation with respect to one another
-  
+!! 
-     nkeys += 1
+!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
-     values(nkeys) = 0.5d0 * c_1 * phase
+!! 
-     keys(1,nkeys) = h2
+!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
-     keys(2,nkeys) = h1
+!!
-     keys(3,nkeys) = h2
+!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
-     keys(4,nkeys) = p1
+!!
-
+!! ispin determines which spin-spin component of the two-rdm you will update 
-     nkeys += 1
+!!
-     values(nkeys) = - 0.5d0 * c_1 * phase
+!! ispin   == 1 :: alpha/ alpha
-     keys(1,nkeys) = h2
+!! ispin   == 2 :: beta / beta
-     keys(2,nkeys) = h1
+!! ispin   == 3 :: alpha/ beta
-     keys(3,nkeys) = p1
+!! ispin   == 4 :: spin traced <=> total two-rdm 
-     keys(4,nkeys) = h2
+!!
-    enddo 
+!! here, only ispin == 2 or 4 will do something
-   endif
+! END_DOC
-  endif
+! implicit none
-  end
+!
-
+! integer, intent(in) :: dim1,ispin
-
+! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
-  subroutine orb_range_off_diag_double_to_two_rdm_aa_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
+! integer(bit_kind), intent(in)  :: det_1(N_int),det_2(N_int)
-  use bitmasks
+! integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
-  BEGIN_DOC
+! integer, intent(in) :: list_orb_reverse(mo_num)
- ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+! double precision, intent(in)   :: c_1
- !
+!
- ! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another
+! integer :: i,j,h1,h2,p1,p2,istate
- ! 
+! integer                        :: exc(0:2,2)
- ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! double precision               :: phase
- ! 
+! logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
- ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+! logical :: is_integer_in_string
- !
+! alpha_alpha = .False.   
- ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+! beta_beta   = .False.
- !
+! alpha_beta  = .False.
- ! ispin determines which spin-spin component of the two-rdm you will update 
+! spin_trace  = .False.
- !
+! if(     ispin == 1)then
- ! ispin   == 1 :: alpha/ alpha
+!  alpha_alpha = .True.
- ! ispin   == 2 :: beta / beta
+! else if(ispin == 2)then
- ! ispin   == 3 :: alpha/ beta
+!  beta_beta   = .True.
- ! ispin   == 4 :: spin traced <=> total two-rdm 
+! else if(ispin == 3)then
- !
+!  alpha_beta  = .True.
- ! here, only ispin == 1 or 4 will do something
+! else if(ispin == 4)then
-  END_DOC
+!  spin_trace  = .True.
-  implicit none
+! endif
- integer, intent(in) :: ispin,sze_buff
+!
- integer(bit_kind), intent(in)  :: det_1(N_int),det_2(N_int)
+! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
- integer, intent(in) :: list_orb_reverse(mo_num)
+! h1 =exc(1,1)
- double precision, intent(in)   :: c_1
+! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
- double precision, intent(out)  :: values(sze_buff)
+! h1 = list_orb_reverse(h1)
- integer         , intent(out)  :: keys(4,sze_buff)
+! h2 =exc(2,1)
- integer         , intent(inout):: nkeys
+! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return
- 
+! h2 = list_orb_reverse(h2)
- 
+! p1 =exc(1,2)
-  integer :: i,j,h1,h2,p1,p2
+! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
-  integer                        :: exc(0:2,2)
+! p1 = list_orb_reverse(p1)
-  double precision               :: phase
+! p2 =exc(2,2)
- 
+! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return
-  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
+! p2 = list_orb_reverse(p2)
-  logical :: is_integer_in_string
+! if(beta_beta.or.spin_trace)then
-  alpha_alpha = .False.   
+!  big_array(h1,h2,p1,p2) += 0.5d0 * c_1* phase 
-  beta_beta   = .False.
+!  big_array(h1,h2,p2,p1) -= 0.5d0 * c_1* phase 
-  alpha_beta  = .False.
+!
-  spin_trace  = .False.
+!  big_array(h2,h1,p2,p1) += 0.5d0 * c_1* phase 
-  if(     ispin == 1)then
+!  big_array(h2,h1,p1,p2) -= 0.5d0 * c_1* phase 
-   alpha_alpha = .True.
+! endif
-  else if(ispin == 2)then
+! end
-   beta_beta   = .True.
-  else if(ispin == 3)then
-   alpha_beta  = .True.
-  else if(ispin == 4)then
-   spin_trace  = .True.
-  endif
-  call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
-  h1 =exc(1,1)
-  if(list_orb_reverse(h1).lt.0)return
-  h1 = list_orb_reverse(h1)
-  h2 =exc(2,1)
-  if(list_orb_reverse(h2).lt.0)return
-  h2 = list_orb_reverse(h2)
-  p1 =exc(1,2)
-  if(list_orb_reverse(p1).lt.0)return
-  p1 = list_orb_reverse(p1)
-  p2 =exc(2,2)
-  if(list_orb_reverse(p2).lt.0)return
-  p2 = list_orb_reverse(p2)
-  if(alpha_alpha.or.spin_trace)then
-    nkeys += 1
-    values(nkeys) = 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h1
-    keys(2,nkeys) = h2
-    keys(3,nkeys) = p1
-    keys(4,nkeys) = p2
-
-    nkeys += 1
-    values(nkeys) = - 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h1
-    keys(2,nkeys) = h2
-    keys(3,nkeys) = p2
-    keys(4,nkeys) = p1
-                                          
-    nkeys += 1
-    values(nkeys) = 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h2
-    keys(2,nkeys) = h1
-    keys(3,nkeys) = p2
-    keys(4,nkeys) = p1
-
-    nkeys += 1
-    values(nkeys) = - 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h2
-    keys(2,nkeys) = h1
-    keys(3,nkeys) = p1
-    keys(4,nkeys) = p2
-  endif
-  end
-
-  subroutine orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-  use bitmasks
-  BEGIN_DOC
- ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
- !
- ! a given couple of determinant det_1, det_2 being a BETA /BETA  DOUBLE excitation with respect to one another
- ! 
- ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
- ! 
- ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
- !
- ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
- !
- ! ispin determines which spin-spin component of the two-rdm you will update 
- !
- ! ispin   == 1 :: alpha/ alpha
- ! ispin   == 2 :: beta / beta
- ! ispin   == 3 :: alpha/ beta
- ! ispin   == 4 :: spin traced <=> total two-rdm 
- !
- ! here, only ispin == 2 or 4 will do something
-  END_DOC
-  implicit none
- 
- integer, intent(in) :: ispin,sze_buff
- integer(bit_kind), intent(in)  :: det_1(N_int),det_2(N_int)
- integer, intent(in) :: list_orb_reverse(mo_num)
- double precision, intent(in)   :: c_1
- double precision, intent(out)  :: values(sze_buff)
- integer         , intent(out)  :: keys(4,sze_buff)
- integer         , intent(inout):: nkeys
- 
-  integer :: i,j,h1,h2,p1,p2
-  integer                        :: exc(0:2,2)
-  double precision               :: phase
-  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
-  logical :: is_integer_in_string
-  alpha_alpha = .False.   
-  beta_beta   = .False.
-  alpha_beta  = .False.
-  spin_trace  = .False.
-  if(     ispin == 1)then
-   alpha_alpha = .True.
-  else if(ispin == 2)then
-   beta_beta   = .True.
-  else if(ispin == 3)then
-   alpha_beta  = .True.
-  else if(ispin == 4)then
-   spin_trace  = .True.
-  endif
- 
-  call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
-  h1 =exc(1,1)
-  if(list_orb_reverse(h1).lt.0)return
-  h1 = list_orb_reverse(h1)
-  h2 =exc(2,1)
-  if(list_orb_reverse(h2).lt.0)return
-  h2 = list_orb_reverse(h2)
-  p1 =exc(1,2)
-  if(list_orb_reverse(p1).lt.0)return
-  p1 = list_orb_reverse(p1)
-  p2 =exc(2,2)
-  if(list_orb_reverse(p2).lt.0)return
-  p2 = list_orb_reverse(p2)
-  if(beta_beta.or.spin_trace)then
-    nkeys += 1
-    values(nkeys) = 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h1
-    keys(2,nkeys) = h2
-    keys(3,nkeys) = p1
-    keys(4,nkeys) = p2
-
-    nkeys += 1
-    values(nkeys) = - 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h1
-    keys(2,nkeys) = h2
-    keys(3,nkeys) = p2
-    keys(4,nkeys) = p1
-                                          
-    nkeys += 1
-    values(nkeys) = 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h2
-    keys(2,nkeys) = h1
-    keys(3,nkeys) = p2
-    keys(4,nkeys) = p1
-
-    nkeys += 1
-    values(nkeys) = - 0.5d0 * c_1 * phase
-    keys(1,nkeys) = h2
-    keys(2,nkeys) = h1
-    keys(3,nkeys) = p1
-    keys(4,nkeys) = p2
-  endif
-  end