A Critical Evaluation of Vibrational Stark Effect (VSE) Probes with the Local Vibrational Mode Theory.

Over the past two decades, the vibrational Stark effect has become an important tool to measure and analyze the in situ electric field strength in various chemical environments with infrared spectroscopy. The underlying assumption of this effect is that the normal stretching mode of a target bond such as CO or CN of a reporter molecule (termed vibrational Stark effect probe) is localized and free from mass-coupling from other internal coordinates, so that its frequency shift directly reflects the influence of the vicinal electric field. However, the validity of this essential assumption has never been assessed. Given the fact that normal modes are generally delocalized because of mass-coupling, this analysis was overdue. Therefore, we carried out a comprehensive evaluation of 68 vibrational Stark effect probes and candidates to quantify the degree to which their target normal vibration of probe bond stretching is decoupled from local vibrations driven by other internal coordinates. The unique tool we used is the local mode analysis originally introduced by Konkoli and Cremer, in particular the decomposition of normal modes into local mode contributions. Based on our results, we recommend 31 polyatomic molecules with localized target bonds as ideal vibrational Stark effect probe candidates.


Introduction
The Stark effect refers to the response of a spectroscopic transition to an applied electric field. It was discovered by Stark in 1913 observing that an applied electric field causes a splitting in the absorption lines of hydrogen [1]. In 1995, Chattopadhyay and Boxer reported that the infrared absorbance of the CN stretching mode in the anisonitrile molecule changed proportionally with the strength of an electric field imposed on the molecule [2]. Since then, the vibrational Stark spectroscopy (VSS) has become an important analytical method, which has been summarized over the past decade in a series of review articles [3][4][5][6][7][8][9][10][11]. The vibrational Stark effect (VSE) describes the perturbation of a vibrational frequency by an electric field [2,3,5,12] according to Equation (1): where ν and ν 0 are the vibrational frequencies of a specific molecular vibrational mode (i.e., the target bond stretching mode in most cases) with (ν) and without (ν 0 ) an external electric field F, respectively; ∆ µ is the difference dipole moment (also known as Stark tuning rate) and ∆α is the difference polarizability in a VSE experiment. The electric field strength is in general below 100 MV/cm; therefore, the quadratic term with regard to ∆α in Equation (1) can be neglected, so that the change in the vibrational frequency ∆ν = ν − ν 0 directly correlates with the change in the strength of the electric field F [5]. This linear relationship between vibrational frequency and electric field has formed the basis for the vibrational Stark spectroscopy. Given a simplified electrostatic description of non-covalent interactions between the vibrational probe and surrounding molecules, the strength of these intermolecular interactions can be assessed by the electric field a target chemical bond feels, as revealed by the VSE [5,13]. The VSE has been extensively applied to study the non-covalent interactions in different types of chemical systems and environments including proteins/enzymes [6][7][8]10,11,[14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33], nucleic acids [34,35], ionic liquids [36,37], biological membranes [38], electrochemical interfaces/surfaces [12,[39][40][41][42][43], and polymers [3,44,45]. Recently, the range of applications has been extended to the investigation of water clusters [46,47] and molecular solids [48].
The normal stretching vibration of a probe bond (e.g., the C=O bond in formaldehyde) is considered to be largely decoupled from rest of the molecule, i.e., its associated normal mode is ideally localized, which is generally not the case [52][53][54][55][56]; 2.
The vibrational frequency shift ∆ν arising from changes in the vicinal environment of the probe molecule can be fully attributed to the external electric field. This is the basic foundation for using the VSE as a tool to characterize non-covalent interactions; 3.
The difference dipole moment ∆ µ in Equation (1) is unaffected by the external electric field F, so that the vibrational frequency shift ∆ν responds to F in a linear fashion; 4.
The linear relationship between vibrational frequency and the electric field, observed for a relatively weak electric field strength (in the order of 1 MV/cm) will also hold for the binding pocket of proteins, where the effective electric field caused by the enzyme environment could be a hundred times stronger.
The first assumption is the most important as the vibrational Stark effect is based on a simplified model assuming that the probe bond stretching vibration encodes all information about the surrounding electric field. However, to the best of our knowledge, no systematic study on the extent to which those commonly applied and/or potential vibrational Stark effect probes can meet this requirement, has been reported so far. To fill this gap, we used in this work as powerful tool the characterization of normal mode (CNM) procedure, which is an important part of the local vibrational mode analysis originally developed by Konkoli and Cremer [57][58][59][60]. CNM determines quantitatively to what extent the local stretching vibrational mode of the probe bond is decoupled from the other local vibrational modes of the probe, and therefore provides a unique measure to assess the qualification of a probe molecule.
This paper is structured in the following way: First, the local vibrational mode theory including the CNM method is summarized and it is discussed how CNM can be applied to evaluate a vibrational Stark effect probe in the Methodology part. Then, the Computational Details are given. In the Results and Discussion part, 68 VSE probes and candidates with C=O, C≡N, S=O and other types of probe bonds are analyzed and scored with the CNM approach. A complete set of 107 VSE probes is given in the Supplementary Information. Furthermore, the sensitivity of the calculated scores with regard to density functional is checked. Lastly, 31 probe molecules with high scores are recommended for experimental verification and application.

Methodology
The VSE and its related spectroscopy require that the normal vibrational mode of the probe bond stretching is decoupled from the local vibrational modes led by other internal coordinates within the probe molecule [5,51], which does not comply with the fact that normal vibrational modes are generally delocalized over several part of a molecule or even the whole molecule because of the mass coupling [52][53][54][55][56]61].
A prominent example is the popular Tolman electron parameter (TEP) which assesses the metal-ligand (ML) bond strength in nickel-tricarbonyls [Ni(CO) 3 L] indirectly using the A 1 -symmetrical CO stretching mode as probe. The TEP rests upon the assumption that the A 1 -symmetrical CO stretching mode is fully localized and does not couple with other local modes [62][63][64]. However, our local mode analysis clearly revealed that this assumption is generally not true [65][66][67]. This indicates that, also in a polyatomic VSE probe, the normal vibrational mode of the probe bond stretching may not be ideally localized, which will impact its qualification for characterizing VSE. Therefore, a method is needed to quantify the local character of the probe bond vibration and it can be easily applied to existing VSE probes and potential probe candidates. In the following, we will review the CNM method, developed by Konkoli, Larsson and Cremer in 1998 [59,60,68,69] within the framework of the local vibrational mode theory [57], which is the perfect tool to determine in a quantitative way to what extent the normal vibration of the probe bond stretching is consisting of pure stretching character.
The harmonic normal vibrational modes and corresponding frequencies for a polyatomic molecular system with N atoms in its equilibrium geometry can be obtained by solving the Wilson equation of vibrational spectroscopy [52,70]: where f x is the Hessian matrix in terms of 3N Cartesian coordinates with the dimension of 3N × 3N. Matrix M is the diagonal mass matrix accounting for all N atoms in x, y, and z directions. The diagonal matrix Λ in the (N vib ×N vib ) dimension contains N vib vibrational eigenvalues λ µ (µ = 1, ..., N vib with N vib = 3N − K) and K equals six or five for nonlinear or linear molecules, respectively. The (3N×N vib ) dimensional matrix L has N vib vibrational eigenvectors l µ as column vectors which are orthonormal to each other. The vibrational frequencies ω µ can be connected with eigenvalue λ µ according to λ µ = 4π 2 c 2 ω 2 µ where c is the speed of light. Equation (2) can be rewritten in terms of normal coordinates Q [52,70] as where K is the Hessian matrix expressed in normal coordinates Q with dimension (N vib ×N vib ) and L T is the transpose of L.
Konkoli and Cremer defined a local vibrational mode via the leading parameter principle [57], which states that a local vibration is initiated by an associated internal coordinate q n via its infinitesimal change. Only the masses of the atoms involved in this internal coordinate q n are kept, the masses of all other atoms are assigned a zero value. As a consequence, the other atoms of the molecule can effortlessly follow the local vibration led by q n as a collection of massless points. The internal coordinate q n can be associated with the Cartesian coordinates x of the molecule via the Wilson B-matrix [52] which collects the partial derivatives of q n with regard to Cartesian coordinates x b n = ∂q n ∂x (4) The local mode vector a n associated with internal coordinate q n is then defined by with d n = b n L [70]. The local mode vector a n , a column vector of length N vib can be transformed into Cartesian coordinates via Equation (6) To each local mode a n , local mode properties can be assigned, such as a local mode force constant, frequency and mass [57]. The local mode force constant k a n of internal coordinate q n is obtained with The local mode force constant k a n was also named adiabatic force constant, where a (adiabatic) as the superscript means "relaxed" and n as the subscript corresponds to the internal coordinate q n leading this local vibration [57].
The local mode mass m a n of mode n is given by where G n,n is the n-th diagonal element of the Wilson G matrix [52,70]. Local mode force constant and mass are needed to determine the local mode frequency ω a n (ω a n ) 2 = 1 4π 2 c 2 k a n G n,n Zou and co-workers demonstrated that, for a complete non-redundant set of N vib local modes, there exists a one-to-one relationship between local and normal vibrational modes that can be verified by an adiabatic connection scheme (ACS), providing the physical fundament for the local vibrational modes [61]. The reason why a complete set of non-redundant parameter set is required in such relation is because this set of local vibrational modes can span the same internal vibration space spanned by N vib normal modes [71][72][73]. In addition, this one-to-one correspondence between the local and normal vibrational modes forms the basis for the CNM method leading to a detailed analysis of a vibrational spectrum and in this way decoding a wealth of information embedded in the spectrum [65,74]. It is also important to note that this analysis can be applied to both calculated and experimentally determined fundamental vibrational frequencies [75].
According to the CNM method [59], any normal vibrational mode l µ can be decomposed into local mode contributions from a non-redundant set of N vib local vibrational modes by calculating the overlap between each local mode vector a x n with this normal mode vector l µ as S nµ according to Equation (10) where (a, b) is a short notation for the scalar product of two vectors of a and b O ij is element within the metric matrix O. In this work, we used the force constant matrix f x as the metric (i.e., O = f x ) to include the influence from the electronic structure. The calculation of S nµ in Equation (10) was simplified via the following steps. If we consider the complete set of N vib normal modes collected in L and the non-redundant set of N vib local modes collected in A x , Equation (10) can be written as where where D = BL collects d n as row vectors, then and Then, Equation (10) can be re-written as S nµ = D 2 nµ k a n K µµ (26) where K µµ is the µ-th diagonal element of matrix K. As long as the molecular system is at a local minimum on the potential energy surface, both k a n and K µµ are positive thus leading to a positive value of S nµ . It can be easily proven that S nµ is independent of the prefactor inside an internal coordinate as a linear combination of a few basic internal coordinates (e.g., S nµ stays the same whatever nonzero value p takes if the internal coordinate q n is defined by q n = p · (q a ±q b )).
Therefore, the contribution of local mode a n to the normal mode l µ can be calculated by In order to evaluate the degree to which a normal vibrational mode l µ involving the stretching of the probe chemical bond is decoupled from all other internal coordinates, a complete non-redundant set of internal coordinates including the one for the probe bond (denoted as q 1 ) has to be constructed and the corresponding local mode vectors a n need to be determined as well as their overlap with the target normal vibrational mode l µ . In this way, one can quantify the extent to which the normal mode has predominantly probe bond local stretching character and determine the actual percentage via C 1µ is a number ranging from 0 to 1. A large C 1µ value indicates that normal vibration mode l µ has predominantly the local vibration led by internal coordinate q 1 and a value of 1 would identify mode l µ as a 100% local probe bond stretching vibration. For simplicity, we coined the term "performance score" being defined as the C 1µ percentage value to evaluate and compare different vibrational Stark effect probes in the remainder of this work. Theoretically speaking, any complete and non-redundant internal coordinate parameter set (including the probe bond) could be used for normal mode decomposition; however, in order to better accommodate the rocking (asymmetric bending) vibration in formaldehyde-like topology (see Figure 1), we employed an antisymmetric combination (denoted as δ) of two angles containing the probe bond (e.g., C=O) when analyzing these probe molecules. In the following, the procedure of calculating the performance score for a VSE probe is demonstrated for formaldehyde (1-1) as an example. This molecule contains four atoms and has as such six normal vibrational modes. We selected as non-redundant set of six internal coordinates the three bond lengths, two bond angles (including δ), and one dihedral angle (τ) shown in Figure 2, which also shows the decomposition of each of the six normal vibrational modes into local mode components in the form of a bar diagram for both experimental and calculated vibrational frequencies.
The two H−C−O angles were anti-symmetrically combined to angle δ. The probe bond vibration is colored in yellow. The normal mode decomposition of formaldehyde (see Figure 2) clearly identifies normal mode 4 as probe bond stretching with the largest local C=O stretching contribution, i.e., C nµ value of 0.8816 and a contribution of 12% from the local H−C−H angle bending mode. As such, the performance score of formaldehyde as the VSE probe is 88.
The potential energy distribution (PED) method [77][78][79][80][81][82] widely used in vibrational spectroscopy has been applied in some scattered investigation to analyze vibrational Stark effect probes [83]. PED takes a step back and it is based on the idea that the potential energy can be expressed as a power series in terms of normal coordinates Q, which can be further decomposed into internal coordinates q n . In contrast, the CNM method as part of the Konkoli-Cremer local vibrational mode theory is directly based on vibrational spectroscopy and local modes can be smoothly transitioned to normal modes via adiabatic connection scheme [67]. In this sense, CNM analysis is superior than PED analysis as the former has better physical fundament in terms of vibrational spectroscopy [58][59][60]84].
One may raise concern on the suitability of calculating the performance score based on an isolated probe molecule model as in real application scenarios a probe molecule may have non-covalent interactions with surrounding molecules. Therefore, we constructed a model system of formaldehyde which is hydrogen bonded with a water molecule as shown in Figure 3. The target normal vibration of C=O stretching consists of the local vibrations of 87.2% R(C1=O4) and 12.6% α(H2−C1−H3). As such, in the hydrogen bonded complex, the performance score for the formaldehyde molecule is 87.2, only somewhat smaller than the score of 88.2 for the single molecule shown in Figure 2 by 1.0 unit. Such minor difference in performance score can be safely ignored. Therefore, in the remainder of this work, we stick to the single molecule approach for the calculation of the performance score of the probe molecules.

