Quantum chemistry applications on quantum computers currently rely heavily on the variational quantum eigensolver algorithm. This hybrid quantum-classical algorithm approximates ground state solutions to the Schroedinger equation using two main parts: (i) a parameterized trial wavefunction, or ansatz, is prepared on a quantum computer and the expectation value of a molecular Hamiltonian is sampled; (ii) a classical computer adjusts the parameters of the ansatz with the goal of minimizing this expectation value. This optimization procedure continues until parameters are found that minimize the expectation value, which then corresponds to the ground state electronic energy of the system at the selected geometry, following the variational principle. This process can be repeated for perturbations to each molecular degree of freedom to generate a Born-Oppenheimer potential energy surface (BOPES) for the molecule. The BOPES can then be used to derive thermodynamic properties, which are often desirable for applications in chemical engineering and materials design. To do this the nuclear Schroedinger equation is solved for some parameterization of the BOPES, where the resulting eigenvalues correspond to molecular vibrational modes. Using the partition function, thermodynamic properties can then be computed. These steps are performed on a classical computer. It is clear from this process that quantum chemistry applications contain a substantial classical computing component in addition to steps that can be performed using a quantum computer. In order to design efficient work-flows that take full advantage of each hardware-type, it is critical to consider the entire process so that the high-accuracy electronic energies possible from quantum computing are not squandered in the process of calculating thermodynamic properties. We will present a summary of the hybrid quantum-classical work-flow to compute thermodynamic properties. This work-flow contains many options that can significantly affect the efficiency and the accuracy of the results, including classical optimizer attributes, number of ansatz repetitions, how the BOPES is parameterized, and how the nuclear Schroedinger equation is solved. We will discuss the effects of these options employing robust statistics along with simulations and experiments on actual quantum hardware. We show that through careful selection of work-flow options, order-of-magnitude speed-ups are possible at equivalent accuracy.