It is shown here how imaginary-time, path-integral methods can be easily extended to the determination of thermal properties of systems in which nuclear motion occurs simultaneously on several potential energy surfaces associated with multiple low-lying electronic states. To do so requires exponentiation of the full matrix of the potential in the diabatic basis., rather than the scalar potential which occurs when motion is limited to a single potential energy surface. In the adiabatic approximation, the partition function and subsequent thermal averages are sums of those associated, separately, with each adiabatic electronic state.