Yes, $f(x\mid\theta)$ is an even function of $x$, but how can the pdf 'change' completely in your calculations? You carried out calculations for an exponential density whereas you are given a Laplace distribution.
Indeed, the sample is drawn from a Laplace distribution with scale parameter $1/\theta$ and location parameter zero. We can conclude that the joint density $f_{\theta}$ belongs to the one-parameter exponential family. Using this, we can find a complete sufficient statistic for $\theta$.
Due to independence, joint density of $\mathbf X=(X_1,X_2,\cdots,X_n)$ is
\begin{align}
f_{\theta}(\mathbf x)&=\prod_{i=1}^n\frac{\theta}{2}e^{-\theta\,|x_i|}
\\&=\left(\frac{\theta}{2}\right)^ne^{-\theta\sum_{i=1}^n|x_i|}
\\&=\exp\left[-\theta\sum_{i=1}^n|x_i|+n\ln\left(\frac{\theta}{2}\right)\right]\quad,\,\mathbf x\in\mathbb R^n\,,\,\theta>0
\end{align}
Clearly, a complete sufficient statistic for the family of distributions $\{f_{\theta}:\theta>0\}$ is $$T(\mathbf X)=\sum_{i=1}^n|X_i|$$
By Lehmann-Scheffe theorem, an unbiased estimator of $\theta$ based on $T$ will be the UMVUE of $\theta$.
It is a simple exercise to show that $|X_i|\sim\mathcal{Exp}(\theta)$ independently for each $i$, where $\theta$ denotes rate of the distribution. As such, we have $T\sim\mathcal{Ga}(\theta,n)$ with density $$g(t)=\frac{\theta^n e^{-\theta t}t^{n-1}}{\Gamma(n)}\mathbf1_{t>0}$$
We find that
\begin{align}
E\left(\frac{1}{T}\right)&=\int_0^{\infty}\frac{1}{t}\frac{\theta^n e^{-\theta t}t^{n-1}}{\Gamma(n)}\,dt
\\&=\frac{\theta^n}{\Gamma(n)}\frac{\Gamma(n-1)}{\theta^{n-1}}
\\&=\frac{\theta}{n-1}
\end{align}
So the UMVUE of $\theta$ is $$\frac{n-1}{T}=\frac{n-1}{\sum_{i=1}^n|X_i|}$$
It turns out that both approaches (my initial attempt and another based on suggestions in the comment section) in my original post give the same answer. I will outline both methods here for a complete answer to the question.
Here, $\text{Gamma}(n,\theta)$ means the gamma density $f(y)=\frac{\theta^n}{\Gamma(n)}e^{-\theta y}y^{n-1}\mathbf1_{y>0}$ where $\theta,n>0$, and $\text{Exp}(\theta)$ denotes an exponential distribution with mean $1/\theta$, ($\theta>0$). Clearly, $\text{Exp}(\theta)\equiv\text{Gamma}(1,\theta)$.
Since $T=\sum_{i=1}^n\ln X_i$ is complete sufficient for $\theta$ and $\mathbb E(X_1)=\frac{\theta}{1+\theta}$, by the Lehmann-Scheffe theorem $\mathbb E(X_1\mid T)$ is the UMVUE of $\frac{\theta}{1+\theta}$.
So we have to find this conditional expectation.
We note that $X_i\stackrel{\text{i.i.d}}{\sim}\text{Beta}(\theta,1)\implies-\ln X_i\stackrel{\text{i.i.d}}{\sim}\text{Exp}(\theta)\implies-T\sim\text{Gamma}(n,\theta)$.
Method I:
Let $U=-\ln X_1$ and $V=-\sum_{i=2}^n\ln X_i$, so that $U$ and $V$ are independent. Indeed, $U\sim\text{Exp}(\theta)$ and $V\sim\text{Gamma}(n-1,\theta)$, implying $U+V\sim\text{Gamma}(n,\theta)$.
So, $\mathbb E(X_1\mid \sum_{i=1}^n\ln X_i=t)=\mathbb E(e^{-U}\mid U+V=-t)$.
Now we find the conditional distribution of $U\mid U+V$.
For $n>1$ and $\theta>0$, joint density of $(U,V)$ is
\begin{align}f_{U,V}(u,v)&=\theta e^{-\theta u}\mathbf1_{u>0}\frac{\theta^{n-1}}{\Gamma(n-1)}e^{-\theta v}v^{n-2}\mathbf1_{v>0}\\&=\frac{\theta^n}{\Gamma(n-1)}e^{-\theta(u+v)}v^{n-2}\mathbf1_{u,v>0}\end{align}
Changing variables, it is immediate that the joint density of $(U,U+V)$ is $$f_{U,U+V}(u,z)=\frac{\theta^n}{\Gamma(n-1)}e^{-\theta z}(z-u)^{n-2}\mathbf1_{0<u<z}$$
Let $f_{U+V}(\cdot)$ be the density of $U+V$. Thus conditional density of $U\mid U+V=z$ is \begin{align}f_{U\mid U+V}(u\mid z)&=\frac{f_{U,U+V}(u,z)}{f_{U+V}(z)}\\&=\frac{(n-1)(z-u)^{n-2}}{z^{n-1}}\mathbf1_{0<u<z}\end{align}
Therefore, $\displaystyle\mathbb E(e^{-U}\mid U+V=z)=\frac{n-1}{z^{n-1}}\int_0^z e^{-u}(z-u)^{n-2}\,\mathrm{d}u$.
That is, the UMVUE of $\frac{\theta}{1+\theta}$ is $\displaystyle\mathbb E(X_1\mid T)=\frac{n-1}{(-T)^{n-1}}\int_0^{-T} e^{-u}(-T-u)^{n-2}\,\mathrm{d}u\tag{1}$
Method II:
As $T$ is a complete sufficient statistic for $\theta$, any unbiased estimator of $\frac{\theta}{1+\theta}$ which is a function of $T$ will be the UMVUE of $\frac{\theta}{1+\theta}$ by the Lehmann-Scheffe theorem. So we proceed to find the moments of $-T$, whose distribution is known to us. We have,
$$\mathbb E(-T)^r=\int_0^\infty y^r\theta^n\frac{e^{-\theta y}y^{n-1}}{\Gamma(n)}\,\mathrm{d}y=\frac{\Gamma(n+r)}{\theta^r\,\Gamma(n)},\qquad n+r>0$$
Using this equation we obtain unbiased estimators (and UMVUE's) of $1/\theta^r$ for every integer $r\ge1$.
Now for $\theta>1$, we have $\displaystyle\frac{\theta}{1+\theta}=\left(1+\frac{1}{\theta}\right)^{-1}=1-\frac{1}{\theta}+\frac{1}{\theta^2}-\frac{1}{\theta^3}+\cdots$
Combining the unbiased estimators of $1/\theta^r$ we obtain $$\mathbb E\left(1+\frac{T}{n}+\frac{T^2}{n(n+1)}+\frac{T^3}{n(n+1)(n+2)}+\cdots\right)=1-\frac{1}{\theta}+\frac{1}{\theta^2}-\frac{1}{\theta^3}+\cdots$$
That is, $$\mathbb E\left(\sum_{r=0}^\infty \frac{T^r}{n(n+1)...(n+r-1)}\right)=\frac{\theta}{1+\theta}$$
So assuming $\theta>1$, the UMVUE of $\frac{\theta}{1+\theta}$ is $g(T)=\displaystyle\sum_{r=0}^\infty \frac{T^r}{n(n+1)...(n+r-1)}\tag{2}$
I am not certain about the case $0<\theta<1$ in the second method.
According to Mathematica, equation $(1)$ has a closed form in terms of the incomplete gamma function. And in equation $(2)$, we can express the product $n(n+1)(n+2)...(n+r-1)$ in terms of the usual gamma function as $n(n+1)(n+2)...(n+r-1)=\frac{\Gamma(n+r)}{\Gamma(n)}$. This perhaps provides the apparent connection between $(1)$ and $(2)$.
Using Mathematica I could verify that $(1)$ and $(2)$ are indeed the same thing.
Best Answer
The Poisson distribution is a one-parameter exponential family distribution, with natural sufficient statistic given by the sample total $T(\mathbf{x}) = \sum_{i=1}^n x_i$. The canonical form is:
$$p(\mathbf{x}|\theta) = \exp \Big( \ln (\theta) T(\mathbf{x}) - n\theta \Big) \cdot h(\mathbf{x}) \quad \quad \quad h(\mathbf{x}) = \coprod_{i=1}^n x_i! $$
From this form it is easy to establish that $T$ is a complete sufficient statistic for the parameter $\theta$. So the Lehmann–Scheffé theorem means that for any $g(\theta)$ there is only one unbiased estimator of this quantity that is a function of $T$, and this is the is UMVUE of $g(\theta)$. One way to find this estimator (the method you are using) is via the Rao-Blackwell theorem --- start with an arbitrary unbiased estimator of $g(\theta)$ and then condition on the complete sufficient statistic to get the unique unbiased estimator that is a function of $T$.
Using Rao-Blackwell to find the UMVUE: In your case you want to find the UMVUE of:
$$g(\theta) \equiv \theta \exp (-\theta).$$
Using the initial estimator $\hat{g}_*(\mathbf{X}) \equiv \mathbb{I}(X_1=1)$ you can confirm that,
$$\mathbb{E}(\hat{g}_*(\mathbf{X})) = \mathbb{E}(\mathbb{I}(X_1=1)) = \mathbb{P}(X_1=1) = \theta \exp(-\theta) = g(\theta),$$
so this is indeed an unbiased estimator. Hence, the unique UMVUE obtained from the Rao-Blackwell technique is:
$$\begin{equation} \begin{aligned} \hat{g}(\mathbf{X}) &\equiv \mathbb{E}(\mathbb{I}(X_1=1) | T(\mathbf{X}) = t) \\[6pt] &= \mathbb{P}(X_1=1 | T(\mathbf{X}) = t) \\[6pt] &= \mathbb{P} \Big( X_1=1 \Big| \sum_{i=1}^n X_i = t \Big) \\[6pt] &= \frac{\mathbb{P} \Big( X_1=1 \Big) \mathbb{P} \Big( \sum_{i=2}^n X_i = t-1 \Big)}{\mathbb{P} \Big( \sum_{i=1}^n X_i = t \Big)} \\[6pt] &= \frac{\text{Pois}(1| \theta) \cdot \text{Pois}(t-1| (n-1)\theta)}{\text{Pois}(t| n\theta)} \\[6pt] &= \frac{t!}{(t-1)!} \cdot \frac{ \theta \exp(-\theta) \cdot ((n-1) \theta)^{t-1} \exp(-(n-1)\theta)}{(n \theta)^t \exp(-n\theta)} \\[6pt] &= t \cdot \frac{ (n-1)^{t-1}}{n^t} \\[6pt] &= \frac{t}{n} \Big( 1- \frac{1}{n} \Big)^{t-1} \\[6pt] \end{aligned} \end{equation}$$
Your answer has a slight error where you have conflated the sample mean and the sample total, but most of your working is correct. As $n \rightarrow \infty$ we have $(1-\tfrac{1}{n})^n \rightarrow \exp(-1)$ and $t/n \rightarrow \theta$, so taking these asymptotic results together we can also confirm consistency of the estimator:
$$\hat{g}(\mathbf{X}) = \frac{t}{n} \Big[ \Big( 1- \frac{1}{n} \Big)^n \Big] ^{\frac{t}{n} - \frac{1}{n}} \rightarrow \theta [ \exp (-1) ]^\theta = \theta \exp (-\theta) = g(\theta).$$
This latter demonstration is heuristic, but it gives a nice check on the working. It is interesting here that you get an estimator that is a finite approximation to the exponential function of interest.