The thermodynamics of hydrogen bonds in aqueous and acidic solutions significantly impacts the kinetics and thermodynamics of acid reaction chemistry. We utilize in this work a multiscale approach, combining density functional theory (DFT) with classical molecular dynamics (MD) to model hydrogen bond thermodynamics in an acidic solution. Using thermodynamic cycles, we split the solution phase free energy into its gas phase counterpart plus solvation free energies. We validate this DFT/MD approach by calculating the aqueous phase hydrogen bond free energy between two water molecules (H2O-···-H2O), the free energy to transform an H3O+ cation into an H5O2+ cation, and the hydrogen bond free energy of protonated water clusters (H3O+-···-H2O and H5O2+-···-H2O). The computed equilibrium hydrogen bond free energy of H2O-···-H2O is remarkably accurate, especially considering the large individual contributions to the thermodynamic cycle. Turning to cations, we find the ion to be more stable than H3O+ by roughly 1-2 kBT. This small free energy difference allows for thermal fluctuation between the two idealized motifs, consistent with spectroscopic and simulation studies. Lastly, hydrogen bonding free energies between either H+ cation and H2O in solution were found to be stronger than between two H2O, though much less so than in vacuum because of dielectric screening in solution. Altogether, our results suggest the DFT/MD approach is promising for application in modeling hydrogen bonding and proton transfer thermodynamics in condensed phases.