Monte Carlo simulations performed on multiple polymer chains have produced accurate relaxation modulus G_s(t) curves which match the experimental G(t) curves of polystyrene reasonably well, over a wide temperature range around the glass transition region. The inter-segmental interactions, defined in terms of ε* (well depth) and σ (monomer size), exert a strong influence on the modulus, the length scale and the relaxation time scale of the system. Judicious selection of these interaction parameters has enabled us to create the whole range of temperature dependence of the thermorheological complexity, from ΔT = 40 °C to ΔT = 0 °C. Near the glass transition temperature, the development of nonergodicity vis-à-vis a crowding effect in the system emerges naturally from the analysis of the G(t) line shapes. The entropic slow mode is well described by the Rouse theory and the energetic fast mode shifts to longer time scales, revealing the generic behavior of the thermorheological complexity. Typical G_s(t) curves, when partitioned into glassy and rubbery components, are shown to obey Inoue–Okamoto–Osaki's modified stress-optical rule, with different stress-optical coefficients for each component. Closer to the glass transition temperature, the distance of the closest monomer shows a considerable increase, suggesting a penetrable resistance to the approach of another monomer. The parameter σ represents the characteristic length scale of the system in the glassy region. The thermorheological complexity incorporates the dynamic length scale of structural relaxation, increasing with the decrease of temperature towards the glass transition point.