Molecular aggregation plays an important role on optical phenomena, and optimizing the latter requires understanding the aggregation behavior as a function of the molecular structures. The formation of molecular aggregates composed of stilbazolium chromophores and small to medium-size anions is studied by carrying out classical molecular dynamics (MD) simulations, shedding light on the dynamics of spontaneous self-aggregation of anion-cation (AC) pairs and estimating their equilibrium aggregation constants. First, analyses based on atom-atom distances and on contact mapping (Voronoi tessellation) highlight that the interactions occur at the level of the methylpyridinium groups of the cations, which enabled monitoring the clustering dynamics. Then, simulations show that clusters containing up to 5 AC pairs dominate the populations. Among the AC pairs, two of them (where the stilbazolium bears a strong donor substituent and when the anion is SCN- or pTS-) aggregate into small clusters even at high concentration, corroborating related experimental results [ChemPhysChem. 11 (2010) 495–507]. Still, these MD simulations do not reveal strict relationships between the aggregate size and the nature of the anion. Nevertheless, AC pairs bearing iodide anions require longer MD simulation times to reach an asymptotic equilibrium tendency of the average number of clusters. Finally, the aggregates are shown to adopt a preferential Λ-shape with contributions from head-to-head and stacked forms. Complementarily, the analyses of the isoperimetric quotient reinforce the AC pairs are not forming spherical aggregates even for clusters including more than 10 AC pairs. So, when the stilbazolium bears a strong donor substituent and when the anion is pTS-, stacked aggregates are formed, at both low and high concentrations, resulting in excellent agreement with the measured aggregation equilibrium constant.