The nucleosome core particle (NCP) is the basic building block of chromatin. Under the influence of multivalent cations, isolated mononucleosomes exhibit a rich phase behaviour forming various columnar phases with characteristic NCP-NCP stacking. NCP stacking is also a regular element of chromatin structure in vivo. Understanding the mechanism of nucleosome stacking and the conditions leading to self-assembly of NCPs is still incomplete. Due to the complexity of the system and the need to describe electrostatics properly by including the explicit mobile ions, novel modelling approaches based on coarse-grained (CG) methods at the multiscale level becomes a necessity. In this work we present a multiscale CG computer simulation approach to modelling interactions and self-assembly of solutions of NCPs induced by the presence of multivalent cations. Starting from continuum simulations including explicit three-valent cobalt(III)hexammine (CoHex(3+)) counterions and 20 NCPs, based on a previously developed advanced CG NCP model with one bead per amino acid and five beads per two DNA base pair unit (Fan et al 2013 PLoS One 8 e54228), we use the inverse Monte Carlo method to calculate effective interaction potentials for a 'super-CG' NCP model consisting of seven beads for each NCP. These interaction potentials are used in large-scale simulations of up to 5000 NCPs, modelling self-assembly induced by CoHex(3+). The systems of 'super-CG' NCPs form a single large cluster of stacked NCPs without long-range order in agreement with experimental data for NCPs precipitated by the three-valent polyamine, spermidine(3+).