We use density functional theory (B3LYP) and the COSMO continuum solvent model to characterize the structure and stability of the hydrated Cu(II) complexes [Cu(MeNH 2 2-x (x ) 1-3) as a function of metal coordination number (4-6) and cluster size (n ) 4-8, 18). The small clusters with n ) 4-8 are found to be the most stable in the nearly square-planar four-coordinate configuration,, which is three-coordinate. In the presence of the two full hydration shells (n ) 18), however, the five-coordinate square-pyramidal geometry is the most favorable for Cu (MeNH 2 ) 2+ (5, 6) and Cu(OH) + (5, 4, 6), and the four-coordinate geometry is the most stable for Cu(OH) 2 (4, 5) and Cu(OH) 3 -(4). (Other possible coordination numbers for these complexes in the aqueous phase are shown in parentheses.) A small energetic difference between these structures (0.23-2.65 kcal/mol) suggests that complexes with different coordination numbers may coexist in solution. Using two full hydration shells around the Cu 2+ ion (18 ligands) gives Gibbs free energies of aqueous reactions that are in excellent agreement with experiment. The mean unsigned error is 0.7 kcal/mol for the three consecutive hydrolysis steps of Cu 2+ and the complexation of Cu 2+ with methylamine. Conversely, calculations for the complexes with only one coordination shell (four equatorial ligands) lead to a mean unsigned error that is >6.0 kcal/mol. Thus, the explicit treatment of the first and the second shells is critical for the accurate prediction of structural and thermodynamic properties of Cu(II) species in aqueous solution.