This chapter addresses the methodological and computational aspects related to the prediction of molecular second-order nonlinear optical properties, i.e., the first hyperpolarizability (β), by using quantum chemistry methods. Both small (reference) molecules and extended push-pull π-conjugated systems are considered, highlighting contrasted effects about (i) the choice of a reliable basis set together with the convergence of β values as a function of the basis set size, (ii) the amplitude of electron correlation contributions and its estimate using wave function and density functional theory methods, (iii) the description of solvent effects using implicit and explicit solvation models, (iv) frequency dispersion effects in off-resonance conditions, and (v) numerical accuracy issues. When possible, comparisons with experiment are made. All in all, these results demonstrate that the calculations of β remain a challenge and that many issues need to be carefully addressed, pointing out difficulties toward elaborating black-box and computationally cheap protocols. Still, several strategies can be designed in order to achieve a targeted accuracy, either for reference molecules displaying small β responses or for molecules presenting large β values and a potential in optoelectronics and photonics.