anchor.diff:        enddo
         call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, doubles(i)), N_int, 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)
  +        v_t(l,l_a) = v_t(l,l_a) + hij * u_t(l,k_a)
  +        v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
           ! same spin => sij = 0 
         enddo
       enddo
  
positive.code:   logical, allocatable           :: is_single_b(:)
    integer                        :: maxab, n_singles_max
    double precision, allocatable  :: u_t(:,:), v_t(:,:), s_t(:,:)
  
    maxab = max(N_det_alpha_unique, N_det_beta_unique)
    allocate( buffer(N_int,maxab),                                     &
  
negative.code:         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 (.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 = 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
  