Highly accurate force fields are a mandatory requirement to generate predictive simulations. Here we present the path for the construction of machine learned molecular force fields by discussing the hierarchical pathway from generating the dataset of reference calculations to the construction of the machine learning model, and the validation of the physics generated by the model. We will use the symmetrized gradient-domain machine learning (sGDML) framework due to its ability to reconstruct complex high-dimensional potential energy surfaces (PES) with high precision even when using just a few hundreds of molecular conformations for training. The data efficiency of the sGDML model allows using reference atomic forces computed with high-level wave-function-based approaches, such as the gold standard coupled-cluster method with single, double, and perturbative triple excitations (CCSD(T)). We demonstrate that the flexible nature of the sGDML framework captures local and non-local electronic interactions (e.g., H-bonding, lone pairs, steric repulsion, changes in hybridization states (e.g., sp2⇌sp3), n → π∗ interactions, and proton transfer) without imposing any restriction on the nature of interatomic potentials. The analysis of sGDML models trained for different molecular structures at different levels of theory (e.g., density functional theory and CCSD(T)) provides empirical evidence that a higher level of theory generates a smoother PES. Additionally, a careful analysis of molecular dynamics simulations yields new qualitative insights into dynamics and vibrational spectroscopy of small molecules close to spectroscopic accuracy.