The electric field-induced second harmonic generation (EFISHG) response has been largely used to describe the first β and the second γ hyperpolarizabilities in solution. Although the EFISHG technique cannot be applied to charged compounds (due to the external static electric field) it can be used to describe ion pairs as neutral complexes. A multiscale computational approach is required to generate representative geometrical configurations of such kind of complexes (using classical force fields), to compute the electronic structure of each configuration (using quantum mechanics methods), and to perform statistical analyses describing the behavior of the nonlinear optical properties. In this work, we target solvated neutral ion pairs complexes, of which the cation is an organic chromophore, and we estimate their EFISHG and hyper-Rayleigh Scattering responses. It is shown that the anion-cation relative spatial distribution determines the permanent dipole moment of the complexes and therefore the relative distance controls the EFISHG response. On the other hand, the β tensor is independent of the dipole moment and it shows a weak linear correlation with the π-electron conjugation length of the cations. The γ contributions in the global EFISHG response ranged from 5% to 15% mainly due to the differences in the μ and β vectors orientation. The applied multiscale approach provides reasonable results compared with the experiment although additional efforts are still required to improve such comparison mainly to consider the possible dissociation effects.