To further underscore the difference in behavior exhibited by the frequencies vs the TDMs, we used the corresponding VSE equations including quadratic electric field contributions (eq. 1a and extension of eq. 3a; see eq. S2 and eq. S1, respectively) to model the DFT-based vibrational observables solely as functions of 𝐹 ⃗ . Towards this goal, we used the polarizable AMOEBA force field to extract the electric field vectors (𝐹 ⃗ ) along the C≡N group for the DFToptimized structures (see Experimental Section for further details). All VSE parameters were allowed to freely vary when fitting the VSE equations against the DFT results (see Tables and). For the TDMs, we found that the VSE modelled the DFT results for purely electrostatic and H-bonding perturbations very well with R 2 > 0.97 (Fig. ). This is consistent with our previous experimental results that TDMs give direct access to the local nitrile electric field in both non-H-bonding and H-bonding environments (Fig. ). Further, this modelling provides a good estimate of the experimental linear field sensitivity of -1.0 𝑚𝐷 𝑀𝑉/𝑐𝑚 (as discussed in SI Section 2). In contrast to the TDMs, the C≡N vibrational frequency shifts are modelled well with eq. 1a for purely electrostatic perturbations but extremely poorly for nitriles with HBs to water molecules (Fig. ). For purely electrostatic perturbations, the correlation between the modelled and DFT frequencies is very good with R 2 of 0.95. This modelling resulted in a Stark tuning rate of |Δ𝜇 ⃗| = 0.22 ), which is impressively close to the experimental value of 0.19