I don't understand how DESeq2 calculates the p-values given by DESeq2::results(). I thought it was

2 * pnorm(abs(log2FoldChange / lfcSE), lower.tail = FALSE)

with log2FoldChange and lfcSE taken from the corresponding columns of the output table given by DESeq2::results(). However, when calculating the p-values this way, I get non-negligible differences (i.e. not only numerical inaccuracies) to the p-values given by DESeq2.

Sorry, just realized my mistake: I was using a threshold for the log fold change.