Computational Details
Geometry optimizations and Hessian evaluations for all probe molecules investigated in this work was carried out in the gas phase using the Gaussian 16 package [85]. All molecules were optimized using the ωB97XD density functional [86] with Dunning's aug-cc-pVDZ basis set [87][88][89] except molecules containing a porphyrin group, which were calculated at the M06L/def2TZVP level [90][91][92].
The decomposition of the normal vibrational modes into local mode contributions was carried out with the COLOGNE2019 package [93]. In order to test the sensitivity of the decomposition results with regard to the density functional employed in this work, all calculations were repeated at the M06-2X/aug-cc-pVDZ level [94] except for molecules containing a porphyrin group.

Results and Discussion
In this section, (i) we evaluate the performance of commonly used VSE probes in comparison with a number of potential vibrational probe candidates, using the CNM method. These polyatomic probe molecules contain C=O, C≡N, S=O or other chemical bonds, whose stretching vibrations are considered as decoupled from other local vibrational modes. (ii) We analyze how the atomic masses influence the performance of selected VSE probes and discuss the feasibility of improving a vibrational probe by isotopic substitution. (iii) In addition, we analyze how the performance score of a VSE probe depends on the density functional used for the calculation. (iv) Practical suggestions on the ideal VSE probes to experimentalists are made with related physicochemical properties (e.g., solubility and reactivity).    Below each structure, a label of n-m is given in blue; n denotes the group number (reflecting the type of VSE probe bond and m is the number of the molecule within this group. Vibrational probe labels with superscripts are taken from the literature (a [95], b [5,96], c [96], d [5,7,25,28], e [5,13,97], f [13], g [13,26], h [14,15,98]). The corresponding performance score as VSE probe is given in purple. The number as superscript in the 2D structure refers to atom index in the molecule. The superscripts are shown only for the atoms whose associated local modes participate in the C=O or C≡O normal mode with more than 5% contribution. Table 1 shows the decomposition of target normal mode of probe bond stretching into corresponding local modes where all contributions greater than 5% are shown. Formaldehyde (1-1) (described in more detail in the Methodology Section) has relatively simple structure and its performance score is 88. As revealed by Equation (26), the overlap S nµ is composed of different terms, k a n characterizes the pure electronic effect, whereas D nµ and K µµ depend on the normal modes, which involve the atomic mass, geometry, and electronic effect. This implies that the performance score is a result of multiple factors. As McKean has shown in his work that the isotopic substitution could result in isolated (i.e., local) modes of CH bond stretching in -CD 2 H groups [99], it would be interesting to see if similar strategy could lead to better performed VSE probes. In order to do so, we calculated the performance score landscape for 1-1 as a function of the atomic masses as shown in Figure 5.

Group 4: Other probes
Mol. a Local mode contributions b 74.6% C1-Na14, 11.5% (C2-C1-C10, C6-C1-C10), 5.8% C2-C1-C6 4-8 99.9% H35-Si34 4-9 97.9% H2-S1 4-10 99.9% H2-S1 4-11 85  Changing the atomic masses of carbon and oxygen atoms (the two atoms of the VSE probe bond) from the most abundant isotopes ( 12 C and 16 O) to higher hypothetical isotopes results only in a trivial change (ca. -3.0) of the performance score. By replacing the hydrogen atoms with deuterium, the performance score drops down to 86.07. The isotope effect seems to play a complex role on the performance score and is not directly intuitive for us to decide which atomic mass to change.
Carbonyl fluoride (1-2), carbonyl chloride (1-3) and carbonyl bromide (1-4) have a similar structure as 1-1 except that the hydrogen atoms are replaced with more electronegative halogen atoms X. The performance score of 83 for 1-2 is the lowest in this series, gradually increasing for the higher homologues, i.e., 93 for 1-3 and 95 for 1-4. With the increase in the atomic number in the series F, Cl and Br, the electronegativity of X decreases, the atomic mass increases, and the C−X bond becomes longer and weaker, which should lead to less coupling with the C=O stretching. This is obviously in line with the calculated performance scores. However, it is oversimplified to conclude that the difference in the performance score for these three molecules can be completely attributed to electronic effects unless the mass effect can be eliminated. Therefore, we performed model calculations using the same atomic mass in order to extract the electronic effect, as shown in Table 2. Table 2. Comparison of the performance scores of molecules X 2 CO, 1-1-1-4 using the atomic masses X = H, F, Cl, and Br. Each column represents the same electronic structure with different masses and a row represents same atomic mass with different electronic structures.  Table 2 shows how the performance scores of these four molecules are affected by the atomic mass of X. When the atomic mass of X is kept the same for all of the four molecules, the performance scores show a consistent pattern: 1-4 > 1-3 > 1-1 > 1-2. As expected, the lowest performance score is found for 1-2, (i.e., shortest and strongest C−X bonds leading to a large mode-mode coupling with the C=O target bond). The highest score is found for 1-4 as its longest and weakest C−X bonds leading to only moderate coupling with the C=O target bond. The electronegative F atoms in 1-2 attract electron and lead to stronger C−X bonds. The joint contribution from the local stretchings of two C-F bonds to the target normal vibration is 10.4%. The less electronegative Br atom in 1-2 is less suited to attract electrons, which is linked to its larger atomic radius leading to longer C−Br bonds, whose local vibration contribution to the target normal vibration is less than 5% in total. Overall, this result reveals that the performance score differences are mainly caused by electronic structure differences and not by a significant mass effect.

Mass of X (amu) H 2 CO (1-1) F 2 CO (1-2) Cl 2 CO (1-3) Br 2 CO (1-4)
Compared to substituent effects, isotopic effects are a less effective choice in tweaking the performance score. One also has to consider that isotopic substitution is often experimentally demanding and time consuming. Therefore, in the remainder of this work, we will focus on substituent effects. The carbonyl group can be incorporated into a ring structure. When the C=O bond is attached to cycloalkyl groups (1-6, 1-7, 1-8 and 1-9), the performance score is relatively low ( 83) due to the large contamination from the local vibrations of the adjacent C−C single bonds. For example, in cyclopropanone (1-6), two local C−C vibrations adjacent to the C=O bond contribute jointly 19% to the target normal vibration mode dominated by the local C=O stretching. Meanwhile, the performance score remains almost the same with increasing ring size. However, the second largest contribution from adjascent C−C−C bond angle to the C=O stretching normal mode in 1-8 (8%) and 1-9 (5%) decreases, which can be related to decreasing ring strain. In 3,4-Dimethyl-2-cyclohexen-1-one (1-12) the C=O bond is attached to a six-membered ring. 1-12 is a simplified motif of an inhibitor which has been frequently used to measure the electric field inside the binding pocket of ketosteroid isomerase [5,22,25]. Its relatively low performance score of 78 compares to the scores of cyclopentanone (1-8) and cyclohexanone (1-9) with similar ring structures. A slight difference in the performance score can be attributed to the extended electron density (i.e., π-conjugation) in 1-12. The CNM analysis shows that the local C−C−C angle bending contributes to 7%.
Acetophenone (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13) has been used to measure the electric field inside various solvents [13,97]. Its performance score is 80, which is slightly lower than that of acetone. While most contamination of the target normal mode of 1-13 is caused by the adjacent C−C−C angle bending (6%), the CC π bonds of the phenyl ring conjugated with the C=O bond contribute 3% to the target normal vibration collectively. Replacement of the methyl group in 1-13 with another phenyl group leading to benzophenone (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14) does not impact the performance score significantly. However, replacing the methyl group of 1-13 with a methoxy group leading to methyl benzoate  decreases the performance score to 76. The CNM analysis identifies for 1-21 the adjacent C−C−O angle bending and C−C local bond stretching vibrations as important contributors to the target normal vibration with 7% and 6%, respectively. Similar contribution patterns to the target normal mode are found for methyl acetate (1-20) with a performance score of 82. In ethyl thioacetate , the performance score is raised to 87. The major contamination arises from the adjacent C−C bond. While the local C−C bond stretching contribution is 5% and the C−S bond contribution is less than 5%.
In 1,4-pentadiyn-3-one (1-15) and oxo-malononitrile (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16), triple bonds are connected to the carbonyl group increasing the performance scores to 89 and 90, respectively. It could be expected that the π conjugation between the triple bond and the C=O bond might lower the performance score as in 1-10 and 1-11. However, due to the linear C−C≡C/N bond angle arising from sp-hybridized carbon, the target normal vibration is then dominated by the local vibration of the C=O bond stretching. In addition, the local C−C single bond stretching and C−C−C angle bending vibrations contribute with less than 5% to the target normal mode.
Carbon monoxide can serve as a diatomic VSE probe [100]. When coordinated to the iron atom in iron porphyrin , the C≡O ligand can be used to characterize the electric field in enzymes [7,10]. According to our analysis, 1-24 has a high performance score of 93. The local Fe−C bond contributes with 7% to the target normal vibration. When an additional imidazole group is coordinated to the iron atom as in 1-25, the performance score increases marginally to 94 as the contribution from the local vibration of C−Fe bond decreases to 5%.
We have observed that the electronic effects play a critical role in determining the performance score compared to the mass effects. To be more specific, increasing ring strain by reducing the ring size and enhancing the π electron delocalization can lower the score.

Groups 2 and 3: C≡N and S=O Probes
In the following, we will discuss the performance scores for 17 molecules with a C≡N probe bond and 6 molecules with a S=O probe bond shown in Figure 6. Nitriles are one of the most commonly used VSE probes [98,101]. According to the CNM analysis, some S=O probes seem to perform even better.   [10,18,24,27], c [24], d [18,50], e [5,10,18,26,27,50], f [18], g [26,50,51], h [102], i [20]). The corresponding performance score as VSE probe is given in purple. The number in superscript refers to atom index in the molecule. The superscripts are shown only for the atoms whose associated local modes participate in the S=O or C≡N normal mode with more than 5% contribution.
Thionyl difluoride (3-1) has a highly localized S=O bond with almost negligible contamination from other local modes. The same applies to thionyl dichloride (3-3) and thioaldehyde (3)(4). In dimethyl sulfoxide (3)(4)(5) and dimethyl sulfite (3)(4)(5)(6), one important contamination comes from the pyramidalization mode, i.e., movement of S atom in the normal direction of the plane formed by adjacent C−O−C plane (5%) and O−O−O plane (6%), respectively. The two equatorial S−F bonds in thionyl tetrafluoride  lead to significant contamination in the normal mode of S=O bond stretching.
We observe that the C≡N and S=O bonds give rise to many potentially good candidate for VSE probes with high performance scores.

