I have not yet had time to read the Clapham paper, but I can perhaps be of use with a related statistic, the number of non-isomorphic graphs on $n$ nodes. There is a wealth of data on this at the OEIS, sequence A000088. Your two colors essentially represent edges being switched on and off, so this is the same statistic.

As you can see from the OEIS entry this is a classic computation, and extremely straightforward. Compute the cycle index of the edge permutation group induced by the $n!$ permutations of the vertex permutation group, substitute with $1+z$ for edges being on and off to get the generating function of non-isomorphic graphs classified by the number of edges. Use Burnside setting all variables to two if only the count rather than a classification is desired. The cycle index of the vertex permutation 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}).$$
Now to turn a vertex permutation into an edge permutation the cycle structure of the former is sufficient, yielding a reduced complexity without which we'd be out of luck. Simply examine and enumerate all possible pairs of vertices and imagine them moving in parallel (i.e., forming an edge) along their respective cycles until the starting configuration is reached. This happens after $\operatorname{LCM}(l_1, l_2)$ steps assuming $l_{1,2}$ are the lengths of the two cycles of the vertex permutation. There are two cases when the two vertices lie on the same cycle, depending on whether the length is even or odd. In the former case there are $l/2$ pairs that return to their origin after $l/2$ steps and the rest take $l$ steps, where $l$ is the length of the cycle. Finally there is some overcounting to take into consideration, as every cycle in the edge permutation is overcounted by a factor equal to its length. That is all.

Here is a list of values to peruse.

01 1
02 2
03 4
04 11
05 34
06 156
07 1044
08 12346
09 274668
10 12005168
11 1018997864
12 165091172592
13 50502031367952
14 29054155657235488
15 31426485969804308768
16 64001015704527557894928
17 245935864153532932683719776
18 1787577725145611700547878190848
19 24637809253125004524383007491432768
20 645490122795799841856164638490742749440
21 32220272899808983433502244253755283616097664
22 3070846483094144300637568517187105410586657814272
23 559946939699792080597976380819462179812276348458981632
24 195704906302078447922174862416726256004122075267063365754368
25 131331393569895519432161548405816890146389214706146483380458576384
26 169487618400693135671278778610295749756246061427513800357039083537927168
27 421260006519643885757174107650953992882782878952295704539600450662355704738816
28 2019295999678571395728328980890810345860807065053566265347514137546665672316696821760
29 18691352722478956041683441055221773100906878077027397169675907651818104181986752359177684992
30 334494316309257669249439569928080028956631479935393064329967834887217734534880582749030521599504384

The Maple code to compute these is as follows.

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;
pet_cycleind_edg :=
proc(n)
option remember;
local all, term, termvars, res, l1, l2, inst1, u, v,
uidx, vidx;
if n=0 or n=1 then return 1; fi;
all := 0:
for term in pet_cycleind_symm(n) do
termvars := indets(term); res := 1;
# edges on different cycles of different sizes
for uidx to nops(termvars) do
u := op(uidx, termvars);
l1 := op(1, u);
for vidx from uidx+1 to nops(termvars) do
v := op(vidx, termvars);
l2 := op(1, v);
res := res *
a[lcm(l1, l2)]
^((l1*l2/lcm(l1, l2))*
degree(term, u)*degree(term, v));
od;
od;
# edges on different cycles of the same size
for u in termvars do
l1 := op(1, u); inst1 := degree(term, u);
# a[l1]^(1/2*inst1*(inst1-1)*l1*l1/l1)
res := res *
a[l1]^(1/2*inst1*(inst1-1)*l1);
od;
# edges on identical cycles of some size
for u in termvars do
l1 := op(1, u); inst1 := degree(term, u);
if type(l1, odd) then
# a[l1]^(1/2*l1*(l1-1)/l1);
res := res *
(a[l1]^(1/2*(l1-1)))^inst1;
else
# a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
res := res *
(a[l1/2]*a[l1]^(1/2*(l1-2)))^inst1;
fi;
od;
all := all + lcoeff(term)*res;
od;
all;
end;
pet_varinto_cind :=
proc(poly, ind)
local subs1, subsl, polyvars, indvars, v, pot;
polyvars := indets(poly);
indvars := indets(ind);
subsl := [];
for v in indvars do
pot := op(1, v);
subs1 :=
[seq(polyvars[k]=polyvars[k]^pot,
k=1..nops(polyvars))];
subsl := [op(subsl), v=subs(subs1, poly)];
od;
subs(subsl, ind);
end;
gf :=
proc(n)
option remember;
local gf;
gf := pet_varinto_cind(1+z, pet_cycleind_edg(n));
expand(gf);
end;
Q :=
proc(n)
option remember;
local cind;
cind := pet_cycleind_edg(n);
subs([seq(v=2, v in indets(cind))], cind);
end;

Here is the GF of non-isomorphic graphs on five vertices classified by the number of edges:

$${z}^{10}+{z}^{9}+2\,{z}^{8}+4\,{z}^{7}+6\,{z}^{6}
+6\,{z}^{5}+6\,{z}^{4}+4\,{z}^{3}+2\,{z}^{2}+z+1.$$

A similar computation can be used to count the number of binary relations on a set of $n$ elements and is documented at this MSE link. This MSE link II shows how to enumerate graphs that permit self-loops.