Hyperfine structure constants have many applications, but are often hard to calculate accurately due to large and canceling contributions from different terms of the hyperfine interaction operator, and also from different closed and spherically symmetric core subshells that break up due to electron correlation effects. In multiconfiguration calculations, the wave functions are expanded in terms of configuration state functions (CSFs) built from sets of one-electron orbitals. The orbital sets are typically enlarged within the layer-by-layer approach. The calculations are energy-driven, and orbitals in each new layer of correlation orbitals are spatially localized in regions where the weighted total energy decreases the most, overlapping and breaking up different closed core subshells in an irregular pattern. As a result, hyperfine structure constants, computed as expectation values of the hyperfine operators, often show irregular or oscillating convergence patterns. Large orbital sets, and associated large CSF expansions, are needed to obtain converged values of the hyperfine structure constants. We analyze the situation for the states of the odd and even configurations in N I, and show that the convergence with respect to the increasing sets of orbitals is radically improved by introducing separately optimized orbital sets targeted for describing the spin- and orbital-polarization effects of the 1s and 2s core subshells that are merged with, and orthogonalized against, the ordinary energy-optimized orbitals. In the layer-by-layer approach, the spectroscopic orbitals are kept frozen from the initial calculation and are not allowed to relax in response to the introduced layers of correlation orbitals. To compensate for this lack of variational freedom, the orbitals are transformed to natural orbitals prior to the final calculation based on single and double substitutions from an increased multireference set. The use of natural orbitals has an important impact on the states of the configuration, bringing the corresponding hyperfine interaction constants in closer agreement with experiment. Relying on recent progress in methodology, the multiconfiguration calculations are based on configuration state function generators, cutting down the time for spin-angular integration by factors of up to 50, compared to ordinary calculations.