  We present our computation of the ${\cal
O}((\alpha_t+\alpha_\lambda+\alpha_\kappa)^2)$ two-loop corrections to the
Higgs boson masses of the CP-violating Next-to-Minimal Supersymmetric Standard
Model (NMSSM) using the Feynman-diagrammatic approach in the gaugeless limit at
vanishing external momentum. We choose a mixed $\overline{\mbox{DR}}$-on-shell
(OS) renormalisation scheme for the Higgs sector and apply both
$\overline{\mbox{DR}}$ and OS renormalisation in the top/stop sector. For the
treatment of the infrared divergences we apply and compare three different
regularisation methods: the introduction of a regulator mass, the application
of a small momentum expansion, and the inclusion of the full momentum
dependence. Our new corrections have been implemented in the Fortran code
NMSSMCALC that computes the Higgs mass spectrum of the CP-conserving and
CP-violating NMSSM as well as the Higgs boson decays including the
state-of-the-art higher-order corrections. Our numerical analysis shows that
the newly computed corrections increase with rising $\lambda$ and $\kappa$,
remaining overall below about 3\% compared to our previously computed ${\cal
O}(\alpha_t(\alpha_t+\alpha_s))$ corrections, in the region compatible with
perturbativity below the GUT scale. The renormalisation scheme and scale
dependence is of typical two-loop order. The impact of the CP-violating phases
in the new corrections is small. We furthermore show that the Goldstone Boson
Catastrophe due to the infrared divergences can be treated in a numerically
efficient way by introducing a regulator mass that approximates the
momentum-dependent results best for squared mass values in the permille range
of the squared renormalisation scale. Our results mark another step forward in
the program of increasing the precision in the NMSSM Higgs boson observables.
