  We present in detail and validate an effective Monte Carlo approach for the
calculation of the nuclear vibrational densities via integration of molecular
eigenfunctions that we have preliminary employed to calculate the densities of
the ground and the excited OH stretch vibrational states in protonated glycine
molecule [C. Aieta et. al. Nat. Commun. 11, 4348 (2020)]. Here, we first
validate and discuss in detail the features of the method on a benchmark water
molecule. Then, we apply it to calculate on-the-fly the ab initio anharmonic
nuclear densities in correspondence of the fundamental transitions of NH and CH
stretches in protonated glycine. We show how we can gain both qualitative and
quantitative physical insight by inspection of different one-nucleus densities
and assign a character to spectroscopic absorption peaks using the expansion of
vibrational states in terms of harmonic basis functions. The visualization of
the nuclear vibrations in a purely quantum picture allows us to observe and
quantify the effects of anharmonicity on the molecular structure, and to
exploit the effect of IR excitations on specific bonds or functional groups,
beyond the harmonic approximation. We also calculate the quantum probability
distribution of bond-lengths, angles and dihedrals of the molecule. Notably, we
observe how in the case of one type of fundamental NH stretching the typical
harmonic nodal pattern is absent in the anharmonic distribution.
