We present a QM/MM ab initio molecular dynamics study of the peptide hydrolysis reaction catalyzed by HIV-1 protease. The QM/MM calculations are based on previous extensive classical MD simulations on the protein in complex with a model substrate (Piana, S.; Carloni, P.; Rothlisberger, U. Protein Sci. 2002, 11, 2393-2402). Gradient-corrected BLYP density functional theory (DFT) describes the reactive part of the active site, and the AMBER force field describes the rest of the protein, the solvent, and the counterions. An unbiased enhanced sampling of the QM/MM free-energy surface is performed to identify a plausible reaction coordinate for the second step of the reaction. The enzymatic reaction is characterized by two reaction free-energy barriers of similar to18 and similar to21 kcal mol(-1) separated by a metastable gem-diol intermediate. In both steps, a proton transfer that involves the substrate and the two catalytic Asp molecules is observed. The orientation and the flexibility of the reactants, governed by the surrounding protein frame, are the key factors in determining the activation barrier. The calculated value for the barrier of the second step is slightly larger than the value expected from experimental data (similar to16 kcal mol(-1)). An extensive comparison with calculations on gas-phase model systems at the Hartree-Fock, DFT-BP, DFT-BLYP, DFT-B3LYP, MP2, CCD, and QM/MM DFTBLYP levels of theory suggests that the DFT-BLYP functional has the tendency to underestimate the energy of the gem-diol intermediate by similar to5-7 kcal mol(-1).