anchor.diff:          l_a += 1
         enddo
   
  +      n_doubles=0
         do while ( l_a < psi_bilinear_matrix_columns_loc(lcol+1) )
           lrow = psi_bilinear_matrix_rows(l_a)
  -        if (is_single_a(lrow)) then
  -          tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
  -
  -          call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij)
  -          call get_s2(tmp_det,tmp_det2,N_int,sij)
  -          do l=1,N_st
  -            v_t(l,k_a) += hij * u_t(l,l_a)
  -            v_t(l,l_a) += hij * u_t(l,k_a)
  -            s_t(l,k_a) += sij * u_t(l,l_a)
  -            s_t(l,l_a) += sij * u_t(l,k_a)
  -          enddo
  +        if (.not.is_single_a(lrow)) then
  +          continue
  +        else
  +          n_doubles = n_doubles+1
  +          doubles(n_doubles) = lrow
  +          idx(n_doubles) = l_a
           endif
  -        l_a += 1
  -
  -      enddo
  +        l_a = l_a+1
  +      enddo
  +
  +      do k=1,n_doubles
  +        lrow = doubles(k)
  +        l_a  = idx(k)
  +        tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
  +        call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij)
  +        call get_s2(tmp_det,tmp_det2,N_int,sij)
  +        do l=1,N_st
  +          v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
  +          v_t(l,l_a) = v_t(l,l_a) + hij * u_t(l,k_a)
  +          s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a)
  +          s_t(l,l_a) = s_t(l,l_a) + sij * u_t(l,k_a)
  +        enddo
  +      enddo
  +
       enddo
   
       ! Diagonal contribution
  
positive.code:       tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
        call i_H_j_mono_spin( tmp_det, tmp_det2, N_int, 1, hij)
        do l=1,N_st
          v_t(l,l_a) += hij * u_t(l,k_a)
          v_t(l,k_a) += hij * u_t(l,l_a)
          ! single => sij = 0 
        enddo
      enddo
  
negative.code:     tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow)
      tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol)
      
      do k=1,singles_a(0,krow)
        is_single_a( singles_a(k,krow) ) = .True.
      enddo
  