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?

Expanding from my comment above, you are not comparing the right quantities.
The null deviance is given by 2(LL(Saturated Model) - LL(Null Model)).
Let's calculate the null deviance manually and compare with mod$null.deviance.
get_null_deviance <- function(fit) {
    
    # Get data from fit object and create new `.idx` column for
    # saturated model
    data <- transform(fit$model, .idx = factor(1:nrow(fit$model)))
    
    # Get response variable from fit object
    # In this case, this returns `"y"`
    resp <- all.vars(attr(fit$model, "terms"))[1]
    
    # Formula for saturated model and fit model
    # In this case, the saturated model is `y ~ .idx`
    fm_saturated <- reformulate(".idx", response = resp)
    fit_saturated <- glm(fm_saturated, data = data, family = poisson)
    
    # Formula for null model and fit model
    # In this case, the null model is `y ~ 1`
    fm_null <- reformulate("1", response = resp)
    fit_null <- glm(fm_null, data = data, family = poisson)
    
    # Return null deviance
    2 * (logLik(fit_saturated) - logLik(fit_null))
    
}

get_null_deviance(mod)
#'log Lik.' 14.10487 (df=12)
mod$null.deviance
#[1] 14.10487

