The chemical bond in external electric fields: Energies, geometries, and vibrational Stark shifts of diatomic molecules

The interaction of molecules with external fields, static or time dependent, can significantly alter the potential energy surfaces (PES) governing chemical reactions, that is, the breaking and making of chemical bonds.^1–15 The simplest PES is that of the dissociation of a diatomic molecule since in this case the potential energy curve depends on only one variable, the internuclear separation. This paper reports a systematic computational study of the effects of external uniform static electric fields (referred to as “fields” in this paper) on diatomics with the goal of relating the field-induced changes in the PES with the accompanying field-induced changes in other molecular, atomic, and bond properties.

Two opposite charges of ±0.5 e⁻ separated by 10 Å in an enzyme active site, for example, generate fields ∼10⁹ V m⁻¹ at their center. These fields can induce observable vibrational Stark shifts in the IR frequencies of a host molecule trapped into the active site and can be exploited to probe the strength and direction of the electric fields in its active site, the host here acting as a “local reporter of its electrostatic environment.”¹⁶ Examples of “local reporter” molecules include carbon monoxide (C=O) attached to the Fe of the porphyrin ring in myoglobin,¹⁷ and the nitrile- (–C≡N) containing substrates in the active site of human aldose reductase enzyme (hALR2).¹⁶ From the changes in the vibrational Stark shifts of the nitrile group, Boxer et al. were able to estimate a change of 10⁸ V m⁻¹ in the electric field accompanying a single point mutation at the active site,¹⁶ but interference of the effect of hydrogen bonding with the Stark effect in enzyme active sites can complicate the interpretation of spectra.^18,19 Zhang et al.²⁰ and Gascón et al.²¹ have recently demonstrated that an accurate modeling of the electric charges in the cavity of hALR2 is essential for the accurate prediction of the mutation-induced Stark shift in the vibrational frequencies of the nitrile group.

Molecules can also be exposed to fields from sources external to the organism, whether applied in an experiment or in the crystal environment of the molecule, for example. In a theoretical study of a set of 20 common nonlinear optical materials, Spackman, Munshi, and Jayatilaka have shown that the crystal fields to which molecules are typically exposed to are in the range ∼(1–3) × 10⁹ V m⁻¹ but can reach ∼10¹⁰ V m⁻¹.²² Fields of strengths ∼10⁹ V m⁻¹ have been shown to quadruple the rate constant of double proton transfer reactions,¹⁵ can be used as “tweezers” to pick a reaction channel when more than one channel link a set of reactants to different sets of products,²³ and are found between the tip and the sample in a scanning tunneling microscope (STM) under normal operating conditions.²⁴ Fields of 10⁸ V m⁻¹ accelerate photosynthetic reactions by an order of magnitude,^25–27 and significantly decrease the enzymatic activity of cytochrome c oxidase.²⁸ 

The maximal strength of macroscopic fields applied to a crystal during an X-ray crystallographic diffraction is around 10⁷ V m⁻¹.²⁹ The field-induced changes in diffraction patterns have been observed for sometime,^30–35 and require for their observation the use of intense synchrotron sources, low temperatures, and efficient data collection methodologies. Tsirelson, Gorfman, and Pietsch^29,36 demonstrated that the shift in atomic positions in the unit cell is ∼100 times more significant than the polarization of the atomic electronic clouds in determining the field-induced perturbation to the diffraction pattern.²⁹ 

In summary, strong external electric fields are common in the environment of molecules and can have drastic effects on their properties and reactivity. This computational study aims at gaining insight into the effects of external fields on the basic properties of the simplest covalent chemical bonds, those in diatomics, with varying degrees of polarity. Each molecule in a set of first to third row homo- and hetero-nuclear diatomics is subjected to fields of different strengths along its C_∞-axis (in two opposite directions in the case of the heterodiatomics). The studied properties include: total energies (E), dipole moments, bond lengths (R), force constants (k), and vibrational frequencies (ν). Field effects on other bond properties, including topological properties, will be discussed elsewhere.

The set of diatomic molecules includes five homo-nuclear diatomic molecules (H₂, N₂, O₂, F₂, and Cl₂) and four hetero-diatomic molecules (^+1.22C–O^−1.22, ^+0.75H–F^−0.75, ^+0.44N–O^−0.44, and ^+0.26H–Cl^−0.26) that cover a range of bond polarity.

Only the ground electronic states of the diatomics were considered: |${}^1\Sigma _g^ + $|$Σg+1$ for closed-shell homo-nuclear dimers, ¹Σ⁺ for closed-shell hetero-nuclear dimers, ²Π_r in the case of NO, and |${}^3\Sigma _g^ - $|$Σg−3$ in the case of O₂.

To increase the confidence in the results, calculations were conducted using two significantly different underlying electronic structure methods, namely, (i) quadratic configuration interaction with single and double excitations (QCISD),³⁷ and (ii) density functional theory (DFT) Becke's hybrid exchange functional³⁸ and Lee-Yang-Parr correlation functional³⁹ (B3LYP). To keep the interpretation of the results simple, a single (large) polarized split-valence basis set augmented with diffuse functions on all atoms [6-311++G(3df,2pd)] was employed in all calculations. The unrestricted formulations of QCISD and of B3LYP have been used for open-shell molecules, atoms, and ions.

