Why null deviance is different from my manual calculations? Let's consider this very simple example with Poisson regression:
y <- c(1, 2, 3, 2, 1, 4, 5, 6, 7, 1, 2, 3)
x <- rnorm(length(y))
z <- runif(length(y))

mod <- glm(y ~ x + z, family = poisson)

As we can observe, the null deviance is equal to:
mod$null.deviance
14.10487

However, my manual calculations (based for example on RPubs) shows something different:
# Estimate null model
mod_null <- glm(y~1, family=poisson)
2 * (logLik(mod) - logLik(mod_null))

'log Lik.' 5.812325 (df=3)

Can you please explain to me what am I doing wrong?

There is already a good answer but I thought we could still clarify some points.
Also the Note at the end shows the input we will use. It is the same as in the question except we added a seed to make the data and models reproducible.
Let $L$ and $D$ be the likelihood and the deviance of our model, respectively, and use a subscript of $0$ or $s$ to describe the corresponding quantities for the null or saturated models.  Then the code in the question calculates the null deviance, i.e. the deviance of the null model, as $2 (L - L_0)$ but the correct formula is $D_0 = -2 (L_0 - L_s)$
In terms of R code we have the following where, in R, gl(n,1) means a factor with n levels assuming n is the length of y.  In the link given in the question it is written as as.factor(1:length(y)) .
n <- length(y)
identical(as.factor(1:length(y)), gl(n, 1))
## [1] TRUE

Now calculating the saturated model and using the correct formula from above we have:
mod <- glm(y ~ x + z, family = poisson)   # model - same as in question

mod_null <- glm(y~1, family = poisson)    # null model - same as in question

n <- length(y)
mod_sat <- glm(y ~ gl(n, 1), family = poisson) # saturated model

mod$null.deviance
## [1] 14.10487

D0 <- c(- 2 * (logLik(mod_null) - logLik(mod_sat))); D0  # null deviance
## [1] 14.10487

Note
set.seed(123)
y <- c(1, 2, 3, 2, 1, 4, 5, 6, 7, 1, 2, 3)
x <- rnorm(length(y))
z <- runif(length(y))