Group 4: Vibrational Probes with Miscellaneous Bonds
Besides the C=O, C≡N and S=O probes, other types of chemical bonds have been applied [18] or could be used for VSS as shown in Figure 7.
As demonstrated above, the nitrile group (-C≡N) leads to VSE probes with high performance scores. Therefore, we tested a series of potential VSE candidates with triple bonds. The Si≡N probe bond of methylnitrilosilane (4-1) has a satisfactory performance score of 92. The local Si−C stretching vibration contaminates the target Si≡N normal vibration with 7%. The related phenylnitrilosilane (4-2) shows a lower performance score of 87 due to larger contamination from the Si−C bond (9%). However, when the nitrogen atom in 4-1 is replaced with phosphorus (4-5), the performance score drops to 61. The major contamination comes from the adjascent C−Si bond of 39%. The C≡S triple bond (4-11) performs well with a score of 86, with a contamination from the local H−O−S angle bending mode (11%).
If the lithium atom in 4-3 is replaced with sodium leading to tert-butylsodium (4-7), the performance score remains the same with contaminations from the three local C−C−C angle bending vibrations. However, if the lithium atom in 4-3 is replaced with hydrogen to isobutane (4-6), the performance score decreases.
The number in superscript refers to the atom index in the molecule. The superscripts are shown only for the atoms whose associated local modes participate in the normal mode of probe bond stretching with more than 5% contribution.
We observe that the Si≡N bond shows similar characteristics as C≡N bond in terms of the performance score. However, replacing N to P significantly lowers the performance score. Including a combination of low and high atomic mass of atoms for a bond is another way to increase the performance score.

