Suppose that $Z \sim \mathcal N(0,1)$ and $W = \max(0,Z)$. Note that $P(W = 0) = \frac{1}{2}$ so $W$ is like a continuous random variable plus a point mass at 0. This means that the CDF $F_W$ of $W$ is
$$
F_W(w) =
\begin{cases}
0 \hspace{15mm} w < 0 \\
\frac{1}{2} \hspace{14mm} w = 0 \\
\Phi(w) \hspace{8mm} w > 0 \\
\end{cases}
$$
with a discontinuity at 0.
We want the pdf of this. Let $\delta_0$ be the probability measure defined by $\delta_0(A) = I(0 \in A)$ where $I$ is the indicator function. Let the Lebesgue measure be $m$.
Consider some measurable set $A \subset \mathbb R$. Note that $(m + \delta_0)(A) = 0 \implies P(W \in A) = 0$ therefore $P_W$ is absolutely continuous with respect to $m + \delta_0$. This means we can obtain the Radon-Nikodym pdf $f_W$ of $W$, denoted by
$$
f_W = \frac{\text d P_W}{\text d(m + \delta_0)}.
$$
Once we have this we can compute the expectation as follows:
$$
E(W) = \int_\mathbb R w \,\text dP_W = \int_\mathbb R w f_W(w) \,\text d(m + \delta_0)
$$
$$
=\int_{\{0\}} w f_W(w) \,\text d\delta_0 + \int_{(0, \infty)} w f_W(w) \,\text dm
$$
$$
= 0 + \int_{(0, \infty)} w f_W(w) \,\text dm
$$
Because $F_W(w) = \Phi(w)$ on $(0, \infty)$ we know that $f_W(w) = \phi(w)$ on this interval therefore
$$
E(W) = \int \limits_0^\infty w \phi(w) \,\text dw = \frac{1}{\sqrt{2\pi}}.
$$
We can confirm this via simulation:
z <- rnorm(1e5)
w <- sapply(z, function(v) max(c(0,v)))
mean(w)
1 / sqrt(2 * pi)
The second moment $E(X^2)$ can be found in a similar fashion. I was not rigorous with much of this, and your particular case is a little more complicated, but I think this shows the overall way to get an analytic solution out of this.
Update: $X \sim \mathcal N(\mu, \sigma^2)$
Let $X \sim \mathcal N(\mu, \sigma^2)$. Then
$$
F_W(w) =
\begin{cases}
0 \hspace{15mm} w < 0 \\
F_X(0) \hspace{6mm} w = 0 \\
F_X(w) \hspace{6mm} w > 0 \\
\end{cases}
$$
so we find that
$$
E(W) = \frac{1}{\sqrt{2\pi\sigma^2}}\int \limits_0^\infty w \exp \left( -\frac{(w - \mu)^2}{2 \sigma^2}\right) \,\text dw.
$$
We can write this in terms of the standard normal CDF $\Phi$: letting $u = \frac{w-\mu}{\sigma}$, we obtain
$$
E(W) = \frac{1}{\sqrt{2\pi\sigma^2}}\int \limits_0^\infty w \exp \left( -\frac{(w - \mu)^2}{2 \sigma^2}\right) \,\text dw =
\frac{1}{\sqrt{2\pi}}\int \limits_{-\mu/\sigma}^\infty (\sigma u + \mu) \exp \left( -u^2/2\right) \,\text du
$$
$$
= \frac1{\sqrt{2\pi}}\left[ \sigma \int \limits_{-\mu/\sigma}^\infty u e^{-u^2/2} \,\text du + \mu \int \limits_{-\mu/\sigma}^\infty e^{-u^2/2} \,\text du \right].
$$
Substituting $z = u^2/2$ in the left integral we find
$$
E(W) = \frac\sigma{\sqrt{2\pi}} e^{-\frac{\mu^2}{2\sigma^2}} + \frac\mu{\sqrt{2\pi}} \left[ \int \limits_{-\infty}^\infty e^{-u^2/2} \,\text du - \int \limits_{-\infty}^{-\mu/\sigma} e^{-u^2/2} \,\text du\right]
$$
$$
= \frac\sigma{\sqrt{2\pi}} e^{-\frac{\mu^2}{2\sigma^2}} + \mu \left( 1- \Phi(-\mu/\sigma) \right).
$$
This expression still has an integral in it but at least it is a standard integral that we are comfortable with.
In R:
x <- rnorm(1e6, mean=2, sd=1.5)
w <- sapply(x, function(v) max(c(0,v)))
## checking CDF
mean(w <= 1)
pnorm(1, 2, 1.5)
## checking expectation
mean(w)
expectation <- function(y, mu, sigma) {
y * exp(-.5 * (y - mu)^2 / sigma^2) / sqrt(2 * pi) / sigma
}
integrate(expectation, 0, Inf, mu=2, sigma=1.5)
## using simpler expression
E_W <- function(mu, sigma)
sigma / sqrt(2 * pi) * exp(-.5 * (mu/sigma)^2) + mu * (1 - pnorm(-mu/sigma))
E_W(mu=2, sigma=1.5)
Update 2
Here's the complete multivariate case. Let $X \sim \mathcal N (\mu, \Sigma)$ and $Y = (Y_1, \dots, Y_n)$ for $Y_i = \max(0, X_i)$. What I have done above is sufficient for $E(Y) = (E(Y_1), \dots, E(Y_n))$ so what remains is the variance $\Sigma_Y := Var(Y)$. The diagonal of $\Sigma_Y$ we can do already: it is just $E(Y_i^2) - E(Y_i)^2$ and we have (albeit ugly) expressions for those.
So in order to fill in the rest of $\Sigma_Y$ we need $Cov(Y_i, Y_j) = E(Y_iY_j) - E(Y_i)E(Y_j)$, so the only missing part at this point is $E(Y_iY_j)$.
The CDF of $(Y_i, Y_j)$ has a form analogous to the univariate CDFs so it turns out that
$$
E(Y_i Y_j) = \int \limits_0^\infty \int \limits_0^\infty y_1 y_2 f_{X_i, X_j}(y_1, y_2)\,\text dy_1 \,\text dy_2.
$$
Confirming this:
library(cubature)
library(mvtnorm)
## generating data
sx = 1; sy = 2; rho = .5
Sigma <- matrix(c(sx^2, sx * sy * rho, sx * sy * rho, sy^2), 2)
mu = c(-1, 2)
x <- MASS::mvrnorm(5e5, mu, Sigma)
y <- t(apply(x, 1, function(v) sapply(v, function(v2) max(c(0,v2)))))
## checking CDF
mean(y[,1] <= .5 & y[,2] <= 1.25)
mean(x[,1] <= .5 & x[,2] <= 1.25)
mvtnorm::pmvnorm(upper=c(.5,1.25), mean=mu, sigma=Sigma)
## checking E(Y_i Y_j)
mean(y[,1] * y[,2])
expectation <- function(vec, mu, sigma){
prod(vec) * mvtnorm::dmvnorm(vec, mu, sigma)
}
adaptIntegrate(expectation, lowerLimit = c(0,0), upperLimit = c(100, 100), mu=mu, sigma=Sigma)$integral
