Monte Carlo methods were used to calculate heat capacities as functions of temperature for classical atomic clusters of aggregate sizes 25 less than or equal toN less than or equal to 60 that were bound by pairwise Lennard-Jones potentials. The parallel tempering method was used to overcome convergence difficulties due to quasiergodicity in the solid-liquid phase-change regions. All of the clusters studied had pronounced peaks in their heat capacity curves, most of which corresponded to their solid-liquid phase-change regions. The heat capacity peak height and location exhibited two general trends as functions of cluster size: for N=25 to 36, the peak temperature slowly increased, while the peak height slowly decreased, disappearing by N=37; for N=30, a very small secondary peak at very low temperature emerged and quickly increased in size and temperature as N increased, becoming the dominant peak by N=36. Superimposed on these general trends were smaller fluctuations in the peak heights that corresponded to "magic number" behavior, with local maxima found at N=36, 39, 43, 46, and 49, and the largest peak found at N=55. These magic numbers were a subset of the magic numbers found for other cluster properties, and can be largely understood in terms of the clusters' underlying geometries. Further insights into the melting behavior of these clusters were obtained from quench studies and by examining rms bond length fluctuations.