Sensitivity of Performance Score to Density Functional
As the performance score for the probe molecules based on the Characterization of Normal Mode (CNM) method is highly dependent on the quality of the Hessian matrix, and it has been tested that the vibrational frequencies calculated for the same molecule by different density functionals [103] can have differences up to 100 cm −1 , one might argue that the performance scores calculated with CNM could be dependent on the employed density functional.
In order to test the sensitivity of performance score to the selected density functional, we repeated all the calculations including geometry optimization and Hessian evaluation with Truhlar's Minnesota functional M06-2X, which is believed to perform well for organic compounds [94]. Then, the performance scores calculated by ωB97XD and M06-2X are compared in Figure 8.  Figure 8 shows that the performance scores calculated by M06-2X and ωB97XD are quite consistent, with an average difference of 1.2 for each probe. Probe molecules with relatively high performance scores above 85 are mostly on the line of y = x, with an average difference of 0.5 excluding the outlier thionyl difluoride. Probes with performance scores below 85 have an average difference of 1.7. The large deviations in thionyl difluoride is due to its different geometry when optimized with M06-2X (non-planar) compared to ωB97XD (planar). The other four outliers (dimethyl sulfite, acetophenone, methyl benzoate, and penta-1,4-dien-3-one) can be attributed to the difference in describing the π conjugation in these four molecules with ωB97XD and M06-2X. This result tells us that the performance score obtained in this work is insensitive to the selected density functional employed in calculation. The performance score is especially reliable for the probes when it is relatively high.