Molecular geometries were optimized with and without external fields within convergence thresholds of 2 × 10⁻⁶ and 1 × 10⁻⁶ hartree/bohr for the residual maximum forces and root-mean square forces on the nuclei, respectively.

Electronic structure calculations, geometry optimizations, and frequency calculations were all conducted using the Gaussian 09 program.⁴⁰ All reported frequencies are those obtained directly from the calculation (unscaled).⁴¹ 

Both levels of theory QCISD and B3LYP generally yield results that agree both absolutely and in trends. Thus, unless mentioned otherwise, the discussion and quoted results are primarily those of the QCISD calculations. The full set of QCISD and DFT-B3LYP results are provided in the supplementary material.⁴² 

There is one exception to the close agreement of QCISD and B3LYP, and this is in the trends in field effects on the frequencies and force constants of NO for which QCISD and B3LYP give widely different trends and absolute values. To determine which of these two levels of theory is more credible in this case, calculations were conducted at four additional levels of theory: coupled clusters^43–45 with single and double excitations (CCSD), Møller-Plesset second order perturbation theory⁴⁶ with full excitation including core electrons, Hartree-Fock (HF) self consistent field (SCF) level,⁴⁷ and the DFT functional mPW1PW91.⁴⁸ 

The atomic unit (a.u.) of electric field strength is defined as the field strength at the first Bohr orbit in the hydrogen atom, |${e\,/\,{4\pi \varepsilon _0 a_0^2 }}$|$e/4πɛ0a02$ ≈ 5.14 × 10¹¹ V m⁻¹. In this work, electric fields from 5.14 × 10⁸ to 5.14 × 10¹⁰ V m⁻¹ (=0.001 to 0.1 a.u.) were applied along the molecular C_∞-axis (the z-axis of the coordinate system, Figure 1), since fields aligned with the bonds maximally affect the vibrational frequencies.⁴⁹ Any field directed obliquely on the internuclear axis can always be decomposed into an (anti)parallel and a perpendicular component.

Disregarding orientation dependence of field effects on non-spherically symmetric electronic systems, a rough estimate of the order of magnitude of the rate of electron tunneling ionization (ω) from a molecule in an external electric field can be obtained from the well-known formula:^50,51

where |$|{\sfb{E}}| = {{\hbox{\sf\itshape E}}}$|$|E|=E$ is the electric field strength, IP is the vertical ionization potential, both quantities in a.u., and ω₀ is the atomic unit of frequency (4 × 10¹⁶ s⁻¹). Tunneling ionization rates calculated from Eq. (1) using the experimental IPs of the nine diatomics are collected in Table S1 of the supplementary material.⁴² 

Chemical reactions generally happen on a time scale of ∼10–10⁴ fs.⁵² One can take the reciprocal of 10⁻⁴ fs⁻¹ as the maximal acceptable average tunneling ionization rate for a molecule in the field. An arbitrary order of magnitude threshold of around ω ∼ 10¹¹ s⁻¹ can thus be considered reasonable. Given these considerations, the values listed in Table S1 suggest a field strength of 5.14 × 10⁹ V m⁻¹ (0.01 a.u.) as an upper ionization threshold with the exception of alkali metal dimers Li₂ and Na₂ which are excluded from this study in view of their low ionization potentials. The results described below show that fields of up to 0.02 a.u. do not induce changes in trends and are, thus, also kept as an upper extreme limiting case, and given the very approximate nature of Eq. (1) especially for non-spherically symmetric molecules.

