Fix $n$, and let $V(x)$ be the polynomial from the question, for which we want to show that it is positive on the open interval $(0,1)$. Set \begin{align*} a(n) &= 1024n^2 - 4096/3n + 320\\ b(n) &= 12288n^3 - 14336n^2 + 4864/3n + 256\\ W(x) &= x^4\left(V(x)-(n+1)^2(2n+1)^2x^{24n-1}(1-x)^4(a(n)x + b(n)(1-x))\right). \end{align*}
For all $n\ge12$ we claim that $(1-x)^6$ divides $W(x)$, and that all the coefficients of $W(x)/(1-x)^6$ are nonnegative. In particular, $W(x)$ is positive for $x>0$. As $a(n)$ and $b(n)$ is positive for all $n\ge1$, we see that $V(x)$ is positive for all $0<x<1$, provided that $n\ge12$. The cases $2\le n\le 11$ are easily checked directly, for instance by noting that the coefficients of $(1+x)^{\deg V}V(1/(1+x))$ are positive.
I'm not sure about the easiest kind to prove the assertion about the coefficients of $W(x)/(1-x)^6$. We compute them via \begin{equation} \frac{W(x)}{(1-x)^6}=W(x)\sum_{k\ge0}\binom{k+5}{5}x^k. \end{equation}
All what follows is assisted/confirmed by the Sage code below.
The exponents of $W(x)$ have the form $2ni+j$ for $0\le i\le 12$ and $0\le j\le 10$. Let $a_{i,j}$ be the corresponding coefficient. It is a polynomial in $n$.
We have \begin{equation} \frac{W(x)}{(1-x)^6} = \sum_{i,j}a_{i,j}(n)x^{2ni+j}\sum_k\binom{k+5}{5}x^k. \end{equation}
With $r\ge0$ and $0\le s\le 2n-1$, the coefficient of $x^{2nr+s}$ in $W(x)\frac{1}{(1-x)^6}$ is the sum of all $a_{i,j}(n)\cdot\binom{nr+s-ni-j+y}{5}$ for all pairs $(i,j)$ where $i<r$ and $0\le j\le 10$, or $i=r$ and $j\le s$.
We distinguish the cases $0\le s\le 9$ from the cases $s\ge10$. In these latter cases, we write $s=10+y$. Recall that $s\le 2n-1$, so $10+y=2n-1-o$ for a nonnegative $o$. Upon replacing $n$ with $(11+y+1)/2$, we get polynomials in $n$ and $o$.
In the former cases, the coefficients are polynomials in $n$, which have positive coefficients upon replacing $n$ with $n+12$. So these coefficients are positive for all $n\ge12$.
In the latter cases it turns out that in all but one case the coefficients are positive. So for nonnegative $y$ and $o$, the values are nonnegative.
For the single exception we see that if we multiply it with a suitable polynomial, the resulting coefficients are positive.
Remark (answering Jason's question): It is a known fact that a polynomial $V$ is positive on $(0,1)$ if and only if it is a nonnegative linear combination of polynomials $x^i(1-x)^j$. The problem is that it may involve terms where $i+j$ is bigger than $\deg V$. (This doesn't happen here, though.) I used an LP solver to play a little with $V$ and $\frac{V}{(1-x^2)^4x^{2n-2}}$. From that a pattern showed up which lead to a solution.
By the way, $V$ actually seems to be a nonnegative linear combination of $x^i(1-x)^{\deg V-i}$. This is equivalent to say that all the coefficients of $(1+x)^{\deg V}V(1/(1+x))$ are nonnegative. It is not hard to compute explicit expressions for these coefficients, but I don't see an easy arguments why they can't be negative.
# Formally compute W
var('X N')
P = (1+X^(4*N-1))*(1+X^(2*N))*(1-X^(2*N+1))
Q = (1+X^(2*N+1))*(1-X^(2*N+2))
V = P*Q^2*P.diff(X,2)+(P*Q.diff(X))^2-P^2*Q*Q.diff(X,2)-(P.diff(X)*Q)^2
a = 1024*N^2 - 4096/3*N + 320
b = 12288*N^3 - 14336*N^2 + 4864/3*N + 256
W = X^4*(V - (2*N + 1)^2 *(N + 1)^2*X^(24*N-1)*(1-X)^4*(a*X + b*(1-X)))
# Check that (1-X)^6 divides W(X)
print all(W.diff(X,i)(X=1).polynomial(QQ) == 0 for i in [0..5])
# Compute the coefficients a_ij for the exponents 2ni+j of W. Somewhat
# clumsy, as I don't know how to deal with polynomials where exponents
# are symbolic expressions.
#
# l is the list of summands of W
l = [z.canonicalize_radical() for z in W.expand().operands()]
K.<n> = QQ[]
aij = {(i,j):K(0) for i in [0..12] for j in [0..10]}
for term in l:
c = term(X=1) # get coefficient of term
e = (term.diff(X)/c)(X=1) # get exponent of term
c = K(c.polynomial(QQ)) # convert c to proper polynomial in n
i, j = ZZ(e.diff(N))//2, ZZ(e(N=0)) # Clumsy method to compute pairs (i,j)
aij[i,j] += c
# Check if coefficients aij[i,j] were correctly computed
Wnew = sum(c(n=N)*X^(2*N*i+j) for (i,j),c in aij.items())
print (W-Wnew).canonicalize_radical().is_trivial_zero()
def bino(k): # binomial coefficient binom(-6,k)=binom(k+5,5)
return prod(k+5-z for z in range(5))/120
# compute the coefficients of W(X)/(1-X)^6
K = K.extend_variables(('y', 'o'))
K.inject_variables()
for r in [0..12]:
for s in [0..10]:
d = 1 if s == 10 else 0
# compute coefficient of X^(2nr+s+d*y)
f = sum(aij[i,j]*bino(2*n*(r-i)+s+d*y-j) for i in [0..r]
for j in [0..s if i == r else 10])
# check non-negativity of the polynomial f
if d == 0:
# Checks if the coefficient of X^(2nr+s), which is a
# polynomial in n, has nonnegative coefficients upon replacing
# n with n+12
f = f(n=n+12)
if min(f.coefficients()+[0]) < 0:
print "False"
else:
# The coefficient of X^(2nr+10+y), which is a polynomial in n
# and y, has to be nonnegative for 0 <= 10+y <= 2n-1, so 10+y
# = 2n-1-o for non-negative o. So upon replacing n with
# (11+y+o)/2, we get a polynomial in y and o which has to be
# nonnegative for all nonnegative y and o. In all but one
# case, this holds because the coeffcients are non-negative.
f = f(n=(11+y+o)/2)
if min(f.coefficients()+[0]) < 0:
c = o^2 + 23*o*y + 1360*y^2 + 99*o + 340*y + 1675
print min(c.coefficients()+(c*f).coefficients()) >= 0