We propose and implement a novel computational method for simulating open-system electronic dynamics and obtaining thermalized electronic structures within an open quantum system framework. The system-bath interaction equation of motion is derived and modeled from the local harmonic oscillator description for electronic density change. The nonequilibrium electronic dynamics in a thermal bath is simulated using first-order kinetics. The resultant electronic densities are temperature-dependent and can take characteristics of the ground and excited states. We present results of calculations performed on H(2) and 1,3-butadiene performed at the Hartree-Fock level of theory using a minimal Slater-type orbital basis set.