You could proceed from $(**)$ as follows: for $x \leq x_0$,
\begin{align*}
&\mathbb P(X_1 - X_{(1)} > x_0 - x)\\
&= 1-\mathbb P(X_1 - X_{(1)} \leq x_0 - x)\\
&= 1-\mathbb P(X_{(1)} \geq X_1 - x_0 + x)\\
&= 1-\mathbb P(X_i \geq X_1 - x_0 + x\text{ for all $i=1,\dots,n$})\\
&= 1-\mathbb P(X_i \geq X_1 - x_0 + x\text{ for all $i=2,\dots,n$})\\
&= 1-\int_{\theta}^\infty \mathbb P(X_i \geq t - x_0 + x\text{ for all $i=2,\dots,n$})f_{X_1}(t)\ dt\\
&= 1- \int_\theta^\infty f_{X_1}(t)\prod_{i=2}^n \mathbb P(X_i \geq t-x_0+x)\ dt \\
&= \dots
\end{align*}
However, there appears to be a problem which needs to be addressed first. When referring to a "complete sufficient statistic for $\theta$", we need to carefully consider what space of parameters $\Theta$ is being implicitly referenced. For the space $\Theta = [0,\infty)$, it is true that $X_{(1)}$ is a complete sufficient statistic, and in this case the identification of $\frac{n-1}ne^{X_{(1)}-x_0}$ as the UMVUE for $e^{-x_0+\theta}$ is correct (for $n>1$). However, this does not provide a UMVUE for $f_X(x_0)$ on this space $\Theta$, as the equation $f_X(x_0) = e^{-x_0+\theta}$ does not hold when $\theta>x_0$. The statistic $T(\mathbf X)$ also fails to be unbiased when $\theta>x_0$.
You may instead be intending to consider the space $\Theta = [0, x_0]$, but in this case $X_{(1)}$ is no longer complete, so you will need a different complete, sufficient statistic ($\min\{X_{(1)}, x_0\}$ should do the trick.) After this is sorted out, you should be able to proceed in a similar way as above.
Update:
Consider the estimator
$$\hat 0 = \bar{X} - cS$$
where $c$ is given in your post. This is is an unbiased estimator of $0$ and will clearly be correlated with the estimator given below (for any value of $a$).
Theorem 6.2.25 from C&B shows how to find complete sufficient statistics for the Exponential family so long as $$\{(w_1(\theta), \cdots w_k(\theta)\}$$ contains an open set in $\mathbb R^k$. Unfortunately this distribution yields $w_1(\theta) = \theta^{-2}$ and $w_2(\theta) = \theta^{-1}$ which does NOT form an open set in $R^2$ (since $w_1(\theta) = w_2(\theta)^2$). It is because of this that the statistic $(\bar{X}, S^2)$ is not complete for $\theta$, and it is for the same reason that we can construct an unbiased estimator of $0$ that will be correlated with any unbiased estimator of $\theta$ that is based on the sufficient statistics.
Another Update:
From here, the argument is constructive. It must be the case that there exists another unbiased estimator $\tilde\theta$ such that $Var(\tilde\theta) < Var(\hat\theta)$ for at least one $\theta \in \Theta$.
Proof: Let suppose that $E(\hat\theta) = \theta$, $E(\hat 0) = 0$ and $Cov(\hat\theta, \hat 0) < 0$ (for some value of $\theta$). Consider a new estimator
$$\tilde\theta = \hat\theta + b\hat0$$
This estimator is clearly unbiased with variance
$$Var(\tilde\theta) = Var(\hat\theta) + b^2Var(\hat0) + 2bCov(\hat\theta,\hat0)$$
Let $M(\theta) = \frac{-2Cov(\hat\theta, \hat0)}{Var(\hat0)}$.
By assumption, there must exist a $\theta_0$ such that $M(\theta_0) > 0$. If we choose $b \in (0, M(\theta_0))$, then $Var(\tilde\theta) < Var(\hat\theta)$ at $\theta_0$. Therefore $\hat\theta$ cannot be the UMVUE. $\quad \square$
In summary: The fact that $\hat\theta$ is correlated with $\hat0$ (for any choice of $a$) implies that we can construct a new estimator which is better than $\hat\theta$ for at least one point $\theta_0$, violating the uniformity of $\hat\theta$ claim for best unbiasedness.
Let's look at your idea of linear combinations more closely.
$$\hat\theta = a \bar X + (1-a)cS$$
As you point out, $\hat\theta$ is a reasonable estimator since it is based on Sufficient (albeit not complete) statistics. Clearly, this estimator is unbiased, so to compute the MSE we need only compute the variance.
\begin{align*}
MSE(\hat\theta) &= a^2 Var(\bar{X}) + (1-a)^2 c^2 Var(S) \\
&= \frac{a^2\theta^2}{n} + (1-a)^2 c^2 \left[E(S^2) - E(S)^2\right] \\
&= \frac{a^2\theta^2}{n} + (1-a)^2 c^2 \left[\theta^2 - \theta^2/c^2\right] \\
&= \theta^2\left[\frac{a^2}{n} + (1-a)^2(c^2 - 1)\right]
\end{align*}
By differentiating, we can find the "optimal $a$" for a given sample size $n$.
$$a_{opt}(n) = \frac{c^2 - 1}{1/n + c^2 - 1}$$
where
$$c^2 = \frac{n-1}{2}\left(\frac{\Gamma((n-1)/2)}{\Gamma(n/2)}\right)^2$$
A plot of this optimal choice of $a$ is given below.
It is somewhat interesting to note that as $n\rightarrow \infty$, we have $a_{opt}\rightarrow \frac{1}{3}$ (confirmed via Wolframalpha).
While there is no guarantee that this is the UMVUE, this estimator is the minimum variance estimator of all unbiased linear combinations of the sufficient statistics.
Best Answer
Since $X_j$ is the number of failures preceding the first success for each $j$, $T=\sum\limits_{j=1}^n X_j$ is the number of failures before the $n$th success. Therefore pmf of $T$ is
$$P(T=t)=\binom{n+t-1}{t}\theta^n(1-\theta)^{t}\,\mathbf1_{t\in\{0,1,2,\ldots\}}$$
Now,
\begin{align} E\left[W\mid T=t\right]&=\frac{P\left[X_1=0,\sum\limits_{i=2}^n X_i=t\right]}{P(T=t)} \\&=\frac{P(X_1=0)P\left[\sum\limits_{j=1}^{n-1} X_j=t\right]}{P(T=t)}\qquad,\,\small j=i-1 \end{align}
So the '$t$' remains as it is; it is a matter of '$n$' and '$n-1$' in the pmf of negative binomial.