Although the OP did not respond, I am answering this to showcase the method I proposed (and indicate what statistical intuition it may contain).
First, it is important to distinguish on which entity is the constraint imposed. In a deterministic optimization setting, there is no such issue : there is no "true value", and an estimator of it. We just have to find the optimizer. But in a stochastic setting, there are conceivably two different cases:
a) "Estimate the parameter given a sample that has been generated by a population that has a non-negative mean" (i.e. $\theta \ge 0$) and
b) "Estimate the parameter under the constraint that your estimator cannot take negative values"(i.e. $\hat \theta \ge 0$).
In the first case, imposing the constraint is including prior knowledge on the unknown parameter. In the second case, the constraint can be seen as reflecting a prior belief on the unknown parameter (or some technical, or "strategic", limitation of the estimator).
The mechanics of the solution are the same, though:The objective function, (the log-likelihood augmented by the non-negativity constraint on $\theta$) is
$$\tilde L(\theta|\mathbf{x})=-\frac n2 \ln(2\pi)-\frac{1}{2}\sum_{i=1}^{n}(x_i-\theta)^2 +\xi\theta,\qquad \xi\ge 0 $$
Given concavity, f.o.c is also sufficient for a global maximum. We have
$$\frac {\partial}{\partial \theta}\tilde L(\theta|\mathbf{x})=\sum_{i=1}^{n}(x_i-\theta) +\xi = 0 \Rightarrow \hat \theta = \bar x+\frac{\xi}{n} $$
1) If the solution lies in an interior point ($\Rightarrow \hat \theta >0$), then $\xi=0$ and so the solution is $\{\hat \theta= \bar x>0,\; \xi^*=0\}$.
2) If the solution lies on the boundary ($\Rightarrow \hat \theta =0$) then we obtain the value of the multiplier at the solution $\xi^* = -n\hat x$, and so the full solution is $\{\hat \theta= 0,\; \xi^*=-n\bar x\}$. But since the multiplier must be non-negative, this necessarily implies that in this case we would have $\bar x\le 0$
(There is nothing special about setting the constraint to zero. If say the constraint was $\theta \ge -2$, then if the solution lied on the boundary, $\hat \theta = -2$, it would imply (in order for the multiplier to have a positive value), that $\bar x \le -2$).
So, if the optimizer is $0$ what are we facing here?
If we are in "constraint type-a", i.e we have been told that the sample comes from a population that it has a non-negative mean, then with $\hat \theta =0$ chances are that the sample may not be representative of this population.
If we are in "constraint type-b", i.e. we had the belief that the population has a non-negative mean, with $\hat \theta =0$ this belief is questioned.
(This is essentially an alternative way to deal with prior beliefs, outside the formal bayesian approach).
Regarding the properties of the estimator, one should carefully distinguish this constrained estimation case, with the case where the true parameter lies on the boundary of the parameter space.
(I am going to fully answer this home study, because it appears that the OP is rather away from any path that could be steered in the right direction by simple hints).
It can be shown that the MLE for $\theta$ is the minimum-order statistic,
$$\hat \theta_n = X_{1:n}$$
In fact, one can derive that $\hat \theta_n$ follows itself a Pareto distribution with scale parameter $\theta$ and (since here the original shape parameter is $1$), with shape parameter $n$, i.e. with distribution function
$$F(\hat \theta_n) = 1- \left(\frac {\hat \theta_n}{\theta}\right)^{-n}$$
So in this case we have available the finite sample distribution of the ML estimator.
Still, the OP is asked to consider asymptotics for the function $\sqrt {n} (\hat \theta _n - \theta)$.
Subtracting $\theta$ does not make the distribution centered at zero or anything like that, since we always have
$$\hat \theta_n = X_{1:n} \geq \min {X}= \theta$$
By this observation alone we know that even if $\sqrt {n} (\hat \theta _n - \theta)$ has a limiting distribution, it won't be the normal.
Now, $(\hat \theta_n - \theta)$ has an exact Lomax distribution, with the same parameters as the distribution of $\hat \theta_n$. The Variance of the Lomax distribution is the same as the variance of the Pareto distribution so in our case
$$\text{Var}[(\hat \theta_n - \theta)] = \frac {\theta^2 n}{(n-1)^2(n-2)}$$
Note that the denominator has leading term $n^3$. So to stabilize this variance as $n$ goes to infinity, so that it neither explodes nor goes to zero, one needs to multiply the variable by $n$, not $\sqrt {n}$.
$$\text{Var}[n(\hat \theta_n - \theta)] = \frac {\theta^2 n^3}{(n-1)^2(n-2)} \to \theta^2$$
So $\sqrt {n} (\hat \theta _n - \theta)$ goes to the constant zero (convergence of the estimator is faster than usual and multiplying by $\sqrt {n}$ "is not enough" to maintain a distribution, we need to multiply by $n$). Standardized by its variance (which goes to zero), it explodes.
Can we obtain the limiting distribution of $Z_n = n(\hat \theta_n - \theta)$?
We have for the distribution function
$$F_{Z_n}(z) = \text{Prob}(Z_n \leq z) = \text{Prob}(n(\hat \theta_n - \theta) \leq z) = \text{Prob}\left(\hat \theta_n \leq \theta + \frac {z}{n} \right)$$
$$ \implies F_{Z_n}(z) = 1- \left(\frac {\theta + \frac {z}{n}}{\theta}\right)^{-n} = 1- \left(1+ \frac {(z/\theta)}{n}\right)^{-n}$$
$$\implies \lim_{n \to \infty}F_{Z_n}(z) = 1 - \exp\{-z/\theta\}$$
which is the CDF of the Exponential distribution, with expected value $\theta$ and variance $\theta^2$.
Best Answer
There is no MLE of $\theta$.
This expression of the likelihood function is incorrect. The left-hand side expresses the likelihood as a function of $\theta$, but the right hand side is a function of $p$.
You can not fix this either because $\theta$ is a non-injective function of $p$ and can not be inverted. Given $\theta$ we have two possible values of $p$
$$p=\frac{1}{2}\pm\sqrt{\frac{1}{4}-\theta}$$
so you can not compute the probability of the data given $\theta$ unless you have a second parameter.
If you would have a valid MLE then you could start with the information matrix of some parameter and apply a scaling according to the square of the derivative of the transformation.
However, note that some transformed variable can be biased and the Fisher Information alone is not an indication of asymptotic variance. See this example: Why the variance of Maximum Likelihood Estimator(MLE) will be less than Cramer-Rao Lower Bound(CRLB)?
The special case of finding the variance of the distribution of the statistic $\hat\theta = \hat{p}(1-\hat{p}) = \hat{p}-\hat{p}^2$ can be done more directly by computing $$\begin{array}{} Var(\hat\theta) &=& E[(\hat{p}-\hat{p}^2)] - E[\hat{p}-\hat{p}^2]^2 \\ &=& E[\hat{p}^4-2\hat{p}^3+\hat{p}^2] - E[\hat{p}-\hat{p}^2]^2\\ &=& E[\hat{p}^4]-2E[\hat{p}^3]+E[\hat{p}^2] - (E[\hat{p}]-E[\hat{p}^2])^2 \end{array}$$
which can be expressed in terms of the raw moments of $\hat{p}$ (a binomial distributed variable, scaled by $1/n$)
$$\begin{array}{} E[\hat{p}] &=& p \\ E[\hat{p}^2] &=& p^2 + \frac{p(1-p)}{n} \\ E[\hat{p}^3] &=& p^3 + 3 \frac{p^2(1-p)}{n} \\ E[\hat{p}^4] &=& p^4 + 6 \frac{p^3(1-p)}{n} + 3 \frac{p^2(1-p)^2}{n^2} \\ \end{array}$$
and the variance can be written as
$$Var(\hat\theta) = \frac{2p^4-4p^3+2p^2}{n^2} + \frac{-4p^4+8p^3-5p^2+p}{n}$$
where the second term becomes dominant for large $n$ and is the same result as using the Fisher Information matrix
$$p(1-p)/n \cdot \left(\frac{\text{d}\theta}{\text{d}p}\right)^2 = \frac{p(1-p)(1-2p)^2}{n}$$
In the case of $p=0.5$ this would lead to zero variance (or an infinite value in the information matrix). In that case you can still use the Delta method with a second order derivative as demonstrated in this question: Implicit hypothesis testing: mean greater than variance and Delta Method