anchor.diff:        is_single_b(N_det_beta_unique),                                &
         idx(maxab), idx0(maxab))
     
  +  do i=1,maxab
  +    idx0(i) = i
  +  enddo
  +
  +  ! Prepare the array of all alpha single excitations
  +  ! -------------------------------------------------
  +
  +  n_singles_max = 0
  +  do i=1,N_det_alpha_unique
  +    spindet(1:N_int) = psi_det_alpha_unique(1:N_int, i)
  +    call get_all_spin_singles(                                       &
  +        psi_det_alpha_unique, idx0, spindet, N_int, N_det_alpha_unique,&
  +        singles, n_singles)
  +    n_singles_max = max(n_singles_max, n_singles)
  +  enddo
  +
  +  allocate (singles_a(0:n_singles_max, N_det_alpha_unique))
  +  do i=1,N_det_alpha_unique
  +    spindet(1:N_int) = psi_det_alpha_unique(1:N_int, i)
  +    call get_all_spin_singles(                                       &
  +        psi_det_alpha_unique, idx0, spindet, N_int, N_det_alpha_unique,&
  +        singles_a(1,i), singles_a(0,i))
  +  enddo
  +
     v_0 = 0.d0
     do k_a=1,N_det
       
  
positive.code: 
    enddo
  
    !----
  
  !  k_a=1
  !  do i=1,N_det_beta_unique
  !
  !    ! Select a beta determinant
  !    ! -------------------------
  !
  !    spindet(1:N_int) = psi_det_beta_unique(1:N_int, i)
  !    tmp_det(1:N_int,2) = spindet(1:N_int)
  !
  !    call get_all_spin_singles(                                       &
  !        psi_det_beta_unique, idx0, spindet, N_int, N_det_beta_unique, &
  !        singles, n_singles )
  !
  !    do j=1,n_singles
  !      is_single_a( singles(j) ) = .True.
  !    enddo
  !
  !    ! For all alpha.beta pairs with the selected beta
  !    ! -----------------------------------------------
  !
  !    kcol = psi_bilinear_matrix_columns(k_a)
  !    do while (kcol < i)
  !      k_a = k_a+1
  !      if (k_a > N_det) exit
  !      kcol = psi_bilinear_matrix_columns(k_a)
  !    enddo
  !
  !    do while (kcol == i)
  !
  !      krow = psi_bilinear_matrix_rows(k_a)
  !      tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
  !
  !      ! Loop over all alpha.beta pairs with a single exc alpha
  !      ! ------------------------------------------------------
  !
  !      l_a = k_a+1
  !      if (l_a > N_det) exit
  !      lrow = psi_bilinear_matrix_rows(l_a)
  !      lcol = psi_bilinear_matrix_columns(l_a)
  !
  !      do while (lrow == krow)
  !
  !        ! Loop over all alpha.beta pairs with a single exc alpha
  !        ! ------------------------------------------------------
  !        if (is_single_a(lrow)) then
  !
  !           tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int,lrow)
  !          
  !           ! Build list of singly excited beta
  !           ! ---------------------------------
  !
  !           m_b = psi_bilinear_matrix_order_reverse(l_a)
  !           m_b = m_b+1
  !           j=1
  !           do while ( (mrow == lrow) )
  !             mcol = psi_bilinear_matrix_transp_columns(m_b)
  !             buffer(1:N_int,j) = psi_det_beta_unique(1:N_int,mcol)
  !             idx(j) = mcol
  !             j = j+1
  !             m_b = m_b+1
  !             if (m_b <= N_det) exit
  !             mrow = psi_bilinear_matrix_transp_rows(m_b)
  !           enddo
  !           j=j-1
  !
  !           call get_all_spin_singles(                                &
  !               buffer, idx, tmp_det(1,2), N_int, j,                  &
  !               doubles, n_doubles)
  !
  !           ! Compute Hij for all doubles 
  !           ! ---------------------------
  !
  !           m_b = psi_bilinear_matrix_order(l_a)+1
  !           mcol = psi_bilinear_matrix_transp_columns(m_b)
  !           do j=1,n_doubles
  !             tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, doubles(j) )
  !             call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij)
  !             do while (mcol /= doubles(j))
  !               m_b = m_b+1
  !               if (m_b > N_det) exit
  !               mcol = psi_bilinear_matrix_transp_columns(m_b)
  !             enddo
  !             m_a = psi_bilinear_matrix_order_reverse(m_b)
  !             v_0(m_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
  !             v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(m_a,1:N_st)
  !           enddo
  !
  !        endif
  !        l_a = l_a+1
  !        if (l_a > N_det) exit
  !        lrow = psi_bilinear_matrix_rows(l_a)
  !        lcol = psi_bilinear_matrix_columns(l_a)
  !      enddo
  !
  !      k_b = k_b+1
  !      if (k_b > N_det) exit
  !      kcol = psi_bilinear_matrix_transp_columns(k_b)
  !    enddo
  !    
  !    do j=1,n_singles
  !      is_single_a( singles(j) ) = .False.
  !    enddo
  !
  !  enddo  
  
      
  end
  
  
  
negative.code: 
            call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij)
            v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
            v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
          endif
          l_a += 1
  
  