I would like to contribute some enrichment material to facilitate additional exploration of this problem. The most important observation is that what we have here is an instance of Power Group Enumeration, with the group acting on the slots where the $n$ colors are placed being the edge permutation group $E$ of the cube and the group acting on the colors being the symmetric group on $n$ elements $S_n$.
We can compute the number of configurations by Burnside's lemma which says to average the number of assignments fixed by the elements of the power group, which has $n!\times |E|!$ elements and $|E|=24$. But this number is easy to compute. Suppose we have a permutation $\alpha$ from $E$ and a permutation $\beta$ from $S_n.$ If we place the appropriate number of complete, directed and consecutive copies of a cycle from $\beta$ on a cycle from $\alpha$ then this assignment is fixed under the power group action, and this is possible iff the length of $\beta$ divides the length of $\alpha$ and there are as many assignments as the length of $\beta.$
To do this we first need the cycle index of $E$, which we now compute. There is the identity, which contributes $a_1^{12}.$ Rotations about an axis passing through a pairs of opposite vertices contribute $4\times 2 a_3 ^4.$ Rotations about an axis passing through opposite faces contribute $3 \times (2 a_4^3 + a_2^6).$ Finally rotations about an axis passing through midpoints of opposite edges contribute $6\times a_1^2 a_2^5.$
This gives the cycle index $Z(E) = \frac{1}{24} \left(a_1^{12} + 8 a_3^4 + 6 a_4^3 + 3 a_2^6 + 6 a_1^2 a_2^5\right).$
Let's verify this cycle index before we proceed. For edge colorings with $n$ colors where the colors are not being permuted we obtain the formula $\frac{1}{24} \left(n^{12} + 8 n^4 + 6 n^3 + 3 n^6 + 6 n^7\right).$ This gives the sequence $1, 218, 22815, 703760, 10194250, 90775566, 576941778, 2863870080,\ldots$ which points us to OEIS A060530, where we find that indeed we have the right cycle index.
Now the Burnside computation is best done with a CAS, here is the 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; cube_edge_cind := 1/24*(a[1]^12 + 8*a[3]^4 + 6*a[4]^3 + 3*a[2]^6 + 6*a[1]^2*a[2]^5); cube_edge_colorings := proc(n) option remember; local idx_cols, res, term_a, term_b, v_a, v_b, inst_a, inst_b, len_a, len_b, p, q; if n = 1 then return 1 fi; idx_cols := pet_cycleind_symm(n); res := 0; for term_a in cube_edge_cind do for term_b in idx_cols do p := 1; for v_a in indets(term_a) do len_a := op(1, v_a); inst_a := degree(term_a, v_a); q := 0; for v_b in indets(term_b) do len_b := op(1, v_b); inst_b := degree(term_b, v_b); if len_a mod len_b = 0 then q := q + len_b*inst_b; fi; od; p := p*q^inst_a; od; res := res + lcoeff(term_a)*lcoeff(term_b)*p; od; od; res; end;
This will produce the following sequence of values for colorings with $n$ colors with the symmetric group acting on the colors: $1, 114, 3891, 29854, 87981, 143797, 170335, 177160,\ldots$
Addendum Jan 29 2019. The above computation is easily adapted to produce a generating function, making it possible not only to count these colorings but also to classify them by the numbers of interchangeable colors used. The technique is the same (cover the cycles of a term from $Z(E)$ by cycles obtained from a permutation belonging to $S_n$), only this time we record the colors being used. With the colors being interchangeable we collect the contributions corresponding to the same partition of the twelve edges according to the colors into one term. This is the result for three colors:
${P_{{1}}}^{12}+29\,{P_{{1}}}^{6}{P_{{2}}}^{6} +38\,{P_{{1}}}^{5}{P_{{2}}}^{7}+27\,{P_{{1}}}^{4}{P_{{2}}}^{8} +282\,{P_{{1}}}^{4}{P_{{2}}}^{4}{P_{{3}}}^{4} +13\,{P_{{1}}}^{3}{P_{{2}}}^{9} \\+1170\,{P_{{1}}}^{3}{P_{{2}}}^{4}{P_{{3}}}^{5} +412\,{P_{{1}}}^{3}{P_{{2}}}^{3}{P_{{3}}}^{6} +5\,{P_{{1}}}^{2}{P_{{2}}}^{10} +370\,{P_{{1}}}^{2}{P_{{2}}}^{5}{P_{{3}}}^{5} \\+600\,{P_{{1}}}^{2}{P_{{2}}}^{4}{P_{{3}}}^{6} +340\,{P_{{1}}}^{2}{P_{{2}}}^{3}{P_{{3}}}^{7} +77\,{P_{{1}}}^{2}{P_{{2}}}^{2}{P_{{3}}}^{8} +P_{{1}}{P_{{2}}}^{11} +236\,P_{{1}}{P_{{2}}}^{5}{P_{{3}}}^{6} \\+170\,P_{{1}}{P_{{2}}}^{4}{P_{{3}}}^{7} +85\,P_{{1}}{P_{{2}}}^{3}{P_{{3}}}^{8} +30\,P_{{1}}{P_{{2}}}^{2}{P_{{3}}}^{9} +5\,P_{{1}}P_{{2}}{P_{{3}}}^{10}.$
The Maple code goes as follows.
with(combinat); cube_edge_cind := 1/24*(a[1]^12 + 8*a[3]^4 + 6*a[4]^3 + 3*a[2]^6 + 6*a[1]^2*a[2]^5); 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_indets2rep := proc(ip) local rep, var, deg, pos; rep := []; pos := 1; for var in indets(ip) do for deg to degree(ip, var) do rep := [op(rep), [seq(s, s=pos..pos+op(1, var)-1)]]; pos := pos + op(1, var); od; od; rep; end; cube_edge_colorings_gf := proc(n) option remember; local idx_cols, rep, res, term_a, term_b, v_a, v_b, inst_a, len_a, len_b, p, q, parts, alldeg, cols; if n = 1 then return P[1]^12 fi; idx_cols := pet_cycleind_symm(n); res := 0; for term_a in cube_edge_cind do for term_b in idx_cols do rep := pet_indets2rep(term_b); p := 1; for v_a in indets(term_a) do len_a := op(1, v_a); inst_a := degree(term_a, v_a); q := 0; for v_b in rep do len_b := nops(v_b); if len_a mod len_b = 0 then q := q + len_b* mul(P[col], col in v_b) ^(len_a/len_b); fi; od; p := p*q^inst_a; od; res := res + lcoeff(term_a)*lcoeff(term_b)*p; od; od; res; parts := 0; for cols in expand(res) do alldeg := sort(map(v -> degree(cols, v), [op(indets(cols))])); parts := parts + lcoeff(cols)* mul(P[v]^alldeg[v], v=1..nops(alldeg)); od; parts; end; cube_edge_colorings := proc(n) option remember; local gf, sl; gf := cube_edge_colorings_gf(n); sl := [seq(v = 1, v in indets(gf))]; subs(sl, gf); end;