  Using the specific model of a bilayer of classical charged particles (bilayer
Wigner crystal), we compare the predictions for energies and pair distribution
functions obtained by Monte Carlo simulations using three different methods
available to treat the long range Coulomb interactions in systems periodic in
two directions but bound in the third one. The three methods compared are: the
Ewald method for quasi-two dimensional systems [D.E. Parry, Surf. Sci.
$\bm{49}$, 433 (1975); \it{ibid.}, $\bm{54}$, 195 (1976)], the Hautman-Klein
method [J. Hautman and M.L. Klein, Mol. Phys. $\bm{75}$, 379 (1992)] and the
Lekner summations method [J. Lekner, Physica A$\bm{176}$, 485 (1991)]. All of
the three method studied in this paper may be applied to any quasi-two
dimensional systems, including those having not the specific symmetry of slab
systems. For the particular system used in this work, the Ewald method for
quasi-two dimensional systems is exact and may be implemented with efficiency;
results obtained with the other two methods are systematically compared to
results found with the Ewald method. General recommendations to implement with
accuracy, but not always with efficiency, the Lekner summations technique in
Monte Carlo algorithms are given.
