How many non-isomorphic binary structures on the set of $n$ elements?

Following up the citations at the OEIS page in my comment, it looks like much is to be learned from Michael A Harrison, The number of isomorphism types of finite algebras, Proc Amer Math Soc 17 (1966) 731-737. The formulas there are elementary but very long, so I won't type any of them out here. I believe the paper is freely available on the AMS website.

If the sources are saying what I think they're saying, for $n=3$ you get 3,330.

EDIT: In a comment, Doug asks about a formula at the OEIS page. I don't think it will fit in a comment, so I'll try to write it out here. Let $${\rm fix\ }A(s_1,s_2,\dots)=\prod_{i,j\ge1}\sum_{d\mid{\rm lcm}(i,j)}(ds_d)^{s_is_j\gcd(i,j)}$$ Then $$a_n=\sum_{s_1+2s_2+\dots=n}{{\rm fix\ }A(s_1,s_2,\dots)\over1^{s_1}s_1!2^{s_2}s_2!\dots}$$

It seems to me important to point out that the above formula iterating over all partition shapes is not the best we can do. What we are looking for here is an enumeration of all $N\times N$ matrices with entries from $1$ to $N$ under the symmetric group acting on the entries and the rows and columns at the same time. We can use Burnside directly and calculate the appropriate cycle index (which represents the action of the symmetric group on ordered pairs) from the cycle index of the symmetric group, as was done at this MSE link. The cycle index of the symmetric group can be computed with a trick that goes back to Lovasz: $$ Z(S_0) = 1 \quad \text{and} \quad Z(S_n) = \frac{1}{n} \sum_{l=1}^n a_l Z(S_{n-l})$$ and we do not need to iterate over all partitions to compute this index. (It must be pointed out however that this cycle index contains terms with disjoint cycle decompositions corresponding to all partitions.)

Now Burnside is quite simple here. For each cycle of a single permutation $P$ of the matrix elements induced by the action of a permutation $Q$ from the symmetric group on $N$ elements we need to determine the assignments to the slots on the cycle that are fixed by $P$. But these are precisely the assignments that consist of cyclical repetitions of complete cycles from $Q$ i.e. we can place any cycle from $Q$ on a cycle from $P$ as long as the length of the former divides the length of the latter. Every such pair produces $r$ possible assignments where $r$ is the length of the cycle from $Q.$ This simple observation suffices to apply Burnside.

We get the following sequence of values: $$ 1, 10, 3330, 178981952, 2483527537094825, 14325590003318891522275680,\\ 50976900301814584087291487087214170039,\\ 155682086691137947272042502251643461917498835481022016,\\ 541851439802559836957713164869818405872834954135521300809902639457510935,\ldots$$ which agrees with the OEIS entry.

The above values were computed with the following Maple code.

pet_cycleind_symm :=
proc(n)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

bs_binop :=
proc(n)
option remember;
local dsjc, varp1, varp2, var, p, q, len,
    l1, l2, res;

    if n=0 then return 1; fi;
    if n=1 then return 1; fi;

    res := 0;

    for dsjc in pet_cycleind_symm(n) do
        p := 1;

        for varp1 in indets(dsjc) do
            l1 := op(1, varp1);

            for varp2 in indets(dsjc) do
                l2 := op(1, varp2);

                len := lcm(l1,l2); q := 0;

                for var in indets(dsjc) do
                    if len mod op(1, var) = 0 then
                        q := q  +
                        op(1, var)*degree(dsjc, var);
                    fi;
                od;

                p := p *
                q^(degree(dsjc, varp1) *
                   degree(dsjc, varp2)
                   *l1*l2/len);
            od;
        od;

        res := res + p*lcoeff(dsjc);
    od;

    res;
end;