Heat capacity curves as a function of temperature are a simple probe for determining whether and where a folding transition takes place. The melting point of a foldamer appears as a peak in heat capacity because thermodynamic work is required to break/make the native non-covalent interactions. With a fixed size simulation box, constant volume heat capacity (C V ) is computed from the derivative of decorrelated potential energy with respect to temperature using a finite difference scheme. The Multistate Bennett Acceptance Ratio (MBAR) method is used to compute the expectation value of heat capacity at temperatures of interest, including arbitrary intermediate temperatures not sampled during the replica exchange simulation. By reweighting the reduced potential energies from a replica exchange simulation to the temperatures of interest, the energies of each uncorrelated configuration sampled at each thermodynamic state can be leveraged in the computation of C V at all states of interest, provided that there is sufficient configurational overlap between the thermodynamic states. Moreover, MBAR provides an asymptotically unbiased estimate of uncertainty at both sampled and unsampled states. In this work we have utilized the MBAR implementation in the Python package pymbar. For all expectation calculations of thermodynamic observables, including heat capacity, native contact fraction, and free energy of folding as a function of T, uncertainties can also be reliably calculated using bootstrap resampling. Uncertainty in detailed features of the heat capacity curve, such as melting point and full-width half maximum, can also be computed in this way. In this work, we implement a variant of standard bootstrapping in which the initial reference frame for slicing decorrelated time series data is also randomly selected in each bootstrap trial. This effectively allows all data points to be utilized in the calculation, leading to smoother mean values.