Gradient of the log likelihood for energy based models Coming from this recent paper Your Classifier is Secretly an Energy Based Model And You Should Treat it Like One, they give the following definition...
$$
p_\theta(\mathbf{x}) = \frac{\exp(-E_\theta(\mathbf{x}))}{Z(\theta)}
$$
with $Z(\theta) = \int_x \exp(-E_\theta(x))dx$, then consequently it says that it can be expressed as
$$
\frac{\partial \log p_\theta(x)}{\partial\theta} = \mathbb{E}_{p_\theta(x')}\Big[\frac{\partial E_\theta(x')}{\partial \theta} \Big] - \frac{\partial E_\theta(x)}{\partial \theta}
$$
I tried getting there myself but I am not sure if this is correct or not
$$
\begin{aligned}
\nabla_\theta \log p_\theta(x) &= \frac{1}{p_\theta(x)} \nabla_\theta p_\theta(x) \\
&= \frac{Z(\theta)}{\exp(-E_\theta(x))} \Big(-\exp(-E_\theta(x))\nabla_\theta E_\theta(x)Z(\theta)^{-1} + \nabla_\theta Z(\theta)^{-1} \exp(-E_\theta(x)) \Big) \\
&= -\nabla_\theta E_\theta(x) + Z(\theta)\nabla_\theta Z(\theta)^{-1} \\
&= -\nabla_\theta E_\theta(x) + Z(\theta)\nabla_\theta \int_x \exp(E_\theta(x)) \\
&= Z(\theta) \int_x \exp(E_\theta(x)) \nabla_\theta E_\theta(x) -\nabla_\theta E_\theta(x) \\
\end{aligned}
$$
The second step comes from the product rule, and everything is just simplified from there. The sign inside the exponential and the exponent of $Z(\theta)$ are the opposite of what they should be. (If it were $Z(\theta)^{-1}$ and $\exp(-E_\theta(x))$ then I would have the expectation perfectly)
I must have gone wrong somewhere, but I do not see what I did. Can anyone spot the mistake?

The issue emerges in the evaluation of the second term in line $(3)$ and $(4)$ of your derivation. Note that
$$\nabla_{\theta} Z(\theta)^{-1} = \nabla_{\theta} \frac{1}{\int_x \exp(-E_{\theta}(x))\, dx} \neq \nabla_{\theta}\int_x \exp(E(_{\theta}(x)) \, dx.$$
Instead, we have
$$\nabla_{\theta} Z(\theta)^{-1} = \nabla_{\theta} \left( \frac{1}{Z(\theta)} \right) = \frac{-\nabla_{\theta}Z(\theta)}{Z(\theta)^2}.$$
Hence,
\begin{align}Z(\theta)\nabla_{\theta} Z(\theta)^{-1} &= \frac{-1}{Z(\theta)}\nabla_{\theta} Z(\theta) \\
&= \frac{-1}{Z(\theta)} \nabla_{\theta} \int_x \exp(-E_{\theta}(x)) \, dx \\
&= \frac{-1}{Z(\theta)} \int_x \nabla_{\theta}  \exp(-E_{\theta}(x)) \, dx \\
&= \frac{1}{Z(\theta)} \int_x \exp(-E_{\theta}(x)) \nabla_{\theta} E_{\theta}(x) \, dx, ​
\end{align}
which yields the expression you seek.
There are some technicalities concerning differentiating under the integral, or commutativity of $\nabla_{\theta}$ and $\int_x$ operators, on which this derivation relies. If the functional form of $E_{\theta}(x)$ is such that $p_{\theta}(x)$ parametrises an exponential family, then there are no issues. If not, then depending on whether the support of $p_{\theta}$ and hence limits of integration in $\int_x$ are finite or infinite, then you will need Leibniz integral rule or Lebesgue's dominated convergence theorem to justify this. See Statistical Inference by Casella and Berger (2004) if you are looking for simple tests of the latter without recourse to analysis results pertaining to Lesbesgue integration.
