The following is just a minor variation of Martin Argerami's proof of the old question. I am even copying his equations and some of his text. If you are +1ing this post, please also +1 his one (if not already done).
Here is a counterexample for $n=2$. Let $\varepsilon $ be a positive real $<
\dfrac{1}{2}$. The matrix
\[
a=
\begin{bmatrix}
1 & 0 & 0 & 1-\varepsilon \\
0 & \varepsilon & 0 & 0 \\
0 & 0 & \varepsilon & 0 \\
1-\varepsilon & 0 & 0 & 1
\end{bmatrix}
\in \mathrm{M}_{4}\left( \mathbb{C}\right)
\]
is positive-definite, but it cannot be written as a sum of tensor products
of nonnegative-semidefinite $2\times 2$-matrices. Here is why:
Assume the contrary. Thus, $a$ is written in the form
$a=\sum_{j} \left[ \begin{matrix}
\alpha _{j} & \overline{\gamma _{j}} \\
\gamma _{j} & \beta _{j}
\end{matrix} \right]
\otimes \left[
\begin{matrix}
\alpha _{j}^{\prime } & \overline{\gamma _{j}^{\prime }} \\
\gamma _{j}^{\prime } & \beta _{j}^{\prime }
\end{matrix} \right]
=\left[
\begin{matrix}
\sum_{j}\alpha _{j}^{\prime }\alpha _{j} & \sum_{j}\alpha _{j}^{\prime }
\overline{\gamma _{j}} & \sum_{j}\overline{\gamma _{j}^{\prime }}\alpha _{j}
& \sum_{j}\overline{\gamma _{j}^{\prime }}\gamma _{j} \\
\sum_{j}\alpha _{j}^{\prime }\gamma _{j} & \sum_{j}\alpha _{j}^{\prime
}\beta _{j} & \ast & \ast \\
\ast & \ast & \sum_{j}\beta _{j}^{\prime }\alpha _{j} & \ast \\
\ast & \ast & \ast & \ast
\end{matrix}\right] $,
where $j$ ranges from $1$ to some positive integer $N$. Since each $
\begin{bmatrix}
\alpha _{j} & \overline{\gamma _{j}} \\
\gamma _{j} & \beta _{j}
\end{bmatrix}
$ is nonnegative-semidefinite, we have $\alpha _{j}\geq 0$, $\beta _{j}\geq 0$,
and $\alpha _{j}\beta _{j}\geq \left\vert \gamma _{j}\right\vert ^{2}$ for
all $j$. Similarly, $\alpha _{j}^{\prime }\geq 0$, $\beta _{j}^{\prime }\geq
0$, and $\alpha _{j}^{\prime }\beta _{j}^{\prime }\geq \left\vert \gamma
_{j}^{\prime }\right\vert ^{2}$ for all $j$.
Now, comparing the entries of $a$ in this equation, we get $\varepsilon
=\sum_{j}\alpha _{j}^{\prime }\beta _{j}$ (from the $\left( 2,2\right) $-th
entry) and $\varepsilon =\sum_{j}\beta _{j}^{\prime }\alpha _{j}$ (from the $
\left( 3,3\right) $-th entry). Taking the arithmetic mean of these two
equations, we get
$\varepsilon =\dfrac{1}{2}\left( \sum_{j}\alpha
_{j}^{\prime }\beta _{j}+\sum_{j}\beta _{j}^{\prime }\alpha _{j}\right)
=\sum_{j}\dfrac{\alpha _{j}^{\prime }\beta _{j}+\beta _{j}^{\prime }\alpha
_{j}}{2}\geq \sum_{j}\sqrt{\alpha _{j}^{\prime }\beta _{j}\beta _{j}^{\prime
}\alpha _{j}}$
(by AM-GM, since we are dealing with nonnegative reals). But
for every $j$, we have
$\sqrt{\alpha _{j}^{\prime }\beta _{j}\beta
_{j}^{\prime }\alpha _{j}}=\sqrt{\alpha _{j}\beta _{j}}\sqrt{\alpha
_{j}^{\prime }\beta _{j}^{\prime }}\geq \left\vert \gamma _{j}\right\vert
\left\vert \gamma _{j}^{\prime }\right\vert $
(since $\alpha _{j}\beta_{j}\geq \left\vert \gamma _{j}\right\vert ^{2}$ and $\alpha _{j}^{\prime }\beta _{j}^{\prime }\geq \left\vert \gamma _{j}^{\prime }\right\vert ^{2}$
), so this becomes
$\varepsilon \geq \sum_{j}\underbrace{\left\vert \gamma
_{j}\right\vert \left\vert \gamma _{j}^{\prime }\right\vert }_{=\left\vert
\overline{\gamma _{j}^{\prime }}\gamma _{j}\right\vert } =
\sum_{j}\left\vert \overline{\gamma _{j}^{\prime }}\gamma _{j}\right\vert
\geq \left\vert \sum_{j}\overline{\gamma _{j}^{\prime }}\gamma
_{j}\right\vert $
(by the triangle inequality). But since $1-\varepsilon
=\sum_{j}\overline{\gamma _{j}^{\prime }}\gamma _{j}$ (by comparing the $
\left( 1,4\right) $-th entry of the matrices in the above equation), this
becomes $\varepsilon \geq \left\vert 1-\varepsilon \right\vert $, what
contradicts the definition of $\varepsilon $.
