If one is familiar with the concepts of sufficiency and completeness, then this problem is not too difficult. Note that $f(x; \theta)$ is the density of a $\Gamma(2, \theta)$ random variable. The gamma distribution falls within the class of the exponential family of distributions, which provides rich statements regarding the construction of uniformly minimum variance unbiased estimators via notions of sufficiency and completeness.
The distribution of a random sample of size $n$ from this distribution is
$$
g(x_1,\ldots,x_n; \theta) = \theta^{2n} \exp\Big(-\theta \sum_{i=1}^n x_i + \sum_{i=1}^n \log x_i\Big)
$$
which, again, conforms to the exponential family class.
From this we can conclude that $S_n = \sum_{i=1}^n X_i$ is a complete, sufficient statistic for $\theta$. Operationally, this means that if we can find some function $h(S_n)$ that is unbiased for $\theta$, then we know immediately via the Lehmann-Scheffe theorem that $h(S_n)$ is the unique uniformly minimum variance unbiased (UMVU) estimator.
Now, $S_n$ has distribution $\Gamma(2n, \theta)$ by standard properties of the gamma distribution. (This can be easily checked via the moment-generating function.)
Furthermore, straightforward calculus shows that
$$
\mathbb{E} S_n^{-1} = \int_0^\infty s^{-1} \frac{\theta^{2n} s^{2n - 1}e^{-\theta s}}{\Gamma(2n)} \,\mathrm{d}s = \frac{\theta}{2n - 1} \>.
$$
Hence, $h(S_n) = \frac{2n-1}{S_n}$ is unbiased for $\theta$ and must, therefore, be the UMVU estimator.
Addendum: Using the fact that $\newcommand{\e}{\mathbb{E}}\e S_n^{-2} = \frac{\theta^2}{(2n-1)(2n-2)}$, we conclude that the $\mathbb{V}ar(h(S_n)) = \frac{\theta^2}{2(n-1)}$. On the other hand, the information $I(\theta)$ from a sample of size one is readily computed to be $-\e \frac{\partial^2 \log f}{\partial \theta^2} = 2 \theta^{-2}$ and so the Cramer-Rao lower bound for a sample of size $n$ is
$$
\mathrm{CRLB}(\theta) = \frac{1}{n I(\theta)} = \frac{\theta^2}{2n} \> .
$$
Hence, $h(S_n)$ does not achieve the bound, though it comes close, and indeed, achieves it asymptotically.
However, if we reparametrize the density by taking $\beta = \theta^{-1}$ so that
$$
f(x;\beta) = \beta^{-2} x e^{-x/\beta},\quad x > 0,
$$
then the UMVU estimator for $\beta$ can be shown to be $\tilde{h}(S_n) = \frac{S_n}{2 n}$. (Just check that it's unbiased!) The variance of this estimator is $\mathbb{V}ar(\tilde{h}(S_n)) = \frac{\beta^2}{2n}$ and this coincides with the CRLB for $\beta$.
The point of the addendum is that the ability to achieve (or not) the CRLB depends on the particular parametrization used and even when there is a one-to-one correspondence between two unique parametrizations, an unbiased estimator for one may achieve the Cramer-Rao lower bound while the other one does not.
Recall that for every random variable $X$ Poisson with parameter $\lambda$ one has $E(X)=\lambda$ and $E(X^2)=\lambda^2+\lambda$ hence $\mathrm{var}(X)=\lambda$.
Thus, if $(Y_k)$ is i.i.d. Poisson $\lambda$ and $\bar Y=\frac1n\sum\limits_{k=1}^nY_k$ then $\bar Y$ is an unbiased estimator of $\lambda$ since $E(\bar Y)=\lambda$.
Likewise, $\mathrm{var}(\bar Y)=\frac1{n}\mathrm{var}(Y_1)=\frac1n\lambda$ hence $E(\bar Y^2)=\frac1n\lambda+\lambda^2$. Solving for $\theta=3\lambda+\lambda^2$ yields $\theta=E(\bar Y^2)+(3-\frac1n)E(\bar Y)$ hence an unbiased estimator of $\theta$ is
$$
\Theta=\bar Y^2+\left(3-\frac1n\right)\bar Y.
$$
Best Answer
Notice that $P(X_n = 0) = e^{-\lambda}$. This motivates the following estimator. Set $Y_n = 1(X_n = 0)$, so that $Y_n \sim \mathrm{Bern}(e^{- \lambda})$ and $\mathbb{E}(Y_n) = P(X_n = 0) = e^{-\lambda}$. The estimator $$\overline{Y}_n = \frac{1}{n}(Y_1 + \ldots + Y_n) = \frac{1}{n} (1(X_1 = 0) + \ldots + 1(X_n = 0))$$ is then an unbiased estimator for $\theta = e^{- \lambda}$.