But MB says the average or mean velocity is sqrt(3kT/M), where the energy of the ions (here 2eV) is not taken into account! I am confused here.
I think you are confusing bulk kinetic energy (i.e., bulk flow) and random kinetic energy (e.g., heat).
We can define moments of the distribution function as expectation values of any dynamical function, $g(\mathbf{x},\mathbf{v})$ (e.g., velocity), as:
$$
\langle g\left( \mathbf{x}, \mathbf{v} \right) \rangle = \frac{ 1 }{ N } \int d^{3}x \ d^{3}v \ g\left( \mathbf{x}, \mathbf{v} \right) \ f\left( \mathbf{x}, \mathbf{v}, t \right) \tag{1}
$$
where $\langle Q \rangle$ is the ensemble average of quantity $Q$.
Then the bulk velocity (i.e., that associated with your 2 eV energy) is given by the first velocity moment:
$$
\mathbf{U}_{s} = \frac{ 1 }{ n_{s} } \int d^{3}v \ \mathbf{v} \ f_{s}\left( \mathbf{x}, \mathbf{v}, t \right) \tag{2}
$$
where $f_{s}(\mathbf{x},\mathbf{v},t)$ is the particle distribution function of species $s$ (e.g., Maxwell-Boltzmann distribution) and $n_{s}$ is the number density (i.e., given by the zeroth velocity moment).
The random kinetic energy or thermal pressure is given by a dyadic product of velocities in the second velocity moment as:
$$
\mathbb{P}_{s} = m_{s} \int d^{3}v \ \left( \mathbf{v} - \mathbf{U}_{s} \right) \left( \mathbf{v} - \mathbf{U}_{s} \right) \ f_{s}\left( \mathbf{x}, \mathbf{v}, t \right) \tag{3}
$$
where $m_{s}$ is the particle mass of species $s$. To relate the pressure tensor to a temperature, $T_{s}$, or thermal speed, we need to assume an equation of state (e.g., ideal gas law). For an ideal gas, we find that:
$$
\langle T_{s} \rangle = \frac{ 1 }{ 3 } Tr\left[ \frac{ \mathbb{P}_{s} }{ n_{s} k_{B} } \right] \tag{4}
$$
where $Tr\left[ \right]$ is the trace operator and $k_{B}$ is the Boltzmann constant.
I have some more notes on velocity moments at https://physics.stackexchange.com/a/218643/59023.
As the ions are not "internally excited", it is in room temperature right?
No, not really. The ion temperature will depend upon how they were generated. For instance, in some cases the gas is heated using electromagnetic radiation until its thermal energy is sufficient for the atoms to begin to lose electrons. In this case, the affected and ionized atoms will have a higher temperature than the ambient or initial temperature.
In this case how the velocity distribution along each axis(x,y and z)? ... In this case how should I assume the distribution of energy along different axis?
The general form of a 3D, anisotropic, multivariate velocity distribution function, assuming uncorrelated velocities, is given by:
$$
f\left( V_{x}, V_{y}, V_{z} \right) = \frac{ A_{x} \ A_{y} \ A_{z} }{ \pi^{3/2} \ \sigma_{x} \ \sigma_{y} \ \sigma_{z} } \ e^{-\frac{1}{2} \left[ \left( \frac{ V_{x} - \mu_{x} }{ \sigma_{x} } \right)^2 + \left( \frac{ V_{y} - \mu_{y} }{ \sigma_{y} } \right)^2 + \left( \frac{ V_{z} - \mu_{z} }{ \sigma_{z} } \right)^2 \right]} \tag{5}
$$
where $V_{j}$ is the jth component of the particle velocity, $A_{j}$ is the jth component of the particle distribution function amplitude (technically, these three would be combined into a single amplitude), $\sigma_{j}^{2}$ is the variance (i.e., related to 2nd velocity moment) of the jth component of the distribution, and $\mu_{j}$ is the mean (i.e., related to 1st velocity moment) of the jth component of the distribution.
