bug fixed in openmp 2 rdms

src/bitmask/core_inact_act_virt.irp.f

@@ -322,6 +322,7 @@ 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_alpha_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_beta_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_alpha_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_beta_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,6 +31,7 @@ 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)
    
@@ -135,6 +136,7 @@ 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)
@@ -443,7 +445,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_alpha_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_beta_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,11 +93,9 @@ 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, m_a, m_b
+   integer                        :: k_a, k_b, l_a, l_b
-   integer                        :: istate
+   integer                        :: krow, kcol
-   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)
@@ -109,7 +107,6 @@ 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
@@ -136,16 +133,11 @@ 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
+   sze_buff = norb ** 3 + 6 * norb
    list_orb_reverse = -1000 
    do i = 1, norb
     list_orb_reverse(list_orb(i)) = i 
@@ -173,13 +165,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)                 &
+       !    !$OMP          ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin)     &
        !    !$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, k8)
+       !    !$OMP          n_singles_b, nkeys, keys, valus, c_average)
    
    ! Alpha/Beta double excitations
    ! =============================
@@ -353,9 +345,17 @@ 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
-!!!!   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)
+       if (nkeys+ 2 * norb .ge. size(values)) then
+         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
-!!!!   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)
+       if (nkeys+4 * norb .ge. size(values)) then
+         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,7 +378,11 @@ 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
-!!!!    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)
+        if (nkeys+4 .ge. size(values)) then
+          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
      
@@ -445,9 +449,17 @@ 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
-!!!!    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)
+        if (nkeys+2 * norb .ge. size(values)) then
+          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
-!!!!    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)
+        if (nkeys+4 * norb .ge. size(values)) then
+          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
      
@@ -469,7 +481,11 @@ 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
-!!!!    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)
+        if (nkeys+4 .ge. size(values)) then
+          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
@@ -509,7 +525,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)
+   deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
    !!$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,istate
+ integer :: i,j,h1,h2
  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,istate
+ integer :: i,j,h1,h2
  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,istate
+ integer :: i,j,h1,h2,p1,p2
  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,istate,p1
+ integer :: i,j,h1,h2,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,istate,p1
+ integer :: i,j,h1,h2,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,istate,p1
+ integer :: i,j,h1,h2,p1
  integer                        :: exc(0:2,2,2)
  double precision               :: phase
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@@ -510,7 +510,6 @@
    p1 =  exc(1,2,2)
    if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
    p1 = list_orb_reverse(p1)
-   do istate = 1, N_states
    do i = 1, n_occ_ab(2)
     h2 = occ(i,2)
     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
@@ -521,7 +520,6 @@
     big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase
     big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase
    enddo 
-   enddo
   endif
  endif
  end
@@ -557,7 +555,7 @@
  integer, intent(in) :: list_orb_reverse(mo_num)
  double precision, intent(in)   :: c_1
 
- integer :: i,j,h1,h2,p1,p2,istate
+ integer :: i,j,h1,h2,p1,p2
  integer                        :: exc(0:2,2)
  double precision               :: phase
 
@@ -590,13 +588,11 @@
  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
 
@@ -631,7 +627,7 @@
  integer, intent(in) :: list_orb_reverse(mo_num)
  double precision, intent(in)   :: c_1
 
- integer :: i,j,h1,h2,p1,p2,istate
+ integer :: i,j,h1,h2,p1,p2
  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,istate
+ integer :: i,j,h1,h2
  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,istate
+ integer :: i,j,h1,h2,p1,p2
  integer                        :: exc(0:2,2,2)
  double precision               :: phase
  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@@ -255,429 +255,553 @@
  endif
  end
 
-! subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
+  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)
-! 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
-!! 
-!! 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 == 3 or 4 will do something
-! END_DOC
-! 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)
+ ! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another
-! 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
+ ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
-! 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)
+ ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
-! call get_single_excitation(det_1,det_2,exc,phase,N_int)
+ !
-! if(alpha_beta)then
+ ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
-!  if (exc(0,1,1) == 1) then
+ !
-!   ! Mono alpha
+ ! ispin determines which spin-spin component of the two-rdm you will update 
-!   h1 = exc(1,1,1)
+ !
-!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+ ! ispin   == 1 :: alpha/ alpha
-!   h1 = list_orb_reverse(h1)
+ ! ispin   == 2 :: beta / beta
-!   p1 = exc(1,2,1)
+ ! ispin   == 3 :: alpha/ beta
-!   if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+ ! ispin   == 4 :: spin traced <=> total two-rdm 
-!   p1 = list_orb_reverse(p1)
+ !
-!    do i = 1, n_occ_ab(2)
+ ! here, only ispin == 3 or 4 will do something
-!     h2 = occ(i,2)
+  END_DOC
-!     if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
+  implicit none
-!     h2 = list_orb_reverse(h2)
+ integer, intent(in) :: ispin,sze_buff
-!     big_array(h1,h2,p1,h2) += c_1 * phase
+ integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
-!    enddo 
+ integer, intent(in) :: list_orb_reverse(mo_num)
-!  else 
+ double precision, intent(in)   :: c_1
-!   ! Mono beta
+ double precision, intent(out)  :: values(sze_buff)
-!   h1 = exc(1,1,2)
+ integer         , intent(out)  :: keys(4,sze_buff)
-!   if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+ integer         , intent(inout):: nkeys
-!   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
  
-! subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
+  integer                        :: occ(N_int*bit_kind_size,2)
-! BEGIN_DOC
+  integer                        :: n_occ_ab(2)
-!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+  integer :: i,j,h1,h2,p1
-!!
+  integer                        :: exc(0:2,2,2)
-!! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another
+  double precision               :: phase
-!! 
-!! 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
  
-! subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
+  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
-! use bitmasks
+  logical :: is_integer_in_string
-! BEGIN_DOC
+  alpha_alpha = .False.   
-!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+  beta_beta   = .False.
-!!
+  alpha_beta  = .False.
-!! a given couple of determinant det_1, det_2 being a BETA  SINGLE excitation with respect to one another
+  spin_trace  = .False.
-!! 
+  if(     ispin == 1)then
-!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+   alpha_alpha = .True.
-!! 
+  else if(ispin == 2)then
-!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+   beta_beta   = .True.
-!!
+  else if(ispin == 3)then
-!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+   alpha_beta  = .True.
-!!
+  else if(ispin == 4)then
-!! ispin determines which spin-spin component of the two-rdm you will update 
+   spin_trace  = .True.
-!!
+  endif
-!! ispin   == 1 :: alpha/ alpha
+ 
-!! ispin   == 2 :: beta / beta
+  call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
-!! ispin   == 3 :: alpha/ beta
+  call get_single_excitation(det_1,det_2,exc,phase,N_int)
-!! ispin   == 4 :: spin traced <=> total two-rdm 
+  if(alpha_beta)then
-!!
+   if (exc(0,1,1) == 1) then
-!! here, only ispin == 2 or 4 will do something
+    ! Mono alpha
-! END_DOC
+    h1 = exc(1,1,1)
-! implicit none
+    if(list_orb_reverse(h1).lt.0)return
-! integer, intent(in) :: dim1,ispin
+    h1 = list_orb_reverse(h1)
-! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
+    p1 = exc(1,2,1)
-! integer(bit_kind), intent(in)  :: det_1(N_int,2),det_2(N_int,2)
+    if(list_orb_reverse(p1).lt.0)return
-! integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
+    p1 = list_orb_reverse(p1)
-! integer, intent(in) :: list_orb_reverse(mo_num)
+     do i = 1, n_occ_ab(2)
-! double precision, intent(in)   :: c_1
+      h2 = occ(i,2)
+      if(list_orb_reverse(h2).lt.0)return
+      h2 = list_orb_reverse(h2)
+      nkeys += 1
+      values(nkeys) = c_1 * phase
+      keys(1,nkeys) = h1
+      keys(2,nkeys) = h2
+      keys(3,nkeys) = p1
+      keys(4,nkeys) = h2
+     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)
+     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) = c_1 * phase
+      keys(1,nkeys) = h1
+      keys(2,nkeys) = h2
+      keys(3,nkeys) = p1
+      keys(4,nkeys) = h2
+     enddo 
+   endif
+  else if(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(2)
+      h2 = occ(i,2)
+      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) = h2
+      keys(2,nkeys) = h1
+      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
  ! 
