anchor.diff:          call get_mono_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int)
           call decode_exc_spin(exc,h1,p1,h2,p2)
           do m=1,N_states
  -          ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
  +          ckl = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(l,m) * phase
             tmp_a(h1,p1,m) += ckl
             tmp_a(p1,h1,m) += ckl
           enddo
  
positive.code:     !$OMP  psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
    allocate(tmp_a(mo_tot_num,mo_tot_num,N_states), tmp_b(mo_tot_num,mo_tot_num,N_states) )
    tmp_a = 0.d0
    tmp_b = 0.d0
    !$OMP DO SCHEDULE(guided)
    do k=1,N_det
      krow = psi_bilinear_matrix_rows(k) 
      kcol = psi_bilinear_matrix_columns(k) 
      tmp_det(:,1) = psi_det_alpha_unique(:,krow)
      tmp_det(:,2) = psi_det_beta_unique (:,kcol)
      call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
      do m=1,N_states
        ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
        do l=1,elec_alpha_num
          j = occ(l,1)
          tmp_a(j,j,m) += ck
        enddo
        do l=1,elec_beta_num
          j = occ(l,2)
          tmp_b(j,j,m) += ck
        enddo
      enddo
  
      if (k == N_det) cycle
      l = k+1
      lrow = psi_bilinear_matrix_rows(l) 
      lcol = psi_bilinear_matrix_columns(l) 
      ! Fix beta determinant, loop over alphas
  
negative.code:    ! Alpha and beta one-body density matrix for each state
     END_DOC
  
     integer                        :: j,k,l,m,k_a,k_b
     integer                        :: occ(N_int*bit_kind_size,2)
     double precision               :: ck, cl, ckl
     double precision               :: phase
  