  An exact approach to compute physical properties for general
multi-electronic-state (MES) systems in thermal equilibrium is presented. The
approach is extended from our recent progress on path integral molecular
dynamics (PIMD) [J. Chem. Phys. 145, 024103 (2016); 147, 034109 (2017)] for
quantum statistical mechanics when a single potential energy surface is
involved. We first define an effective potential function that is numerically
favorable for MES-PIMD, and then derive corresponding estimators in MES-PIMD
for evaluating various physical properties. Its application to several
representative one-dimensional and multi-dimensional models demonstrates that
MES-PIMD in principle offers a practical tool in either of the diabatic and
adiabatic representations for studying exact quantum statistics of
complex/large MES systems when the Born-Oppenheimer approximation, Condon
approximation, and harmonic bath approximation are broken.