Suggestion on Ideal Vibrational Probes
From a systematic evaluation of the polyatomic vibrational Stark effect probes with our CNM method shown above, we pick out the molecules with scores higher than 85 as ideal probes to be recommended to vibrational Stark spectroscopy practitioners. Note that we do not deny that probes with lower performance score can be used as VSE probes; however, choosing well-performed probes will significantly reduce the chances of errors caused by the contamination of other local modes.
In Table 3, we have listed the performance score and calculated frequency for the target normal vibration for each of the ideal probes. As in a real "wet lab" when the vibrational probes are to be put into use, we have to consider the solubility of the probes into a specific solvent and its reactivity; therefore, we try to collect all related information for these probes as a reference for experimentalists. The bold label refers to novel molecular probes introduced in this work. b Physiochemical properties are taken from PubChem database [104] and literature [105]. c The vibrational frequencies are computed at M06-2X/aug-cc-pVDZ level of theory and then scaled by an empirical factor of 0.9500 [106].

Conclusions
In this work, we have employed the characterization of normal mode (CNM) method from Konkoli-Cremer local vibrational mode theory to evaluate more than 68 different vibrational Stark effect probes by quantifying to what extent the probing normal vibrational mode is indeed localized to the local stretching of the probe bond. The quantification was realized by using a performance score in the range of 0∼100 for each vibrational probe for comparison. The probe bonds investigated in this work include C=O/C≡O, C≡N, S=O, Si≡N, C−Li, Si≡P, C−Na, S≡C, N=N + =N -, S−H, N=O, P=O and C−H. In general, we found probe molecules with double or triple bond tend to score higher than those with single probe bond. However, in several cases of probe molecules with single probe bond (e.g., S−H in a tetrahedral geometry) the performance score can also be high.
We have shown that isotopic substitution for a vibrational probe is a much less efficient approach for obtaining better performed probe candidates compared with tweaking the structure via chemical modification. The validation of the calculated performance score results with different density functionals consolidates the reliability of the CNM method employed in this work.
It is important to note that we did not attempt to explore exhaustively the chemical space by trying all possible combinations of different elements and functional groups for a vibrational probe molecule. However, we have included those representative probes as a guidance for vibrational Stark spectroscopy practitioners. As it has been tested in this work, we can expect that the performance score of a probe molecule will only have a slight change if the chemical modification is distant from the probe bond. Furthermore, we provide an online service (https://vse-server.github.io/) for interested readers to evaluate their vibrational probe candidate molecules with our CNM method.