We present the first complete treatment for calculating theoretical estimates of free energies of formation of macromolecule−ligand complexes with molecular dynamics simulations, as the free energy for transforming the ligand into a non-interacting state by gradually diminishing the forces between macromolecule (plus solvent) and ligand. The calculations become possible due to the introduction of a specially designed potential (“molecular tweezers”) which restrains the spatial position and orientation of the ligand molecule and is gradually applied as the transformation proceeds from complexed to non-interacting components. The binding of benzene to a mutant T4 lysozyme (Morton et al. Biochemistry 1995, 34, 8564−8575) has been used as a test case. The simulations reproduce the value of the free energy of binding (−5.19 kcal/mol if the standard state of benzene is a 1 M aqueous solution) within the sum of experimental and statistical error. Another series of such simulations with rigid protein models provi...