Glutaredoxins are small enzymes that catalyze the oxidation and reduction of protein disulfide bonds by the thiol-disulfide exchange mechanism. They have either one or two cysteines in their active site, resulting in different catalytic reaction cycles that have been investigated in many experimental studies. However, the exact mechanisms are not yet fully known, and to our knowledge, no theoretical studies have been performed to elucidate the underlying mechanism. In this study, we investigated a proposed mechanism for the reduction of the disulfide bond in the protein HMA4n by a mutated monothiol Homo sapiens glutaredoxin and the co-substrate glutathione. The catalytic cycle involves three successive thiol-disulfide exchanges that occur between the molecules. To estimate the regioselectivity of the different attacks, classical molecular dynamics simulations were performed and the trajectories analyzed regarding the sulfur-sulfur distances and the attack angles between the sulfurs. The free energy profile of each reaction was obtained with hybrid quantum mechanical/molecular mechanical metadynamics simulations. Since this required extensive phase space sampling, the semi-empirical density functional tight-binding method was used to describe the reactive cysteines. For an accurate description, we used specific reaction parameters fitted to B3LYP energies of the thiol-disulfide exchange and a machine learned energy correction that was trained on coupled-cluster single double perturbative triple [CCSD(T)] energies of thiol-disulfide exchanges. Our calculations show the same regiospecificity as observed in the experiment, and the obtained barrier heights are about 12 and 20 kcal/mol for the different reaction steps, which confirms the proposed pathway.