I've written a Python/Sage script that computes the discriminant of a general element in $L = \mathbb{Q}(\sqrt{-1+2\sqrt{2}})$.
#Set up the field L
K.<a> = QuadraticField(2)
S.<y> = K[]
L.<b> = K.extension(y^2 - (-1 + 2*a))
L.<c> = L.absolute_field()
#Set up the general integer theta
T.<i,j,k,l> = QQ[]
U.<x> = T[]
U.<x> = U.quotient_ring(U.ideal(L.polynomial()))
int_basis = map(lambda omega:U(omega.lift()), L.integral_basis())
theta = sum(map(mul,zip(int_basis,T.gens()))) #dot product
#Compute the matrix determinant
m = matrix([(theta*x^i).lift().coefficients() for i in range(4)])
traces = [(m^i).trace() for i in range(7)]
trace_matrix = matrix([traces[i:i+4] for i in range(4)])
disc = trace_matrix.det()
The variable disc
should be $-448 = -2^6 \cdot 7$ (the discriminant of $L$) times an integer polynomial in $i,j,k,l$. This doesn't prove that $L = \mathbb{Q}(\sqrt{-1+2\sqrt{2}})$ isn't monogenic, but if we could show that the polynomial is divisible by 2 that would suffice.
Here's how it works. After setting up the field $L$, we find an integral basis for its integers and express a general element $\theta$ in terms of it. In the definition of int_basis
we switch from the number $c$ to the variable $x$ so that Sage treats $\sqrt{-1+2\sqrt{2}}$ and the integral basis as polynomials in $x$; this makes it easier to collect coefficients and compute the matrix of $\theta$ acting on $L$ using the coefficients()
method. In the last block we compute the matrix $(\text{Tr }\theta^{i+j-2})_{i,j}$ and its determinant.
This script could easily be generalized to handle other number fields.
See this question and its answer for a successful application of this method.