To demonstrate that finding the optimal bound is hard, I'll work out the case $n=3$. We already found out that the optimum will occur at a point of the form $(x,x,y)$ with $x>y$. I'll normalize $G=1$, so $y=1/x^2$. Then
$$A = \frac{2x + x^{-2}}{3} \quad H = \frac{3}{2x^{-1} + x^2}.$$
We want to optimize
$$\theta = \frac{1-H}{A-H} = \frac{3 x^2 (2+x)}{2(1+x+x^2)^2}.$$
Macavity's bound, which is very slick and awesome, is to replace $H$ by $G^3/A^2 = 1/A^2$ and consider $\phi = \frac{1-1/A^2}{A-1/A^2}$ instead.
Here is a plot of both $\theta$ and $\phi$ for $1 \leq x \leq 3$ ($\theta$ is blue, $\phi$ is red). Note that $\phi > \theta$, $\phi$ is decreasing and $\phi(1) = 2/3$, as Macavity computes. However, $\max \theta < 2/3$.

We have
$$\frac{d \theta}{dx} = \frac{3 x (4 + 3 x - 3 x^2 - x^3)}{2 (1 + x + x^2)^3}$$
The global maximum is at the root of $4+3x-3x^2-x^3$ which is roughly $1.36147$. The corresponding optimum for $\theta$ is roughly $0.52605$. (Of course, you can get closed form answers using the cubic formula, but they aren't enlightening.)
I still think that finding nice functions $p_n$ which work is a good project, but I don't think you'll get a closed form for the optimal answer.
Some more thoughts: Repeating this computation with $n$ in place of $3$, we have
$$\frac{d \theta}{d x} = \frac{ n x^{n-2} \cdot f(x)}{(n-1) (1-x^n)^3}$$
where
$$f(x) := (n^2-2n+1) - n^2 x + (n^2+2n-2) x^n - n^2 x^{n+1}.$$
By Descartes' Rule of signs, $f$ only has four positive roots. We can compute that $f(x)$ has a triple root at $x=1$, which cancels the triple root of $(1-x^n)^3$, so there is only one positive root to care about.
I set $r=x/y$, since this is independent of our choice of normalization (so $r=x^n$ if we normalize $G=1$.) Computations for $3 \leq n \leq 100$ suggest that $r$ grows something like $n \log n$. I've copy pasted the numbers below in the format $\{ n, r \}$ case someone wants to play with them:
{{3, 2.52361}, {4, 4.34916}, {5, 6.40716}, {6, 8.65639}, {7, 11.0686}, {8, 13.6232}, {9, 16.3041}, {10, 19.0987}, {11, 21.9967}, {12, 24.9895}, {13, 28.0699}, {14, 31.2317}, {15, 34.4694}, {16, 37.7785}, {17, 41.1546}, {18, 44.5942}, {19, 48.094}, {20, 51.6509}, {21, 55.2622}, {22, 58.9256}, {23, 62.6388}, {24, 66.3997}, {25, 70.2064}, {26, 74.0572}, {27, 77.9504}, {28, 81.8847}, {29, 85.8585}, {30, 89.8706}, {31, 93.9197}, {32, 98.0047}, {33, 102.124}, {34, 106.278}, {35, 110.464}, {36, 114.683}, {37, 118.932}, {38, 123.212}, {39, 127.521}, {40, 131.859}, {41, 136.225}, {42, 140.618}, {43, 145.039}, {44, 149.485}, {45, 153.957}, {46, 158.454}, {47, 162.975}, {48, 167.521}, {49, 172.09}, {50, 176.681}, {51, 181.296}, {52, 185.933}, {53, 190.591}, {54, 195.271}, {55, 199.971}, {56, 204.693}, {57, 209.434}, {58, 214.196}, {59, 218.976}, {60, 223.776}, {61, 228.595}, {62, 233.433}, {63, 238.289}, {64, 243.163}, {65, 248.054}, {66, 252.963}, {67, 257.889}, {68, 262.833}, {69, 267.792}, {70, 272.769}, {71, 277.761}, {72, 282.77}, {73, 287.794}, {74, 292.834}, {75, 297.89}, {76, 302.96}, {77, 308.046}, {78, 313.146}, {79, 318.261}, {80, 323.39}, {81, 328.533}, {82, 333.691}, {83, 338.862}, {84, 344.047}, {85, 349.246}, {86, 354.458}, {87, 359.683}, {88, 364.922}, {89, 370.173}, {90, 375.437}, {91, 380.714}, {92, 386.003}, {93, 391.305}, {94, 396.619}, {95, 401.945}, {96, 407.283}, {97, 412.634}, {98, 417.995}, {99, 423.369}, {100, 428.754}}
If $r \approx n \log n$ is correct, and I choose $(1,1,1,\ldots, 1, 1/(n \log n))$ as my representative, then I get
$$A \approx 1 - 1/n, \ G \approx 1-\log n /n,\ \ H \approx 1/\log n \
\mbox{and} \ \theta \approx 1-\log n/n.$$
But finding explicit inequalities that achieve this behavior doesn't sound fun to me.