The number of families of subsets of $\{1,2,...,n\}$ whose union is not the whole set.

I wrote a code in Mathematica that returns the number of families (collections) of subsets of $\{1,2,...,n\}$ whose union is not the whole set. The code can only return values for $n = 1,2,3,4$. The respective values are $2,6,40, 1376$. These are the first four terms of Sloane's A051185 which counts the number of intersecting families. Is this a coincidence or is there some reason why these two counts are equal?

A051185 is the number of (pairwise-) intersecting families. Two subsets of $\{1,\ldots, n\}$ intersect iff the union of their complements is not $\{1,\ldots,n\}$. But you want not just the pairwise unions but the union of all sets in your family to not be the whole set. So you should get a different result for $n=3$: the family $\{\{1,2\}, \{1,3\}, \{2,3\}\}$ is pairwise intersecting, so it should be included in A051185, but its complements $\{\{3\},\{2\},\{1\}\}$ should not be included in your count.
In fact, I get $38$, not $40$, for $n=3$ and $942$, not $1376$, for $n=4$.

If I understand correctly what you're trying to count, your code is wrong. I think the sequence you're after is this one: https://oeis.org/A005530

The problem is equivalent to count the number of covers of a finite set.

It is mentioned at the beginning of the paper Minimal Covers of Finite Sets by Hearne and Wagner that the number is given as the following (where "$n=1$" should be a typo):

Alternatively, the number of possible covers for a set of $N$ elements are

$$ C(N)=\frac{1}{2}\sum_{k=0}^N(-1)^k\binom{N}{k}2^{2^{N-k}}, $$ the first few of which are $1, 5, 109, 32297,\cdots$ (OEIS A003465: Number of ways to cover an n-set. ).

Here is an approach using the Polya Enumeration Theorem (PET) using the same technique as at this MSE link.

We start with the complementary problem counting families that cover all of $[n].$ We get the following generating function for these subsets:

$$\prod_{q=1}^n (1+A_q).$$

Selecting a subset from these of size $k$ is done with the unlabeled set operator $\mathfrak{P}:$

$$Z(P_k)\left(\prod_{q=1}^n (1+A_q)\right).$$

This operator has OGF

$$Z(P_k) = [w^k] \exp\left(\sum_{d\ge 1} (-1)^{d+1} a_d \frac{w^d}{d}\right)$$

The substitution into the cycle index uses

$$a_d = \prod_{q=1}^n (1+A_q^d).$$

We use inclusion-exclusion to compute the terms where not all $A_q$ a present. These are obtained by first subtracting those where a single $A_q$ or more are not present, i.e. by setting that $A_q$ to zero. Then we add in those where a pair of two $A_q$ or more are not present, again by setting these $A_q$ to zero and so on. The remaing $A_q$ are set to one. Hence if $p$ or more of the $A_q$ were not present the substitution becomes

$$a_d = 2^{n-p}.$$

We thus obtain for the substituted cycle index

$$[w^k] \exp\left(\sum_{d\ge 1} (-1)^{d+1} 2^{n-p} \frac{w^d}{d}\right) \\ = [w^k] \exp\left(2^{n-p} \sum_{d\ge 1} (-1)^{d+1} \frac{w^d}{d}\right) \\ = [w^k] \exp\left(2^{n-p} \log(1+w) \right) = [w^k] (1+w)^{2^{n-p}} = {2^{n-p}\choose k}.$$

This means we have by inclusion-exclusion for the problem of the set being covered

$$\sum_{p=0}^n {n\choose p} (-1)^p {2^{n-p}\choose k}.$$

Summing over all possible $k$ (can be anywhere from zero to $2^n$) we finally obtain

$$\sum_{p=0}^n {n\choose p} (-1)^p \sum_{k=0}^{2^n} {2^{n-p}\choose k} = \sum_{p=0}^n {n\choose p} (-1)^p 2^{2^{n-p}}.$$

This yields for the original problem whose answer is being sought

$$\bbox[5px,border:2px solid #00A000]{ 2^{2^n} - \sum_{p=0}^n {n\choose p} (-1)^p 2^{2^{n-p}}.}$$

which is the sequence OEIS A005530

$$2, 6, 38, 942, 325262, 25768825638, \\ 129127208425774833206, \ldots$$

as observed in a companion post.

The Maple code for this (which will compute the first four values) was like this (shown here to clarify what interpretation of the problem is being used):

with(combinat);

ENUM :=
proc(n)
option remember;
local set, res, covered;

    res := 0;

    for set in powerset(powerset(n)) do
        covered := `union`(op(set));

        if nops(covered) < n then
            res := res + 1;
        fi;
    od;

    res;
end;

X := n-> 2^(2^n)
-add(binomial(n,p)*(-1)^p*2^(2^(n-p)), p=0..n);

Addendum. Some thought produces the following simple Perl script which adhereres to first principles (problem definition), yet is capable of computing the first five values of the sequence in less than three seconds, which appears to be the natural boundary of what is possible. We solve the complementary problem as above.

#! /usr/bin/perl -w
#

sub recurse {
    my ($src, $seen, $pos, $n, $cref) = @_;

    my $chk;
    for($chk = 0; $chk < $n; $chk++){
        last if $seen->[$chk] == 0;
    }

    if($chk == $n){
        $$cref += 1 << ((1<<$n)-$pos);
        return;
    }

    return if $pos >= (1 << $n);

    recurse($src, $seen, $pos+1, $n, $cref);

    foreach my $el (@{ $src->[$pos] }){
        $seen->[$el]++;
    }
    recurse($src, $seen, $pos+1, $n, $cref);
    foreach my $el (@{ $src->[$pos] }){
        $seen->[$el]--;
    }

    1;
}

MAIN : {
    my $n = shift || 4;

    my $srcset = [];

    for(my $idx = 2**$n-1; $idx >= 0; $idx--){
        my @data;

        for(my ($bitpos, $idx2) = (0, $idx); 
            $bitpos < $n; $bitpos++){
            my $bit = $idx2 % 2;

            push @data, $bitpos if $bit == 1;
            $idx2 = ($idx2-$bit)/2;
        }

        push @$srcset, \@data;
    }

    my (@vis) = (0) x $n;

    my $result = 0; 
    recurse($srcset, \@vis, 0, $n, \$result);
    $result = 2**(2**$n) - $result;

    print "$result\n";

    1;
}