Mammalian metallothioneins (MTs) comprise a Zn3Cys9 cluster in the β domain and a Zn4Cys11 cluster in the α domain. They play a crucial role in storing and donating Zn(2+) ions to target metalloproteins and have been implicated in several diseases, thus understanding how MTs release Zn(2+) is of widespread interest. In this work, we present a strategy to compute the free energy for releasing Zn(2+) from MTs using a combination of classical molecular dynamics (MD) simulations, quantum-mechanics/molecular-mechanics (QM/MM) minimizations, and continuum dielectric calculations. The methodology is shown to reproduce the experimental observations that (1) the Zn-binding sites do not have equal Zn(2+) affinity and (2) the isolated β domain is thermodynamically less stable and releases Zn(2+) faster with oxidizing agents than the isolated α domain. It was used to compute the free energies for Zn(2+) release from the metal cluster in the absence and presence of the protein matrix (protein architecture and coupled protein-water interactions) to yield the respective disulfide-bonded product. The results show the importance of the protein matrix as well as protein dynamics and coupled conformational changes in accounting for the differential Zn(2+)-releasing propensity of the two domains with oxidizing agents.