Usual Woods-Saxon single-particle levels with BCS pairing are not able to reproduce the experimental occupation probabilities of the proton and neutron levels 1p 3/2 , 1p 1/2 , 0f 5/2 , and 0g 9/2 in the double-β decay system 76 Ge to 76 Se. Shifting down the 0g 9/2 level by hand can explain the data, but it is not satisfactory. Here it is shown that a self-consistent Hartree-Fock + BCS approach with experimental deformations for 76 Among double-β decay-with-two-neutrino processes, the best-studied case is that of 76 Ge going to 76 Se. Within a spherically symmetric description of these nuclei, the relevant proton and neutron single-particle levels involved in the double-β decay process are in the valence shells 1p 3/2 , 1p 1/2 , 0f 5/2 , and 0g 9/2 . Their occupation probabilities were measured recently by Schiffer et al.[1] (for neutrons) and by Kay et al. [2] (for protons). These data are not reproduced by previous mean-field calculations using Woods-Saxon potential and pairing within BCS approximation. By increasing the binding energy of the 0g 9/2 level, the agreement is better, because its occupation probability increases and eventually reaches the experimental one. This energy shifting has been performed by hand in previous works [3][4][5] to obtain the experimental occupation probabilities, and the resulting ground-state structure has been used to compute neutrinoless double-β matrix elements. Shell-model calculations [6] have also been carried out with occupation probabilities in good agreement with experimental ones, but the single-particle energies used are not available in the literature.In this work we use a deformed Skyrme Hartree-Fock (HF) mean-field with pairing correlations within BCS approximation. At each HF iteration the BCS equations are solved to get new single-particle wave functions i , energies e i , and occupation probabilities v 2 i until convergence is achieved. After convergence, the proton-neutron quasiparticle randomphase approximation (pnQRPA) equations are solved on each deformed ground-state basis for 76 Ge and 76 Se to get their Gamow-Teller (GT) strength distributions and to compute the two-neutrino double-β decay matrix element. The formalism is similar to that used in previous works [7][8][9] with an important difference. In this work we allow the Skyrme force parameters to vary and search for a parametrization that provides better agreement with the experimental occupations of the previously mentioned valence shells in 76 Ge and 76 Se, as well as in their neighbors. We focus on the effect of the spin-orbit strength, W 0 , which we find to be of primary importance for the desired description of proton and neutron valence-shell occupations.Experimental data and HF + BCS calculations suggest prolate deformations for the 76 Ge and 76 Se ground states. The most precise data on the quadrupole moments of these nuclei come from Coulomb excitation reorientation measurements [10], which, together with experimental charge radii [11], give a deformation parameter β = 0.10 for 76 Ge ...