Molecular dynamics simulations over 20 ns with the Gromacs all-atom force field, coupled with cluster analyses of the trajectories, have been applied to the sixteen charge states of glutathione (GSH) in order to examine the distribution conformations in aqueous solution as a function of pH. The simulations show that GSH is very flexible and does not adopt a strongly preferred conformation at any pH. Comparison with limited conformational data deduced from NMR analyses shows little agreement. Contrary to the NMR results which found essentially equal populations of the three rotameric forms, the simulations reveal lower populations for gauche rotameric forms of the Cys and Glu side chains. In most species, the lowest populations were found for the sterically hindered gauche-gauche orientation about the C-alpha-C-beta bond in both residues, except when electrostatic attraction between oppositely charged ends, unshielded by intervening water, was found to dominate. In the majority of the enzyme-bound GSH structures extracted from the Brookhaven Protein Data Bank, the bound GSH has a conformation that is either the same or similar to that free in solution. In some cases, as in the case of solid-state GSH and the oxidized form, GSSG, crystal packing and intermolecular H-bonding interactions force the GSH skeleton into a conformation that is not seen in solution. The distribution of the separation of the Cys S atom from the Glu C-alpha-H bond was monitored over the course of the 20 ns simulations to deduce conditions under which H atom transfer may occur from the Glu C-alpha-H bond to a thiyl radical of the Cys moiety of GSH, as has been observed experimentally.