  Monte Carlo methods were used to calculate heat capacities as functions of
temperature for classical atomic clusters of aggregate sizes $25 \leq N \leq
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.