Hetero-nuclear molecules are subjected to fields in the two opposite directions: parallel and antiparallel to the C_∞-axis. The molecules are oriented in the coordinate system by placing the least electronegative atom (Pauling's scale)⁵³ at the origin and the second atom along the positive z-axis. The Cartesian coordinate system and orientations of the diatomics are displayed in Figure 1. Whenever a field is mentioned without its magnitude what is being referred to is the field with the strongest magnitude. Hence, |${\sf{E}}_{+}$|$E+$ or |$\overrightarrow {{\sf{E}}_ + } $|$E+⃗$ mean a field of magnitude 1.03 × 10¹⁰ V m⁻¹ = 2.0 × 10⁻² a.u. oriented to point in the positive direction of the z-axis shown in Figure 1 while |${\sf{E}}_ -$|$E−$ or |$\overleftarrow {{\sf{E}}_ - } $|$E−←$ is a field of the same strength but pointing in the opposite direction.

In this paper, the dipole moment vector is directed according to the “physicists convention,”^54,55 i.e., it originates at the negative pole and points to the positive pole, and can be denoted by |$^ + \overleftarrow {\bm \mu } ^ -$|$μ←−+$ where μ symbolizes the dipole moment (e.g., ^δ+H ← Cl^δ−). On the other hand, all electric fields originate at positive charges (sources) and end at negative charges (sinks). (Note that the Gaussian 09 program orients electric fields in the reverse direction.) All fields in this paper have the same direction and sign as the z-axis along which they are aligned (Figure 1) and are uniform (as in an infinite parallel-plate capacitor). These conventions ensure that a dipole is in its most stable direction (the direction that minimizes its energy in the field) when parallel to the external field since |$E_{\scriptsize{\sfb {E}}} = - {\bm \mu} \cdot {\sfb{E}}.$|$EE=−μ·E.$⁵⁴ With these conventions, a stabilizing (energy lowering) orientation of the dipole is parallel to that of the external field and can be given the symbol , where the arrow below the brackets is the direction of the field. On the other hand, symbolizes an antiparallel and destabilizing (energy raising) orientation of the dipole with respect to the external field.

All studied molecules are oriented in the coordinate system so that one of the two atoms is at the origin and the other lying in the positive side of the z-axis. In the case of heterodiatomics, the least electronegative atom is the one placed at the origin as shown in Figure 1. This convention is independent of the direction of the permanent molecular dipole. Thus, in two instances the dipole points at the origin (H ← F and H ← Cl) and in the other two it points away from the origin (C → O and N → O) as depicted in the figure [the arrows between the atomic symbols indicate the direction of the permanent (field-free) molecular dipole which changes in magnitude in external fields and can also change direction, as described below].

We remind the reader that while the permanent molecular dipole moment generally points in the direction expected on the basis of the different electronegativities of the bonded atoms due to the inter-atomic transfer of charge as in ^+0.75H ← F^−0.75 or ^+0.26H ← Cl^−0.26, there are known exceptions where the positive end of the dipole points at the atom bearing a net negative charge such as ^+1.22C → O^−1.22 and ^+0.44N → O^−0.44. The reasons for this unexpected direction of the dipole moment (unexpected on the basis of electronic charge transfer) has been discussed in Ref. 56 and is due to a large and opposite atomic polarization term,^57–59 that is, a distortion of the charge cloud of an atom in a molecule from spherical symmetry that can cancel and even overwhelm and dominate the charge transfer dipole. The interplay of the charge transfer dipole and the atomic polarization dipole has been recently studied over the entire potential energy surface of the reactive collision of halogens with methane.⁶⁰ 

Table I summarizes the molecular properties in the field free case and in the presence of the strongest field strength (0.02 a.u. = 1.03 × 10¹⁰ V m⁻¹), in both directions in the case of hetero diatomic molecules. The properties investigated as functions of the applied field strength and direction in this work include the total energy (E), the molecular dipole moment (⁠|$\bm\mu$|$μ$⁠), the equilibrium bond length or internuclear separation (R), the force constant (k), and the harmonic vibrational frequency (ν). The table is organized so as to list, for every molecular property, three values: The field-free value in the middle row flanked by the values under the strongest studied field strength in the two opposite directions. Since for homo-nuclear diatomics the two field directions are equivalent, the fields are assigned a positive sign.

The expansion of the energy [|${\displaystyle{E({\sfb{E}},r)\equiv E_{{\sfb{\scriptsize E}}}}}$|$E(E,r)≡EE$] of a diatomic molecule in an external homogenous electric field (E) as a power series takes the form⁶¹ 

where E₀, μ₀, α_//0, and α_⊥0 are the field-free energy, (permanent) dipole moment, and parallel and perpendicular polarizability tensor components, respectively, and where the functional dependence of all these quantities on the internuclear separation (r) is emphasized; |${\hbox{\sf\itshape E}} = |{\sfb{E}}|$|$E=|E|$ is the magnitude of the external electric field that makes an angle θ with the permanent molecular dipole; and the last term collects higher order terms.

Similarly, the ith cartesian component of the total molecular dipole moment under the field, |$\mu _i ({\sfb{E}},r)$|$μi(E,r)$⁠, may be expressed at the sum of the permanent (field free) component, μ_i0(r), plus the induced dipole represented by the remaining terms in Eq. (3):

Plots of the molecular dipole moments |$\bm\mu$|$μ$ of the nine studied diatomic molecules as functions of the external fields are displayed in Figure 2. From the figure it is clear that, within the range of the electric field strengths considered here, the molecular dipole moment |$\bm\mu$|$μ$ of each molecule is directly proportional to the applied field strength over the entire range of field strengths and orientations. The Pearson linear regression coefficient relating each molecular dipole moment to the external field is unity to three decimals, and this includes the doubled range in the case of heterodiatomics due to the flipping of the sign of the field from parallel to antiparallel. Because of this proportionality of the dipole moment and the electric field strength, the higher terms in Eq. (2) can be ignored.⁶² 

Further, only fields that are co-linear with the molecular axis (parallel θ = 0, or antiparallel θ = π radians) are considered, and hence Eq. (2) can be simplified and rearranged to define the field-induced change in the energy (stabilization or destabilization), ΔE, as

where the dipole term assumes a negative sign for parallel (stabilizing, energy lowering) fields and a positive sign for antiparallel field. From now on, the subscript “0” will be dropped from the symbols for dipole moment and polarizability when it is clear from the context that the parameter of interest is the field-free parameter.

Table II lists the experimental and calculated polarizabilities and permanent dipole moments of the molecules considered in this work, sorted in the order of increasing polarizability within each of the two groups. The selected molecules cover a range of polarizabilities and of dipole moments. The homo-nuclear diatomics molecular set includes H₂ for which α ≈ 0.8 ų and up to the highly polarizable Cl₂ with α ≈ 4.6 ų. On the other hand, the hetero-nuclear diatomics include combinations of polarizability and permanent dipole moment. Thus, HF has a significant dipole moment of 1.8 debye and a polarizability as low as that of H₂; CO and NO have small permanent dipole moments (∼0.1 debye) and sizeable polarizabilities (∼1.8 ų); and finally HCl has both a sizable dipole moment (∼1.1 debye) and a sizable polarizability (∼2.5 ų). The listing in the table also shows a good agreement between measured and calculated parameters.

The change in the total energy as a function of the applied field is displayed for the nine diatomics in Figure 3 where plots for homo-nuclear molecules are collected at the top of the figure and those for the hetero-nuclear diatomics at the bottom. Equation (4) and the relative values of dipole moment and polarizability for each molecule lead one to anticipate its response to an external field as long as the response of the dipole moment remains linear (which is the case for the nine molecules over the range of fields considered). For homo-nuclear diatomics, the dipolar term in Eq. (4) vanishes and hence: (a) all fields are energy lowering to all molecules due to the nonlinear term which is always negative and (b) the energy lowering is proportional to the square of the field strength.

From the plots in Figure 3 (top), the molecule that is least affected by the field is H₂, the one with the smallest polarizability among the homo-nuclear diatomics, while the molecule with the largest polarizability, Cl₂, is the most stabilized by the field. The curves representing the response of the remaining homo-nuclear diatomics fall into their relative places according to their respective polarizabilities. A glance at Table I shows that at the highest studied field (1.03 × 10¹⁰ V m⁻¹) a sorting of the homonuclear molecules X₂ in order of increasing polarizability parallels a corresponding sorting in the magnitude of the stabilization energy. Thus, (α_// in ų/ΔE in eV): H₂ (1.0176/−0.0354) < F₂ (1.7655/−0.0634) < O₂ (2.2099/−0.0838) ≈ N₂ (2.2167/−0.0800) < Cl₂ (6.1443/−0.2215). The mild discrepancy in the observed trend in the case of O₂ (slightly smaller polarizability than N₂ but slightly larger magnitude of field stabilization) might be related to the different nature of its ground electronic state (open shell, triplet).

A numerical fit through the origin of |$\Delta E = a{\hbox{\sf\itshape E}\,}^2$|$ΔE=aE2$⁠, where a is the constant of proportionality, yields, when ΔE is expressed in eV and |${\hbox{\sf\itshape E}}$|$E$ in multiples of 10⁹ V m⁻¹, a = −2.0948 for the molecule exhibiting the most pronounced stabilization in response to the field (Cl₂). For the least pronounced case (H₂), and in the same units, a = −0.3345. The corresponding values of a for the remaining molecules are: F₂ (−0.5997), N₂ (−0.7574), and O₂ (−0.7931). Each one of the five regression correlations yield r² = 1.0000 which indicates a perfect fit of the energetic stabilization to a quadratic dependence on the field strength.

The lower plot of Figure 3 displays the field-induced energy change of the hetero-diatomic molecules. The plot can be divided into three behavioral modes depending on the relative magnitude of the permanent dipole moment and polarizability. The first response type is that of HF, which is practically linear with the field strength over the entire range and which can be fit to a linear regression equation |$\Delta E = b{\hbox{\sf\itshape E}}$|$ΔE=bE$⁠, where b = 0.0383 (in multiples of 1.602 × 10⁻²⁸ C.m = 48.02 debye, when ΔE is expressed in eV and |${\hbox{\sf\itshape E}}$|$E$ in multiples of 10⁹ V m⁻¹) is the constant of proportionality with r² = 0.9953 indicating a perfectly linear correlation. This linear dependence of the energy of HF on the field strength results from the combination of a high dipole moment [μ(expt./calc.) = 1.8262/1.8358 debye] and low polarizability [α_// = 0.8658 ų] (Table II), i.e., the field-response in this case is dominated by the linear (dipolar) term in Eq. (4).

The second behavioral mode of the heterodiatomics is that exhibited by CO and NO, two molecules with feeble permanent dipole moments [μ_CO(expt./calc.) = 0.1098/0.0846 debye and μ_NO(expt./calc.) = 0.1587/0.1247 debye] but significant polarizabilities [α_//,NO = 2.2085 ų and α_//,CO = 2.2911 ų]. As a consequence, these two molecules exhibit trends not too dissimilar to those of the homo-nuclear diatomics in that they are generally stabilized by both field directions. At the weaker studied field strengths however, CO and NO are first slightly destabilized in antiparallel fields due to their weak dipole moment until the field strength reaches approximately 2.57 × 10⁹ and 3.60 × 10⁹ V m⁻¹, respectively, when the total energy is lowered again [Figure 3 (bottom)] as the external field becomes sufficiently strong to reverse the direction of the molecular dipole moment into alignment inducing a parallel dipole moment (Figure 2).

A regression of the field stabilization energies of CO and NO against the squared electric field, |$\Delta E = a{\hbox{\sf\itshape E}\,}^2$|$ΔE=aE2$⁠, yields for CO r² = 0.9258 and a = −8.158 × 10⁻⁴ (when ΔE is in eV, and |${\hbox{\sf\itshape E}}$|$E$ is in multiples of 10⁹ V m⁻¹); while for NO r² = 0.8308 and a = −7.703 × 10⁻⁴ (in the same units). The strength of the linear correlation as judged from r² is stronger for CO, the molecule with the smaller permanent dipole moment of 0.0846 debye, and hence a smaller first term in Eq. (4), approaching more closely a homo-nuclear diatomic than for NO which has a stronger dipole moment of 0.1247 debye (Table II).

Finally, and between these two extremes one finds the curve representing the field-response of HCl. This molecule is considerably polarizable (α_// = 2.5958 ų) and simultaneously has a strong permanent dipole moment (⁠|$\bm\mu$|$μ$ = 1.0949 debye) (Table II). Therefore, for HCl the two terms in Eq. (4) have important contributions. The field-induced energy change of HCl shown in Figure 3 exhibits a destabilization in antiparallel fields (first quadrant) that decrease in slope with the magnitude of the field and a rising stabilization in parallel field (third quadrant). These changes can be understood by examining the rate of change of the stabilization energy with the field. Taking the derivative of the expression of ΔE [Eq. (4)] with respect to the field strength:

In the third quadrant (parallel fields for HCl) both terms in Eq. (5) have the same (negative) sign and hence the slope is not constant and decreases with the field strength (as can be seen in Figure 3). At the origin (when E = 0), the tangent to the curve equals the field-free dipole μ₀. Upon moving to the first quadrant (antiparallel fields) the dipole moment enters in the expression for the slope with a positive sign and thus at small magnitudes of the electric field the slope is positive (the curve is rising) and steadily decreasing with |E| (tapering curve) until α_//0|${\hbox{\sf\itshape E}}$|$E$ = μ₀ when the slope become zero (the tangent to the curve at this point is a horizontal line). At field strengths beyond this point, the slope changes sign and the curve takes a downturn. Our data points did not reach fields of enough strength to bring the curve to its turning point, but the curve in Figure 3 clearly exhibits the expected trend.

Figure 4 and Table I show the trends in ΔR and in the actual equilibrium bond length under the external fields |$R_{\sf{E}}$|$RE$⁠, respectively, for the nine molecules. Following a convention similar to the one in Eq. (4), ΔR is defined:

where R₀ is the field-free bond length. The plot for homo-nuclear diatomics (top) shows that increasing the external electric field strength always lengthens the bond. The response of hetero-nuclear diatomics (bottom), however, depends on the relative field-molecule orientation.

The field-induced response of the bond length is clearly nonlinear for all studied molecules (Figure 4). On the bases of the correlation coefficients within the studied range of fields strengths (Table S2 of the supplementary material),⁴² the equilibrium bond length can be fitted closely to an exponential dependence on |${\hbox{\sf\itshape E}}$|$E$⁠:

where γ is a constant of dimensions reciprocal of |${\hbox{\sf\itshape E}}$|$E$⁠. In addition to the exponential regression model (Model I), Table S2 also lists regression equations according to a linear model (Model II). From the listed values of the correlation coefficients in the table, 0.862 (O₂) ≤ r² ≤ 0.881 (H₂), one can see that the exponential relation captures reasonably well the bond stretching under the field for homo-diatomics, with r² values always higher than the linear model for these molecules. The same conclusion applies for the four hetero-diatomic molecules.

For HCl and HF, all parallel fields increasingly stretch the bond with increasing field strength, antiparallel fields compress the bond but to a lesser extent. From Table I, parallel |${\hbox{\sf\itshape E}}_ -$|$E−$ stretches HF and HCl by 0.0066 and 0.0092 Å, respectively. On the other hand, the antiparallel |${\hbox{\sf\itshape E}}_ +$|$E+$ compresses HF and HCl by only 0.0045 and 0.0034 Å, respectively.

The bond lengths of the two heterodiatomics with very small permanent dipole moment (CO and NO) behave in a peculiar manner that differ from the other heterodiatomics that possess large permanent dipole moments. Surprisingly, parallel (rather than antiparallel) fields compress these two molecules quite significantly, by 0.0051 and 0.0032 Å for CO and NO. Antiparallel fields stretch the bond to a lesser (rather that to a larger) extent than their compression by the parallel field, the stretching being, respectively, 0.0071 and 0.0052 Å.

Delley has developed a simple but important model to account for the observed field effects on the bond lengths and on the harmonic frequencies of molecules subjected to external electric fields.⁶³ This author proposed V(r) = D_e(1 − ζ)² + Fr as a generalized Morse potential⁶⁴ that includes a field-perturbation term (Fr) where F = F(E) has the dimensions of force. The form of the equations that followed from this proposal⁶³ are correct and hence achieve excellent statistical fits to direct numerical results.⁶³ Delley's perturbed Morse potential goes to ±infinity [V(∞) = ±∞] for Fr > 0 or Fr < 0, respectively, rendering the interpretation of bond dissociation energy under the field impossible.

Following a similar path, but with modification, we propose here a perturbed Morse potential whereby Fr is replaced by the explicit expression for the field-induced energy change [Eq. (4)]. Further, and in addition to the parallel orientation (the only orientation allowed in Delley's paper), we also consider the antiparallel orientation. Thus, we propose the following field-perturbed Morse potential:

where

and where r is the internuclear separation, |$a = \sqrt {{k}/2D_e} $|$a=k/2De$ is a constant of dimensions of [length]⁻¹ related to the curvature (or force constant k) of the potential energy curve at the field-free equilibrium distance R₀ and which can be taken as the reciprocal of the width of the potential well, and D_e is the (bottom of well) dissociation energy (a positive quantity for a bound state), and where parallel orientations are expressed with a negative sign in front of the dipolar term and antiparallel orientations with a positive sign. The proposed potential [Eq. (8)] has a number of desirable features: (a) It incorporates the explicit expression for the energy change in the field; (b) it takes into account the two relative field-dipole moment orientations (via the ∓ sign) and which can be generalized straightforwardly to any oblique field orientation with respect to the internuclear axis by retaining the terms containing cosθ and sinθ in Eq. (2); and (c) it has a constant and finite asymptotic limit at which μ₀(r → ∞) = 0 and α₀(r → ∞) = constant, the polarizability of the system at r → ∞ being the sum of the polarizabilities of the separate atoms in their respective ground states.

Rewriting Eq. (9) to express the equilibrium bond length under the external field yields

where ζ is obtained by applying the equilibrium condition, ∂V/∂r = 0, to Eq. (8), which gives

where

Only the positive root of the quadratic equation |$\zeta ^2 - \zeta\break - \frac{1}{4} ( { \mp c_1 {\hbox{\sf\itshape E}} - c_2 {\hbox{\sf\itshape E}\,}^2 } ) = 0$|$ζ2−ζ−14(∓c1E−c2E2)=0$ is retained since otherwise ζ would vanish in the field-free case, which is clearly unphysical. It is important to note that for homo-nuclear diatomics, c₁ vanishes since ∂μ₀/∂r = 0, the reason for their IR inactivity. Thus for homo-nuclear diatomics Eq. (11) has the simpler form |$\zeta = \frac{1}{2}( {1 + \sqrt {1 - c_2 {\hbox{\sf\itshape E}\,}^2 } } )$|$ζ=12(1+1−c2E2)$ and the response of the bond length to the field is exclusively dependent on the slope of the polarizability function, upon which the intensity of the Raman signal depends. In the double harmonic approximation, the IR intensity is proportional to the square of the rate of change of the dipole moment with respect to the normal coordinate of the given fundamental⁶⁵ and the Raman intensity to the square of the derivative of the polarizability and that of its anisotropy.⁶⁶ 

Since we did not explicitly calculate the derivatives in Eqs. (12) in this preliminary work, and in order to test the validity of the mathematical form of Eqs. (10)–(12), we resort to a numerical fitting of the bond lengths calculated directly from the Gaussian 09 program (brute force calculation) to Eq. (10) (Table III). For homo-nuclear diatomics, the fitting of c₂ includes, in addition to the field-free value, six values of the bond length up to a maximum field intensity of 1.03 × 10¹⁰ V m⁻¹, i.e., 7 data points, with a data-to-fitting parameters ratio of 7:1. In the case of the hetero-nuclear diatomics, both c₁ and c₂ are non-zero and the external fields range from –1.03 × 10¹⁰ V m⁻¹ to +1.03 × 10¹⁰ V m⁻¹ (13 data points) yielding a data/parameters ratio of 13:2.

The values of the squared correlation coefficients and those of the root mean square deviations (RMSD) in Table III demonstrate that the fitted and directly calculated values of the bond length as a function of E are in perfect agreement. Table IV lists the directly calculated equilibrium bond lengths under different field strengths and direction for the set of the nine studied molecules in comparison with those obtained from the fitted parameters and Eq. (10). The table shows that in addition to a universal agreement of trends (shortening or lengthening of bonds), the quantitative agreement between the fitted and the brute force values is remarkable in absolute and relative (%) terms, differences being only discernable in the fifth decimal.

External fields stretch all studied homo-nuclear diatomic molecules but can stretch or compress hetero-nuclear diatomics depending on the orientation of the molecule in the field. Generally, stretching weakens a chemical bond since a stretching to infinity is tantamount to bond breaking. A measure of bond strength⁶⁷ (or more precisely, of its stiffness) is the force constant, k, which depends on the atomic numbers of the bonded atoms and not on their isotopic composition in contrast to the corresponding vibrational frequency, ν, which depends on the atomic masses as well.

The field-induced change in the force constant and in the harmonic frequency is defined, respectively:

and

where the subscript 0 refers to the field free value and E to the value under the perturbing field.

Figure 5 displays plots of Δk and Δν for all molecules against the strength and direction of the field. As the figure and Table I demonstrate, external fields always reduce k for all of the five studied homo-nuclear diatomic molecules, an effect which is maximal in the case of O₂ (k₀ = 26.26 mdynes Å⁻¹, |${k}_{{\sfb{\scriptsize E}}}$|$kE$ = 25.62 mdynes Å⁻¹). The reduction in the force constant of homo-nuclear diatomics is not minimal for N₂, a molecule bound by one of the strongest chemical bonds known [k₀(N≡N) = 47.70 mdyne Å⁻¹]. Instead, bonds that are only 1/4th as stiff (H–H and F–F) exhibit the smallest field-induced softening: In mdynes Å⁻¹, the force constants of these two bonds drop from 11.52 to 11.35 for H₂ and from 10.94 to 10.78 for F₂. In other words, the magnitude of the field-induced softening of homo-nuclear chemical bonds bears no simple relation to its stiffness in field-free conditions.

The ordering of the magnitude of the Stark-shift is different than that of the force constants due to its dependence on the reduced mass (M_red.). Thus, N₂ and F₂ show the smallest Stark-shift, while H₂ – the molecule with the smallest reduced mass – the largest (Figure 5 and Table I). A glance at Figure 5 shows that, generally, both Δk and Δν are not linear functions of the field strength, the nonlinearity being generally more striking for homo-nuclear diatomics.

In the case of the two halogen halides, parallel fields decrease k (soften the bond) the same relative fields orientation that induce the opposite effect (stiffen the bond) in the case of CO and NO (Table I). The same observation can be reached by noting the rising positive values of Δk in the first quadrant of the lower right plot of Figure 5. HCl exhibits the least pronounced change in response to the field, in mdyne Å⁻¹: k₀ = 5.57, |${k}_{{{\sfb{\scriptsize E}}} + }$|$kE+$ (antiparallel, stiffer) = 5.65; |${k}_{{{\sfb{\scriptsize E}}} - }$|$kE−$ (parallel, softer) = 5.27. In contrast, NO (a bond with considerable stiffness) exhibits the largest change, in mdyne Å⁻¹ (CCSD level of theory): k₀ = 34.91, |${k}_{{{\sfb{\scriptsize E}}} + }$|$kE+$ (parallel, stiffer) = 37.01; |${k}_{{\scriptsize{\sfb{E}}} - }$|$kE−$ (antiparallel, softer) = 33.29.

Interestingly, the (U)QCISD/6-311++G(3df,2pd) level of theory, the principal level of theory utilized in this paper, fails to reproduce physically meaningful trends of the Δk and Δν field-responses of the NO molecule (supplementary material).⁴² 

The field-response of the force constant and of the harmonic frequency of this molecule calculated at the UQCISD/6-311++G(3df,2pd) level of theory, each, lies on a parabolic curve. This parabolic relationship is in disagreement with the linear response of k to the field observed for the remaining eight molecules including CO (which is similar to NO in its dipole moment, polarizability, field-free force constant, and charge separation). The parabolic curve is also inconsistent with the linear response of k and ν of NO obtained in response to the perturbing fields at the B3LYP/6-311++G(3df,2pd) level of theory.

In order to resolve the inconsistency, we have tested the field-effect on the force constant of this molecule at four additional commonly used computational levels of theory using the same basis (Figure 6). The details of these computations are outlined in Sec. II B.

A glance at Figure 6 shows that the various levels of theory yield very different field-free force constants for NO, these values in mdyne Å⁻¹ are: H–F (43.33), MP2 (88.66), B3LYP (34.32), mPW1PW91 (36.17), CCSD (34.91), and QCISD (29.09). The consensus value appears to be around the CCSD value, i.e., ∼35 mdyne Å⁻¹. MP2 appears to dramatically overestimate this value while QCISD underestimates it.

The experimental (field-free) frequency for this molecule is reported to be 1904 cm⁻¹,⁶⁸ while the theoretical methods, in order of increasing absolute deviation from the experimental value, are [Method (ν_calc., deviation from experimental value; all in cm⁻¹)]: B3LYP (1979, +75); QCISD (1822, –82); CCSD (1996, +92); mPW1PW91 (2032, +128); HF (2224, +320); MP2 (3181, +1277). The unscaled B3LYP frequency of NO is thus the closest in agreement with experiment, outperforming even CCSD in this case.

The Stark shift is consistent at all levels of theory (except QCISD) and is reflected in the slopes of the lines representing k as a function of E displayed in Figure 6. The slopes of these curves (in 10⁻⁹ mdyne Å⁻¹ V⁻¹ m) and the corresponding squared correlation coefficients (r²) of the regression lines, respectively, obtained at the various levels of theory are: H–F (0.178/0.999), MP2 (3.323/0.996), B3LYP (0.121/0.995), mPW1PW91 (0.126/0.995), CCSD (0.174/0.997). Thus, five different levels of theory, including the highly accurate CCSD, predict a linear response of the force constant (and hence of the frequency) of NO to the external field strength, just as any of the other eight molecules. It is concluded that QCISD while reproducing well the field-free frequency is incapable of correctly reproducing the associated Stark-shift, at least with the used basis set. Hence, the values of k and ν listed for NO in all tables and figures in this paper are obtained from UCCSD/6-311++G(3df,2pd) calculations.

Following similar step as those of Delley⁶³ but starting from the field-perturbed potential [Eq. (8)], the model is elaborated further to rationalize the observed trends in the field-response of k and of ν.

All the bonds we are investigating in this study are of the strong covalent type, accurately described by the harmonic approximation near the bottom of the well, whereby

where M_red. is the reduced mass.

The second derivative of the field-free Morse potential yields

which when evaluated at r = R₀ gives

which is the field free force constant. But since the field alters the equilibrium bond length, i.e., |$R_0 {\mathop{\hbox to 2pc{\rightarrowfill}}\limits^{\hbox{\sf\itshape E}}}$|$R0→E$|${R}_{\scriptsize{\sfb{E}}}$|$RE$⁠, we keep ζ and solve for |${k}_{{\sfb{\scriptsize E}}}$|$kE$ to obtain

Converting the harmonic force constant into the frequencies using the last equality in Eq. (15):

where

From Eqs. (19) and (20) it is clear that the field-perturbed frequency |$\nu _{\scriptsize{\sfb{E}}}$|$νE$ is given by the field-free value ν₀ multiplied by the dimensionless factor κ which depends on both the magnitude and direction of the external field. When κ(E) > 1 this corresponds to a hypsochromic shift, while a bathochromic shift corresponds to κ(E) < 1. Obviously κ(E = 0) ≡ κ₀ = 1 yields the field-free frequency. This formula should therefore explain the trends in Figure 5, whether the Stark shift is to higher or lower frequencies.

Thus, the necessary quantities to evaluate the frequency under any given field are the field-free frequency ν₀ weighted by κ. The latter, κ, can be evaluated, at least in principle, from field-free IR⁶⁵ and Raman⁶⁶ through its dependence on c₁ and c₂ defined in Eq. (12). For homo-nuclear diatomics, c₁ = 0, and hence Eq. (20) reduces to

Regression estimates of c₁ and c₂ were obtained from the nonlinear fitting of Eq. (10) (Table III), as explained above. Such nonlinear regression with more than one fitting parameter can have multiple local minima for the function being optimized with no guarantees that the numerical minimization of residuals will locate the “global minimum,” i.e., a particular solution depends on the initial guess.

Table V lists the field-free and Stark-shifted frequencies calculated directly from the Gaussian 09 program against those estimated using the fitted values of the c₁ and c₂ parameters (listed in Table III) by substitution into Eqs. (19)–(21). The table also lists the signed percent errors in the estimated frequencies. The table reveals that the correct trends of bathochromic or hypsochromic Stark shifts as a function of the field strength and direction are well reproduced. Further, it reveals that the absolute errors are insignificant especially at relatively low-to-intermediate field strengths, and that larger absolute and relative errors are generally encountered only at the strongest field intensities. The agreement is remarkable given the approximations of the model and the statistical “local minimum” nature of the fitted solution.

A number of molecular and bond properties of a set of five homo-nuclear and four hetero-nuclear diatomic molecules are studied as functions of external homogenous static electric fields. The homo-nuclear diatomics include H₂, N₂, O₂, F₂, and Cl₂ and the hetero-nuclear diatomics include HF, HCl, CO, and NO. The selected molecules cover ranges and combinations of permanent dipole moments and polarizabilities: (1) The homo-nuclear diatomics, lacking a permanent dipole moment, range in polarizability from H₂ (0.8 ų) to Cl₂ (4.6 ų); (2) the hetero-nuclear diatomics include those with a large permanent dipole moment and a small polarizability such as HF (μ = 1.8 debye, α = 0.8 ų), those with a negligible dipole moment but sizable polarizability such as CO (μ = 0.1 debye, α = 2.0 ų), and those with a large dipole moment and a large polarizability such as HCl (μ = 1.1 debye, α = 2.6 ų). The response of the total energy to the field ΔE is shown to depend on the relative sizes of the dipolar and polarizability terms in the energy expression of each individual molecule, consistent with the explicitly calculated trends.

Delley proposed a Morse potential augmented with a term to account for the external field-effect.⁶³ The field term proposed by Delley leads to infinities in energy at infinite internuclear separations. We propose to replace that term with the explicit expression for the energy change of the molecule in an external field and also to allow for both the parallel and antiparallel orientations. Following Delley's derivations but starting from our field-perturbed Morse potential we have refined a simple predictive model that correctly accounts for the responses of the energies, bond lengths, vibrational frequencies, and force constants of diatomics to externally applied homogenous electric fields.

A statistical fitting of one parameter, c₂, proportional to the rate of change of the polarizability with internuclear separation (hence related to Raman intensities), is shown to capture the essential physics of the bond length elongation in an external field and the associated bathochromic Stark shift of the homo-nuclear diatomic molecules. In the case of the hetero-nuclear diatomics, which are both Raman- and IR-active, an additional parameter proportional to the rate of change of the dipole moment with the bond length (hence related to IR intensities), c₁, is necessary to explain both the bond length changes in the fields and the concomitant bathochromic or hypsochromic Stark shifts, depending on the orientation of the molecule in the field.

The model presented here demonstrates the predictability of the field-responses from the field-free molecular properties with only a parametric dependence on the external electric field. It is also emphasized that the prediction of equilibrium bond lengths under perturbing fields and the accompanying Stark-shift require knowledge of quantities that are all measurable, at least in principle, from field-free experiments. These measurable quantities include the bond dissociation energy, the Raman- and IR-intensities, the equilibrium bond length, the vibrational frequency, and the curvature of the potential energy surface at the equilibrium bond length.

This paper reports some of the results of an M.Sc. thesis submitted by S.S.-H. to Saint Mary's University in partial fulfillment of the requirements for the Master's Degree in Applied Science.

The authors are indebted to Professor André D. Bandrauk (Université de Sherbrooke), Professor Bernard Delley (Paul Scherrer Institut), Professor Lou Massa (Hunter College, City University of New York), and Professor Francesc Illas (Universitat de Barcelona & Parc Científic de Barcelona) for illuminating discussions. The financial support of this work was provided by the Natural Sciences and Engineering Research Council of Canada (NSERC), Canada Foundation for Innovation (CFI), Mount Saint Vincent University, and Saint Mary's University.