Quantum dots with three-dimensional isotropic harmonic confining potentials and up to 60 electrons are studied. The Dirac-Coulomb Hamiltonian serves as a framework, so that relativistic effects are included, and electron correlation is treated at a high level by the Fock-space coupled cluster method, with single and double excitations summed to all orders. Large basis sets composed of spherical Gaussian functions are used. Energies of ground and excited states are calculated. The orbital order is 1s, 2p, 3d, 3s, 4f, 4p, 5g, ... , and closed-shell structures appear for 2, 8, 18, 20, 34, 40, and 58 electrons. Relativistic effects are negligible for low strengths of the harmonic potential and increase rapidly for stronger potentials. Breit contributions, coming from the lowest order relativistic correction to the interelectronic repulsion terms, are also studied. Correlation effects are significant for these systems, in particular for weak confining potentials and for small systems, where they constitute up to 6% of the total energies. Their relative weight goes down (although they increase in absolute value) for larger systems or confining potentials. Planned applications to quantum dots with impurities are discussed briefly.