Context. Spectropolarimetric interpretation of chromospheric lines requires solving the radiative transfer problem under non-local thermodynamic equilibrium (NLTE) conditions. The atomic level populations must be computed self-consistently with the radiation field. Traditional inversion codes employ 1.5D approximations for computational tractability, neglecting the horizontal radiative transfer that becomes significant near magnetic structures and dynamic chromospheric features. Aims. We present a method for solving the 3D statistical equilibrium (SE) populations of any atom using graph neural networks (GNNs), extending prior work in 1.5D to the full 3D domain. The goal is to develop a fast surrogate model that fully accounts for horizontal and vertical radiative coupling while being robust to changes in the underlying grid, enabling future 3D NLTE inversions. Methods. We discretized the solar atmosphere as a directed graph, where nodes encode local physical properties (temperature, velocity, magnetic field, electron density, etc.) and edges encode geometric distances between spatial points. An encode-process-decode GNN architecture propagates information across the 3D domain to capture the radiative coupling efficiently. The network is trained on a Bifrost radiation magnetohydrodynamic (MHD) simulation using Ca II populations computed with Multi3D as ground truth. Results. The trained GNN accurately predicts populations of the five-level Ca II atom plus continuum. Correlations exceed 0.99 in the photosphere and mid-chromosphere; errors increase in the upper chromosphere, where non-local effects dominate but remain unbiased with Gaussian distributions. Inference is ∼106 faster than traditional iterative SE solvers. Spectral synthesis of the Ca II 8542 Å line yields intensity profiles with ≲2% mean residuals relative to the full 3D solution, opening up the option of more realistic inversions of strong chromospheric lines. Conclusions. The 3D GNN framework bypasses the computational bottleneck of iterative SE solvers while preserving essential NLTE physics, including horizontal radiative transfer. The method is readily extensible to other atoms and paves the way toward routine 3D NLTE inversions.