  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 the Maxwell 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.
