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:       lcol = psi_bilinear_matrix_columns(l) 
      enddo
  
      l = psi_bilinear_matrix_order_reverse(k)
      if (l == N_det) cycle
      l = l+1
      ! Fix alpha determinant, loop over betas
      lrow = psi_bilinear_matrix_transp_rows(l) 
      lcol = psi_bilinear_matrix_transp_columns(l) 
      do while ( lrow == krow )
        tmp_det2(:) = psi_det_beta_unique (:, lcol)
        call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int)
        if (degree == 1) then
          call get_mono_excitation_spin(tmp_det(1,2),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_transp_values(l,m) * phase
            tmp_b(h1,p1,m) += ckl
            tmp_b(p1,h1,m) += ckl
          enddo
  
negative.code:   enddo
    !$OMP END DO NOWAIT
    !$OMP CRITICAL
    one_body_dm_mo_beta(:,:,:)  = one_body_dm_mo_beta(:,:,:)  + tmp_b(:,:,:)
    !$OMP END CRITICAL
  
    deallocate(tmp_b)
    !$OMP END PARALLEL
  
  END_PROVIDER
  