An efficient parallelization approach to simulate optical properties of ensembles of quantum emitters in realistic electromagnetic environments is considered. It relies on balancing computing load of utilized processors and is built into three-dimensional domain decomposition methodology implemented for numerical integration of Maxwell's equations. The approach employed enables directly accessing dynamics of collective effects as the number of molecules in simulations can be drastically increased. Numerical experiments measuring speedup factors demonstrate the efficiency of the proposed methodology. As an example, we consider dynamics of nearly 700,000 diatomic molecules with ro-vibrational degrees of freedom explicitly accounted for coupled to electromagnetic radiation crafted by periodic arrays of split-ring resonators and triangular nanoholes. As an application of the approach, dissociation dynamics under strong coupling conditions is scrutinized. It is demonstrated that the dissociation rates are significantly affected near polaritonic frequencies.