To determine the reaction free energies and the corresponding minimum free energy pathway of the mechanism, we performed finite temperature string calculations. We selected 19 coordinates that define the space for the reaction mechanism (see SI), which include all bond breaking and forming distances for the proton transfers and phosphate cleavage. We performed 60 iterations during the string free energy calculations, running 100 steps of QM/MM dynamics at each iteration using Langevin dynamics with a timestep of 1 fs. Each coordinate was restrained using a force constant of 150 kcal mol -1 Å 2 . The pathway was defined via 26 windows obtained from the QM/MM minimizations. After each iteration, a new pathway was calculated by fitting a polynomial of order 8 at each coordinate, and dividing the resulting 1 dimensional path within the 19 dimensional coordinate space into 26 equal arclengths. The initial coordinates were swapped using a Metropolis Monte Carlo Hamiltonian replica exchange algorithm after each iteration. MD simulations of the product complex: MD simulations of the product state were setup with the CHARMM-GUI and used the CHARMM 36 force field due to its better treatment of nucleic acids . The product state structure from the QM/MM calculations of HIV-1 RT was solvated in a cubic water box with 138 Å in each dimension and 10 Å distance to the protein. The phosphate groups were modelled as completely deprotonated. Periodic boundary conditions were applied. The system was neutralized with 150 mM potassium chloride which has been shown to be optimal for catalysis in Hs RNase H and in a subsequent step either 0 mM, 6 mM or 50 mM magnesium chloride was added. The concentration of 6 to 10 mM MgCl2 has been shown to be optimal for catalysis in vitro and the concentration of 50 mM is considered inhibitory. After minimizing for 10 000 steps, the system was equilibrated for 5 ns with restrained protein backbone atoms to allow for relaxation of side chains and to heat up the system in NVT conditions. The restraints were then lifted and the trajectories for analysis were collected over 100 ns in NPT conditions. A time-step of 2 fs was used. Nonbonded interactions were modeled with scaled 1-4 setting in NAMD2 and a cutoff of 12 Å with switching distance of 10 Å. The Particle Mesh Ewald method was used to evaluate long range electrostatic interactions and the system was coupled to a Langevin thermostat at a temperature of 303.15 K and a Langevin barostat at standard pressure. All systems were run in triplicates. NAMD 2.9 was used as the MD engine.