  We present a code to compute the relic density of dark matter in the
Next-to-Minimal Supersymmetric Standard Model (NMSSM). Dominant corrections to
the Higgs masses are calculated with NMHDECAY as well as theoretical and
collider constraints. All neutralino annihilation and coannihilation processes
are then computed with an extended version of micrOMEGAs, taking into acount
higher order corrections to Higgs vertices. We explore the parameter space of
the NMSSM and consider in particular the case of a bino LSP, of a mixed
bino-higgsino LSP and of a singlino LSP. As compared to the MSSM, neutralino
annihilation is often more efficient as it can take place via (additional)
Higgs resonances as well as annihilation into light Higgs states. Models with a
large singlino component can be compatible with WMAP constraints.
