Conceptual question on estimation : How to calculate the variance of estimation error EDIT/ UPDATE:
I have understood CRLB & why we need it. But my problem is something else. In book Introduction to statistical Signal processing with Applications, by M.D. Srinath, P.K. Rajasekaran and R. Viswanathan,
Chap 5, Sec5.5 Properties of the estimators Eq(5.5.10):
$\textrm{Var}[\hat{\theta}_{ml} −θ]=E(\hat{\theta}_{ml} −θ)^2 \ge \frac{1}{I}$,  the R.H.S in this expression represents lower bound on the variance of any unbiased estimate of $\theta$  & if satisfied with an equality then is called an efficient estimator. 
So, for an estimator to be efficient, the variance of the estimation error must be less than the Cramer - Rao Lower bound (CRLB): 
\begin{equation}
\textrm{Var}[\hat{\theta_{ml}} - \theta] \ge -\frac{1}{E\left[ \displaystyle  \frac{\partial^2 \ln p(z|\theta)}{\partial \theta^2} \right]}  \tag{1}
\end{equation}
For the example : $Z_i = \theta + V_i$, $i =1,..,N$ where
$p(z|\theta) = \displaystyle \prod_{i=1}^{N} \frac{1}{\sqrt{2\pi \sigma^2_v}} \exp\left({-\frac{(z_i - \theta)^2}{2\sigma_v^2}}\right)$
I found out that $$\hat{\theta_{ml}} = \frac{1}{N}\sum_{i=1}^{N}z_i  \tag{2}$$ and is an efficient estimator. 
The value of
\begin{equation}E\bigg[\frac{\partial^2 \ln p(z|\theta)}{\partial \theta^2} \bigg]  = -\frac{N}{\sigma_v^2} \tag{3}\end{equation} 
Question: 
My problem is that in order to check if it is a minimum variance (MVU) estimator, I need to calculate the LHS of $(1)$. In book, they directly equate $(1)$ without calculating  what $\textrm{Var}[\hat{\theta}_{ml} - \theta]$ actually is. How do I calculate this (LHS of $(1)$) quantity.
My confusion is should we not calculate in the first place what the variance of the estimator is in order to check if the expression of the variance is equal to 1/I ?I do not know the technique of calculating $\textrm{Var}[\hat{\theta}_{ml} −\theta]$  of the example I posted in question & hence unable to determine whether it is $1/I$ . How will I calculate separately the expression for $\textrm{Var}[\hat{\theta}_{ml} −\theta]$ . I could calculate the expression for $1/I$  as given in Eq(3). This will be helpful in cases where pdf and regularization is not straight forward as a plug in.
I tried calculating what the value of the L.H.S in $(1)$ will be but got stuck. In general 
$\textrm{Var}(X) = E[X^2] - (E[X])^2$ where $X = \hat{\theta}_{ml} - \theta$.
$\textrm{Var}(\hat{\theta}_{ml}) - \textrm{Var}(\theta) = \textrm{Var}\left( \displaystyle \frac{1}{N}\sum_{i=1}^{N}z_i \right) - \textrm{Var}(\theta)$. 
Then what should I do next so that I get the expression $\textrm{Var}(\hat{\theta}_{ml} - \theta) = \sigma_v^2/N$
This may sound too trivial, but in other problem exercises I am unable to calculate the variance. Can somebody please show how I can calculate the variance of the estimation error?
 A: If the regularity condition is satisfied prior to $(1)$ above, that is if 
$E\left[\displaystyle \frac{\partial\ln p(\mathbf{z}|\theta)}{\partial\theta}\right] = 0$  for all $\theta$. 
Then an unbiased estimator attaining CRLB may be found if and only if
\begin{equation}
\displaystyle \frac{\partial\ln p(\mathbf{z}|\theta)}{\partial\theta} = I(\theta)\left(g(\mathbf{z})-\theta\right) \tag{4}
\end{equation}
Then the MVU estimator is $\hat{\theta} = g(\mathbf{z})$ and the minimum variance is $\displaystyle \frac{1}{I(\theta)}$ where $I(\theta)$ is the Fisher information.
Back to your question, a way of doing it is to express the partial derivative of the $PDF$ is the form in $(4)$. Or you can do partial differentiation twice and use the negative reciprocal in $(1)$ as you have it in $(3)$. But let's stick to the first method because if the CRLB is attained then $\textrm{Var}(\hat{\theta}) = \displaystyle \frac{1}{I(\theta)}$.
Allow me the change of notation from capital to small. And I'm assuming $v$ is WGN with variance $\sigma_v^2$. 
\begin{align}
p(\mathbf{z}|\theta) = &\displaystyle \prod_{i=1}^N\frac{1}{\sqrt{2\pi\sigma_v^2}}\exp\left[-\frac{1}{2\sigma_v^2}\left(z_i - \theta\right)^2\right]\\
=&\frac{1}{\left(2\pi\sigma_v^2\right)^{\frac{N}{2}}}\exp\left[-\frac{1}{2\sigma_v^2}\sum_{i=1}^N\left(z_i - \theta\right)^2\right]\\
\Rightarrow \frac{\partial\ln p(\mathbf{z}|\theta)}{\partial\theta}=&\frac{1}{\sigma_v^2}\sum_{i=1}^N\left(z_i - \theta\right) = \frac{N}{\sigma_v^2}\left(\bar{z} - \theta\right)
\end{align}
Comparing with $(4)$, the minimum variance is ${\sigma_v^2}/N$.
EDIT:
If you still want to calculate the variance the usual way as you put it, proceed with the equation you started above. Following your notations:
\begin{align}
\textrm{Var}(\hat{\theta}_{ml}) - \textrm{Var}(\theta) = & \displaystyle \textrm{Var}\left(\frac{1}{N}\sum_{i=1}^N z_i\right) - \textrm{Var}(\theta)\\
= &\frac{1}{N^2}\sum_{i=1}^N \textrm{Var}(z_i) - \textrm{Var}(\theta)\\
= & \frac{1}{N^2}(N\sigma_v^2) - \textrm{Var}(\theta)\\
= & \frac{\sigma_v^2}{N}
\end{align}
