An exhaustive search approach was used to establish all possible rotamers of α- and β-d-galactopyranose using DFT at the B3LYP/6-311+G** and M06-2X/6-311+G** levels, both in vacuum calculations, and including two variants of continuum solvent models as PCM and SMD to simulate water solutions. Free energies were also calculated. MM3 was used as the starting point for calculations, using a dielectric constant of 1.5 for vacuum modeling, and 80 for water solution modeling. For the vacuum calculations, out of the theoretically possible 729 rotamers, only about a hundred rendered stable minima, highly stabilized by hydrogen bonding and scattered in a ca. 14 kcal/mol span. The rotamer with a clockwise arrangement of hydrogen bonds was the most stable for the α-anomer, whereas that with a counterclockwise arrangement was the most stable for the β-anomer. Free energy calculations, and especially solvent modeling, tend to flatten the potential energy surface. With PCM, the total range of energies was reduced to 9–10 kcal/mol (α-anomer) or 7–8 kcal/mol (β-anomer). These figures fall to 4.5–6 kcal/mol using SMD. At the same time, the total number of possible rotamers increases dramatically to about 300 with PCM, and to 400 with SMD. Both models show a divergent behavior: PCM tends to underestimate the effect of solvent, thus rendering as the most stable many common rotamers with vacuum calculations, and giving underestimations of populations of β-anomers and gt rotamers in the equilibrium. On the other hand, SMD gives a better estimation of the solvent effect, yielding correct populations of gt rotamers, but more β-anomers than expected by the experimental values. The best agreement is observed when the functional M06-2X is combined with SMD. Both DFT models show minimal geometrical differences between the optimized conformers.