Python 3 (PyPy) (n=4 repeated 40 times in 0.43 seconds, n=6 in 36 seconds on TIO)

import itertools
import functools
from collections import Counter
import time

@functools.lru_cache(None)
def gen_rectangles(n):
    x = []
    for t in itertools.product([
        (0,0,0),
        (1,1,1),
        (0,0,1),
        (1,1,0),
        (0,1,0),
        (1,0,1),
    ], repeat=n):
        a = b = c = 0
        for a1, b1, c1 in t:
            a, b, c = a*2+a1, b*2+b1, c*2+c1
        if a < b < c:
            x.append((a, b, c, a^b^c))
    return x

def mask(n, k):
    return ((1<<(1<<n))-1)//((1<<(1<<k))+1)

def swap_parts(n, k1, k2):
    m1 = mask(n, k1)
    m2 = mask(n, k2)
    m12 = m1 ^ m2
    ma =  m1 & ~m2
    mb = ~m1 &  m2
    d = 2**k2 - 2**k1
    return m12, ma, mb, d

def swap(x, m12, ma, mb, d):
    xs = x & m12
    a = x & ma
    b = x & mb
    return x ^ xs ^ (a>>d) ^ (b<<d)

def flip_parts(n, k):
    m = mask(n, k)
    k1 = 1 << k
    return m, k1

def flip(x, m, k1):
    a = x & m
    b = (x >> k1) & m
    return (a << k1) | b

@functools.lru_cache(None)
def parts(n):
    return [flip_parts(n, k) for k in range(n)], [swap_parts(n, k1, k2) for k2 in range(n) for k1 in range(k2)]

def reduce(x, n):
    fp, sp = parts(n)
    for p in fp:
        x1 = flip(x, *p)
        if x1 < x:
            return x1, 0, p
    for p in sp:
        x1 = swap(x, *p)
        if x1 < x:
            return x1, 1, p

def make_order(d):
    d = {a:b[0] for a,b in d.items()}
    C = Counter(d.values())
    empty = d.keys() - d.values()
    l = []
    while empty:
        i = empty.pop()
        j = d.get(i, None)
        if j is None: continue
        l.append(i)
        k = C[j] - 1
        if k == 0:
            empty.add(j)
        else:
            C[j] = k
    return l

@functools.lru_cache(None)
def gen_cubes(n):
    if n == 0:
        return [0, 1]
    if n == 1:
        return [0, 1, 2, 3]
    prev = gen_cubes(n-1)
    half = 2 ** (n-2)
    k = 2 ** (n-1)
    full = 2 ** k - 1
    n1 = n - 1
    bad = {}
    good = []
    for x in prev:
        bc = bin(x).count('1')
        if bc > half: 
            bad[x] = full ^ x, 2, full
            continue
        t = reduce(x, n1)
        if t:
            bad[x] = t
            continue
        good.append(x)
    order = make_order(bad)
    compatible = {}
    masks = []
    for a,b,c,d in gen_rectangles(n):
        if b >= k or c < k: continue
        a, b, c, d = 1 << a, 1 << b, 1 << (c-k), 1 << (d-k)
        masks.append((a|b, c|d, a, d, b, c))
    for b in good:
        l = prev
        for ma, mb, pa1, pb1, pa2, pb2 in masks:
            q = b & mb
            if q == pb1:
                pa = pa1
            elif q == pb2:
                pa = pa2
            else:
                continue
            l = [a for a in l if a&ma != pa]
        compatible[b] = l
    for x in reversed(order):
        pre, disc, p = bad[x]
        l = compatible[pre]
        if disc == 0:
            l = [flip(i, *p) for i in l]
        elif disc == 1:
            l = [swap(i, *p) for i in l]
        else:
            l = [p ^ i for i in l]
        compatible[x] = l
    out = []
    for x, l in compatible.items():
        x <<= k
        out.extend(x|i for i in l)
    return out


if __name__ == '__main__':
    t = time.time()
    for i in range(7):
        print(i, len(gen_cubes(i)))
    print(f'{time.time() - t:.1f} seconds 0-6')
    t = time.time()
    for _ in range(40):
        gen_cubes.cache_clear()
        parts.cache_clear()
        gen_rectangles.cache_clear()
        for i in range(5):
            gen_cubes(i)
    print(f'{time.time() - t:.3f} seconds 0-4 x40')

Try it online!

Based on @Level River St's Ruby answer. The main difference is that this answer exploits additional symmetry.

As a simple example, consider a single square of the hypercube. There are 3 cases:

• All four vertices are colored the same (2 cases like this (flip))
• A single vertex is colored differently (8 cases like this (flip/rotate))
• 2 adjacent vertices are colored 1 (4 cases like this (rotate only))

This answer considers the symmetries of the n-1 hypercube, drastically reducing the amount of work necessary (in the case of n=6, it reduces 94572 cases to only 63 cases)