-! integer                        :: occ(N_int*bit_kind_size,2)
+ ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
-! 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
  ! 
+ ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
  !
-! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
+ ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
-! call get_single_excitation(det_1,det_2,exc,phase,N_int)
-! if(beta_beta.or.spin_trace)then
-!  if (exc(0,1,1) == 1) then
-!   return
-!  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 istate = 1, N_states
-!    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(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase
  !
-!     big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase
+ ! ispin determines which spin-spin component of the two-rdm you will update 
-!     big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase
+ !
-!     enddo 
+ ! ispin   == 1 :: alpha/ alpha
-!   enddo
+ ! ispin   == 2 :: beta / beta
-!  endif
+ ! ispin   == 3 :: alpha/ beta
-! endif
+ ! ispin   == 4 :: spin traced <=> total two-rdm 
-! end
+ !
+ ! 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
  
  
-! subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
+  call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
-! use bitmasks
+  call get_single_excitation(det_1,det_2,exc,phase,N_int)
-! BEGIN_DOC
+  if(beta_beta.or.spin_trace)then
-!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+   if (exc(0,1,1) == 1) then
-!!
+    return
-!! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another
+   else
-!! 
+    ! Mono beta
-!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+    h1 =  exc(1,1,2)
-!! 
+    if(list_orb_reverse(h1).lt.0)return
-!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+    h1 = list_orb_reverse(h1)
-!!
+    p1 =  exc(1,2,2)
-!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+    if(list_orb_reverse(p1).lt.0)return
-!!
+    p1 = list_orb_reverse(p1)
-!! ispin determines which spin-spin component of the two-rdm you will update 
+    do i = 1, n_occ_ab(2)
-!!
+     h2 = occ(i,2)
-!! ispin   == 1 :: alpha/ alpha
+     if(list_orb_reverse(h2).lt.0)return
-!! ispin   == 2 :: beta / beta
+     h2 = list_orb_reverse(h2)
-!! ispin   == 3 :: alpha/ beta
+     nkeys += 1
-!! ispin   == 4 :: spin traced <=> total two-rdm 
+     values(nkeys) = 0.5d0 * c_1 * phase
-!!
+     keys(1,nkeys) = h1
-!! here, only ispin == 1 or 4 will do something
+     keys(2,nkeys) = h2
-! END_DOC
+     keys(3,nkeys) = p1
-! implicit none
+     keys(4,nkeys) = h2
-! 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
 
-! subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin)
+     nkeys += 1
-! use bitmasks
+     values(nkeys) = - 0.5d0 * c_1 * phase
-! BEGIN_DOC
+     keys(1,nkeys) = h1
-!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+     keys(2,nkeys) = h2
-!!
+     keys(3,nkeys) = h2
-!! a given couple of determinant det_1, det_2 being a BETA /BETA  DOUBLE excitation with respect to one another
+     keys(4,nkeys) = p1
-!! 
+  
-!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+     nkeys += 1
-!! 
+     values(nkeys) = 0.5d0 * c_1 * phase
-!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+     keys(1,nkeys) = h2
-!!
+     keys(2,nkeys) = h1
-!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+     keys(3,nkeys) = h2
-!!
+     keys(4,nkeys) = p1
-!! ispin determines which spin-spin component of the two-rdm you will update 
+
-!!
+     nkeys += 1
-!! ispin   == 1 :: alpha/ alpha
+     values(nkeys) = - 0.5d0 * c_1 * phase
-!! ispin   == 2 :: beta / beta
+     keys(1,nkeys) = h2
-!! ispin   == 3 :: alpha/ beta
+     keys(2,nkeys) = h1
-!! ispin   == 4 :: spin traced <=> total two-rdm 
+     keys(3,nkeys) = p1
-!!
+     keys(4,nkeys) = h2
-!! here, only ispin == 2 or 4 will do something
+    enddo 
-! END_DOC
+   endif
-! implicit none
+  endif
-!
+  end
-! 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)
+  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)  :: orb_bitmask(N_int)
+  use bitmasks
-! integer, intent(in) :: list_orb_reverse(mo_num)
+  BEGIN_DOC
-! double precision, intent(in)   :: c_1
+ ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
-!
+ !
-! integer :: i,j,h1,h2,p1,p2,istate
+ ! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another
-! integer                        :: exc(0:2,2)
+ ! 
-! double precision               :: phase
+ ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
-! logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
+ ! 
-! logical :: is_integer_in_string
+ ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
-! alpha_alpha = .False.   
+ !
-! beta_beta   = .False.
+ ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
-! alpha_beta  = .False.
+ !
-! spin_trace  = .False.
+ ! ispin determines which spin-spin component of the two-rdm you will update 
-! if(     ispin == 1)then
+ !
-!  alpha_alpha = .True.
+ ! ispin   == 1 :: alpha/ alpha
-! else if(ispin == 2)then
+ ! ispin   == 2 :: beta / beta
-!  beta_beta   = .True.
+ ! ispin   == 3 :: alpha/ beta
-! else if(ispin == 3)then
+ ! ispin   == 4 :: spin traced <=> total two-rdm 
-!  alpha_beta  = .True.
+ !
-! else if(ispin == 4)then
+ ! here, only ispin == 1 or 4 will do something
-!  spin_trace  = .True.
+  END_DOC
-! endif
+  implicit none
-!
+ integer, intent(in) :: ispin,sze_buff
-! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
+ integer(bit_kind), intent(in)  :: det_1(N_int),det_2(N_int)
-! h1 =exc(1,1)
+ integer, intent(in) :: list_orb_reverse(mo_num)
-! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
+ double precision, intent(in)   :: c_1
-! h1 = list_orb_reverse(h1)
+ double precision, intent(out)  :: values(sze_buff)
-! h2 =exc(2,1)
+ integer         , intent(out)  :: keys(4,sze_buff)
-! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return
+ integer         , intent(inout):: nkeys
-! h2 = list_orb_reverse(h2)
+ 
-! p1 =exc(1,2)
+ 
-! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
+  integer :: i,j,h1,h2,p1,p2
-! p1 = list_orb_reverse(p1)
+  integer                        :: exc(0:2,2)
-! p2 =exc(2,2)
+  double precision               :: phase
-! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return
+ 
-! p2 = list_orb_reverse(p2)
+  logical                        :: alpha_alpha,beta_beta,alpha_beta,spin_trace
-! if(beta_beta.or.spin_trace)then
+  logical :: is_integer_in_string
-!  big_array(h1,h2,p1,p2) += 0.5d0 * c_1* phase 
+  alpha_alpha = .False.   
-!  big_array(h1,h2,p2,p1) -= 0.5d0 * c_1* phase 
+  beta_beta   = .False.
-!
+  alpha_beta  = .False.
-!  big_array(h2,h1,p2,p1) += 0.5d0 * c_1* phase 
+  spin_trace  = .False.
-!  big_array(h2,h1,p1,p2) -= 0.5d0 * c_1* phase 
+  if(     ispin == 1)then
-! endif
+   alpha_alpha = .True.
-! end
+  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(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