Acyclic orientations of an n-dimensional cube

Question

The goal of this challenge is to check and extend the OEIS sequence A334248: Number of distinct acyclic orientations of the edges of an n-dimensional cube.

Take an n-dimensional cube (if n=1, this is a line; if n=2, a square; if n=3, a cube; if n=4, a hypercube/tesseract; etc), and give a direction to all of its edges so that no cycles are formed. The terms of the sequence are the number of (rotationally and reflectively) distinct orientations that can be obtained.

Example

For n=2, there are 14 ways to assign directions to the edges of a square with no cycles formed:

However, only three of these are distinct if they are rotated and reflected:

So when n=2, the answer is 3.

---

The winner of this fastest-code challenge will be the program that calculates the most terms when run on my computer (Ubuntu, 64GB RAM) in one hour. If there is a tie, the program that computes the terms in the shortest time will be the winner.

Answer 1

Wolfram Language (Mathematica), n = 4, <1 second on my computer

(* Elements of the hyperoctahedral group,represented by matrices *)
HyperoctahedralGroup[n_] := 
  Join @@ Table[
    Permute[IdentityMatrix@n, perm]*sign, {perm, 
     Permutations@Range@n}, {sign, Tuples[{1, -1}, n]}];

(* Action of group elements on a vertex of the hypercube,
where each vertex is represented by an integer from 1 to 2^n *)
ActOn[n_, g_, v_] := 
  FromDigits[(1 + (IntegerDigits[v - 1, 2, n]*2 - 1) . g)/2, 2] + 1;

(* Orbits of vertices under such action *)
Orbits[n_, g_] := 
  ConnectedComponents@Table[v <-> ActOn[n, g, v], {v, 2^n}];

(* Whether any orbit contains adjacent vertices *)
OrbitConatinsAdjacentVertices[edges_, orbits_] := 
  Or @@ Or @@@ Outer[SubsetQ, orbits, edges, 1];

(* Merging vertices in each orbit into a single vertex *)
MergeVertices[edges_, orbits_] := 
  Graph@DeleteDuplicates[
    Sort /@ (edges /. (Alternatives@##2 -> #1 & @@@ orbits))];

(* Number of acyclic orientations fixed by some group element *)
InvariantAcyclicOrientations[n_, edges_, g_] := 
  With[{orbits = Orbits[n, g]}, 
   If[OrbitConatinsAdjacentVertices[edges, orbits], 0, 
    Abs[ChromaticPolynomial[MergeVertices[edges, orbits], -1]]]];

(* OEIS sequence A334248 *)
A334248[n_] := 
  With[{edges = List @@@ EdgeList@HypercubeGraph@n}, 
   Sum[InvariantAcyclicOrientations[n, edges, g], {g, 
      HyperoctahedralGroup[n]}]/(2^n*n!)];

Do[Print[n -> A334248[n]], {n, 1, 4}]

Try it online!

Not really optimized. Takes less than 1 seconds on my computer to calculate the first four terms, but takes forever for the fifth term. It seems that Mathematica cannot compute the chromatic polynomial of HypercubeGraph[5] in a reasonable time.

Based on Peter Kagey's deleted answer (which did not consider rotations and reflections), and uses Burnside's lemma.