The pair distribution function g(r) of ions in body-centered-cubic and face-centered-cubic Coulomb crystals is calculated using the harmonic-lattice (HL) approximation in a wide temperature range, from the high-temperature classical regime (T ω p ,ω p is the ion plasma frequency) to the low-temperature quantum regime (T ω p ). The radial pair distribution function g(r) is calculated by averaging g(r) over orientations of r. In the classical limit, g(r) is also obtained from extensive Monte Carlo (MC) simulations. MC and HL results are shown to be in good agreement. With decreasing temperature T, the correlation peaks of g(r) and g(r) become narrower and finally freeze at T ω p being solely determined by zero-point ion vibrations.