Your present working appears to be taking all the $x_i$ values to be the same, which is not sufficient generality to properly describe your problem. In your initial setup for the problem there is no stipulation that these values must be the same. Taking the values $\mathbf{x} = (x_1,...,x_n)$ to be fixed explanatory variables, your response variables are independent (but not identically distributed) geometric random variables:
$$y_i | \mathbf{x} \sim \text{Geom}(\theta_i)
\quad \quad \quad
\theta_i \equiv \frac{\beta x_i}{1 + \beta x_i}.$$
As you can see, this generality makes the problem much more complicated than how you are treating it. If the $x_i$ values are different to each other then the response variables $y_i$ are not identically distributed, which will give you a more complicated likelihood function than what you are using. As shown below, this leads to an MLE that is defined (in general) only by an implicit function.
Finding the MLE: To get you back on track, I will show you how to derive the implicit equation for the MLE. For the observed data $\mathbf{y} = (y_1,...,y_n)$ we have the log-likelihood function:
$$\begin{aligned}
\ell_{\mathbf{y}}(\beta)
&= \sum_{i=1}^n \ln \text{Geom} \bigg( y_i \Bigg| \frac{\beta x_i}{1 + \beta x_i} \bigg) \\[6pt]
&= \sum_{i=1}^n \Big[ \ln(\beta) + \ln(x_i) - y_i \ln(1+\beta x_i) \Big] \\[6pt]
&= n \ln(\beta) + \sum_{i=1}^n \ln(x_i) - \sum_{i=1}^n y_i \ln(1+\beta x_i). \\[6pt]
\end{aligned}$$
So you have the score function and Hessian function:
$$\begin{aligned}
s_\mathbf{y}(\beta)
\equiv \frac{\partial \ell_{\mathbf{y}}}{\partial \beta} (\beta)
&= \frac{n}{\beta} - \sum_{i=1}^n \frac{x_i y_i}{1+\beta x_i}, \\[6pt]
H_\mathbf{y}(\beta)
\equiv \frac{\partial^2 \ell_{\mathbf{y}}}{\partial \beta^2} (\beta)
&= - \frac{n}{\beta^2} + \sum_{i=1}^n \frac{x_i^2 y_i}{(1+\beta x_i)^2}. \\[6pt]
\end{aligned}$$
At any critical point $\hat{\beta}$ (satisfying the critical point equation $s_\mathbf{y}(\hat{\beta}) =0$) you have:
$$\begin{aligned}
H_\mathbf{y}(\hat{\beta})
&= - \frac{n}{\hat{\beta}^2} + \sum_{i=1}^n \frac{x_i^2 y_i}{(1+\hat{\beta} x_i)^2} \\[6pt]
&= - \Bigg[ \frac{1}{n} \bigg( \sum_{i=1}^n \frac{x_i y_i}{1+\hat{\beta} x_i} \bigg)^2 - \sum_{i=1}^n \frac{x_i^2 y_i}{(1+\hat{\beta} x_i)^2} \Bigg] \\[6pt]
&= - \Bigg[ \frac{1}{n} \bigg( \sum_{i=1}^n \sum_{j=1}^n \frac{x_i y_i}{1+\hat{\beta} x_i} \frac{x_j y_j}{1+\hat{\beta} x_j} \bigg) - \sum_{i=1}^n \frac{x_i^2 y_i}{(1+\hat{\beta} x_i)^2} \Bigg] \\[6pt]
&= - \sum_{i=1}^n \frac{x_i y_i}{1+\hat{\beta} x_i} \bigg( \frac{1}{n} \sum_{j=1}^n \frac{x_j y_j}{1+\hat{\beta} x_j} - \frac{x_i}{1+\hat{\beta} x_i} \bigg). \\[6pt]
\end{aligned}$$
The maximising point is defined implicitly by the critical point equation:
$$\frac{n}{\hat{\beta}} = \sum_{i=1}^n \frac{x_i y_i}{1+\hat{\beta} x_i}.$$
In the special case where $x \equiv x_1 = \cdots = x_n$ you have $\hat{\beta} = 1/(x (\bar{y}_n-1))$, but this explicit result does not hold in general.
Finding the asymptotic distribution of the MLE: If you want to find the asymptotic variance of the MLE, there are a few ways to do it. The complicated way is to differentiate the implicit function multiple times to get a Taylor approximation to the MLE, and then use this to get an asymptotic result for the variance of the MLE. That is quite a bit of work, and I will leave you to do it if you want to, but the above should get you started by giving you the correct equation for the MLE. You can read more about dealing with implicitly defined random variables using Taylor series and the delta method in Benichou and Gail (1989).
The simpler way to get the MLE is to rely on asymptotic theory for MLEs. Under some regularity conditions, you have the asymptotic distribution:
$$\sqrt{n}(\hat{\beta} - \beta)\overset{\rightarrow}{\sim} \text{N} \bigg( 0, \frac{1}{\mathcal{I}(\beta)} \bigg),$$
where $\mathcal{I}$ is the expected Fisher information for a single observation. Since $\mathbb{E}(\theta_i Y_i) = 1$ and $\mathbb{E}(\theta_i^2 Y_i^2) = 2-\theta_i$, in this case we have the Fisher information function for $n$ observations as being:
$$\begin{aligned}
n\mathcal{I}(\beta)
&= \mathbb{E}(s_\mathbf{Y}(\beta)^2) \\[6pt]
&= \mathbb{E} \Bigg( \bigg( \frac{n}{\beta} - \sum_{i=1}^n \frac{x_i Y_i}{1+\beta x_i} \bigg)^2 \Bigg) \\[6pt]
&= \frac{n^2}{\beta^2} \cdot \mathbb{E} \Bigg( \bigg( 1 - \frac{1}{n} \sum_{i=1}^n \theta_i Y_i \bigg)^2 \Bigg) \\[6pt]
&= \frac{n^2}{\beta^2} \cdot \mathbb{E} \Bigg( 1 - \frac{2}{n} \sum_{i=1}^n \theta_i Y_i + \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n \theta_i \theta_j Y_i Y_j \Bigg) \\[6pt]
&= \frac{n^2}{\beta^2} \cdot \mathbb{E} \Bigg( 1 - \frac{2}{n} \sum_{i=1}^n \theta_i Y_i + \frac{1}{n^2} \sum_{i=1}^n \theta_i^2 Y_i^2 + \frac{1}{n^2} \sum_{i \neq j} \theta_i \theta_j Y_i Y_j \Bigg) \\[6pt]
&= \frac{n^2}{\beta^2} \cdot \Bigg( 1 - 2 + \frac{1}{n^2} \sum_{i=1}^n (2-\theta_i) + \frac{n-1}{n} \Bigg) \\[6pt]
&= \frac{n^2}{\beta^2} \cdot \Bigg( 1 - 2 + \frac{2}{n} - \frac{1}{n^2} \sum_{i=1}^n \theta_i + 1 - \frac{1}{n} \Bigg) \\[6pt]
&= \frac{n^2}{\beta^2} \cdot \Bigg( \frac{1}{n} - \frac{1}{n^2} \sum_{i=1}^n \theta_i \Bigg) \\[6pt]
&= \frac{n(1 - \bar{\theta}_n)}{\beta^2}. \\[6pt]
\end{aligned}$$
Here is an alternative derivation using the regularity conditions:
$$\begin{aligned}
n\mathcal{I}(\beta)
&= -\mathbb{E}(H_\mathbf{Y}(\beta)) \\[6pt]
&= \mathbb{E} \Bigg( \frac{n}{\beta^2} - \sum_{i=1}^n \frac{x_i^2 Y_i}{(1+\beta x_i)^2} \Bigg) \\[6pt]
&= \frac{n}{\beta^2} \cdot \mathbb{E} \Bigg( 1 - \frac{1}{n} \sum_{i=1}^n \theta_i^2 Y_i \Bigg) \\[6pt]
&= \frac{n}{\beta^2} \cdot \Bigg( 1 - \frac{1}{n} \sum_{i=1}^n \theta_i \Bigg) \\[6pt]
&= \frac{n(1-\bar{\theta}_n)}{\beta^2}. \\[6pt]
\end{aligned}$$
So in this case you have the asymptotic variance:
$$\frac{1}{\mathcal{I}(\beta)} = \frac{\beta^2}{(1 - \bar{\theta}_\infty)}$$
where $\bar\theta_\infty = \lim_{n\rightarrow\infty} \bar\theta_n$. Based on this, for large but finite sample sizes, a large-sample variance for your estimator $\hat\beta$ is given by
$$\mathbb{V}(\hat{\beta}) \approx \frac{\hat\beta^2}{n(1 - \bar{\theta}_n)}$$
