Folding free energy is the fundamental thermodynamic quantity characterizing the stability of a protein. Yet, its accurate determination based on computational techniques remains a challenge in physical chemistry. A straightforward brute-force approach would be to conduct molecular dynamics simulations and to estimate the folding free energy from the equilibrium population ratio of the unfolded and folded states. However, this approach is not sensible at physiological conditions where the equilibrium population ratio is vanishingly small: it is extremely difficult to reliably obtain such a small equilibrium population ratio due to the low rate of folding/unfolding transitions. It is therefore desirable to have a computational method that solely relies on simulations independently carried out for the folded and unfolded states. Here, we present such an approach that focuses on the probability distributions of the effective energy (solvent-averaged protein potential energy) in the folded and unfolded states. We construct these probability distributions for the protein villin headpiece subdomain by performing extensive molecular dynamics simulations and carrying out solvation free energy calculations. We find that the probability distributions of the effective energy are well-described by the Gaussian distributions for both the folded and unfolded states due to the central limit theorem, which enables us to calculate the protein folding free energy in terms of the mean and the width of the distributions. The computed protein folding free energy (-2.5 kcal/mol) is in accord with the experimental result (ranging from -2.3 to -3.2 kcal/mol depending on the experimental methods).