We have applied molecular dynamics umbrella-sampling simulations and ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT) to explore the free energy surface, reaction paths, and dideuterium kinetic isotope effect (KIE) of the β-oxidation of butyryl-coenzyme A by short-chain acyl-CoA dehydrogenase. The potential energy surface is obtained by combined quantum mechanics-molecular mechanics (QM/MM) with specific reaction parameters and a simple valence bond term. The calculations include determination of the potential of mean force (PMF) in both two dimensions (2D) and one dimension (1D) by using the weighted histogram analysis method. The 2D PMF indicates that the hydride transfer is the rate-limiting step of the mechanism, and the 1D PMFs are used to calculate rate constants. We include full molecular dynamics of a 7824-atom reaction zone with stochastic boundary conditions in a first approximation of the quasiclassical rate constant. This first approximation is then corrected by transmission coefficients that account for dynamical recrossing and tunneling. For the calculation of transmission coefficients, the system is divided into a 43-atom primary zone, consisting of atoms directly involved in the enzyme reaction, and a 23 234-atom secondary zone, called the bath. In both approximations considered, the primary zone atoms are allowed to have kinetic energy in all steps. For calculation of the transmission coefficient in the static-secondary-zone (SSZ) approximation, the bath is frozen for each member of an ensemble of transition state configurations and associated reaction paths, which are then averaged, whereas in the equilibrium-secondary-zone (ESZ) approximation 7781 atoms of the bath are fully relaxed (with stochastic boundary conditions) along each reaction path. We found that the SSZ simulations underestimate the KIE, whereas the ESZ simulations overestimate the KIE.