How can a cosmological constant be a candidate for dark energy if the universe is flat?

In a word, No: a cosmological constant does not guarantee your space is not flat.

Take the metric $$ ds^2 = -dt^2 + e^{2Ht}(dx^2 + dy^2 + dz^2). $$

This is the flat slicing of De Sitter space (you can find it in "S.W. Hawking and G.F.R. Ellis, The large scale structure of space-time (Cambridge University Press, Cambridge, England, 1973)" pag 125 and following): De Sitter space is partitioned into two flat regions along a diagonal cut.

You can read a good discussion of this in "Steady-State Eternal Inflation" by Anthony Aguirre and Steven Gratton (https://arxiv.org/pdf/astro-ph/0111191.pdf). You also find a good picture there: FIG. 1 on page 2.

The point is that you can always see a CC as a new type of "exotic field", which, then, is brought on the right side of EFE and, if it cancels out neatly with the energy-stress tensor of matter (in our universe, it looks like it does), makes null the curvature.

This particular spacetime (the flat slicing of De Sitter space) was used by Hoyle and Narlikar as the setting of the steady state model: it is not geodesically complete, but Aguirre and Gratton argue that the two flat regions cannot communicate without causality violation.

In Friedmann cosmology, the sign of spatial curvature is an independent model parameter, whereas the sign of spacetime curvature is determined by the energy contents of the universe (including dark energy as a proxy for a cosmological constant) as well as the equation of state.

Note that by spatial curvature, we mean the curvature of the hyperplanes of constant cosmological time where the energy distribution is homogeneous (remember, there generally is no canonical space/time decomposition).

From tracing over the Einstein equation, spacetime curvature is given by $$ \mathrm{R} = 8\pi(1-3w)\rho $$ where we assumed a perfect fluid with equation of state $\rho = wP$.

Getting an expression for the spatial curvature is a bit more involved as you'll actually have to compute the Ricci tensor, but you should arrive at $$ {}^{(3)}\mathrm{R} = \frac{6k}{R_0^2a^2} $$ where $k$ is a model parameter giving its sign.

In the chosen conventions with dimensionless $k$ and $a$, the Friedmann equations should read $$ \dot a^2 = \frac {8\pi}3 \rho a^2 - \frac k{R_0^2} \\ \ddot a = -\frac {4\pi}3 (1+3w) \rho a $$ where energy density and scale factor are related by $$ \rho = \rho_0 a^{-3(1+w)} $$

A matter dominated universe corresponds to $w=0$, yielding $\mathrm R > 0$ and $\ddot a < 0$. For a universe dominated by radiation or ultra-relativistic particles, we have $w=1/3$ and thus $\mathrm R=0$ and still $\ddot a < 0$.

In a dark energy dominated universe, the situation is a bit different as we can have a positive as well as negative energy density and its sign (corresponding to the sign of the cosmological constant) directly determines the sign of $\mathrm R$ as well $\ddot a$. As you mentioned, this yields de Sitter and anti de Sitter space respectively. De Sitter space can be sliced according to any choice of $k$ (cf Wikipedia for the flat, hyperbolic and spherical slicings). In contrast, for a negative comological constant and non-negative $k$, the right-hand side of the first Friedmann equation will be negative and you won't find real-valued solutions.

In conclusion: Adding dark energy to a universe will a priori not affect its spatial geometry, but its expansion rate. However, a dominant negative cosmological constant is only possible for hyperbolic geometries.