The prediction of stabilities of ion pairs relative to the corresponding neutral pairs is a challenging computational problem but very important in biological systems where the occurrence of such pairs is often mechanistically significant. Here we have computed the relative free energy for the protonated and neutral forms of the mechanism-based substrate 8-methylpterin bound to dihydrofolate reductase (DHFR) in a ternary complex with cofactor nicotinamide adenine dinucleotide phosphate (NADPH). The free energy components required to calculate relative affinities of the protonated substrate hydrogen (H) bonded to the unprotonated form of the conserved active-site carboxylate residue (Glu30 for chicken DHFR) and neutral substrate H bonded lo the neutral (carboxylic acid) form of Glu30 were obtained using both ab initio QM methods and combined semiempirical (AM1) QM and MM (QM/MM) methods. The energy difference for the H-bonded systems was first calculated in a vacuum using ab initio QM methods for a model system consisting of the Glu30 side chain and the substrate. The free energy of interaction of this system with the protein/solvent surroundings was then computed using a coordinate-coupled free energy perturbation (FEP) method implemented within the molecular dynamics (MD) simulation scheme. The ab initio QM results show the neutral form is more stable in a vacuum, but that this relative stability is reduced in the enzyme due to the more favorable electrostatic interactions made by the ion pair. The neutral pair is calculated to be approximately 5 +/-2 kcal/mol (Delta Delta G(bind)) more stable than the ion pair, but the results are suggestive of a small free energy gap and the possibility of a low energy barrier for the proton transfer. The results are in qualitative agreement with other theoretical studies on similar types of H-bonded systems but appear to contradict our own experimental data for these DHFR ligands, Reasons for the apparent discrepancy are discussed.