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

It is shown that the response of molecular properties of diatomics such as the total energy, the bond length, and the vibrational Stark shift to an external homogenous electric field (EF) can be predicted from field-free observable properties such as the equilibrium bond length, the bond dissociation energy, the polarizability and dipole moment functions, and the vibrational frequency. Delley [J. Mol. Struct.: THEOCHEM 434, 229 (1998)] suggested to approximate the potential energy surface under an EF by a Morse function augmented with a EF term proportional to the internuclear separation. In this work, this term is replaced by the expression of the field-induced energy change which yields a field-perturbed Morse potential that tends to a constant asymptotic limit when the EF term itself become proportional to the sum of the polarizabilities of the separated atoms. The model is validated by comparison with direct calculations on nine diatomics, five homo-nuclear (H2, N2, O2, F2, and Cl2) and four hetero-nuclear (HF, HCl, CO, and NO), covering a range and combinations of dipole moments and polarizabilities. Calculations were conducted at the quadratic configuration interaction with single and double excitations (QCISD) and density functional theory (DFT)-B3LYP levels of theory using the 6-311++G(3df,2pd) basis set. All results agree closely at the two levels of theory except for the Stark effect of NO which is not correctly predicted by QCISD calculations as further calculations, including at the coupled cluster with single and double excitation (CCSD) level of theory, demonstrate.


I. INTRODUCTION
][3][4][5][6][7][8][9][10][11][12][13][14][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 9 V m −1 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." 16Examples of "local reporter" molecules include carbon monoxide (C=O) attached to the Fe of the porphyrin ring in myoglobin, 17 and the nitrile-(-C≡N) containing substrates in the active site of human aldose reductase a) E-mail: cherif.matta@msvu.ca.Telephone: +1(902)-457-6142.enzyme (hALR2). 16From the changes in the vibrational Stark shifts of the nitrile group, Boxer et al. were able to estimate a change of 10 8 V m −1 in the electric field accompanying a single point mutation at the active site, 16 but interference of the effect of hydrogen bonding with the Stark effect in enzyme active sites can complicate the interpretation of spectra. 18,19 ang et al. 20 and Gascón et al. 21have 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 9 V m −1 but can reach ∼10 10 V m −1 . 22Fields of strengths ∼10 9 V m −1 have been shown to quadruple the rate constant of double proton transfer reactions, 15 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, 23 and are found between the tip and the sample in a scanning tunneling microscope (STM) under normal operating conditions. 24Fields of 10 8 V m −1 accelerate photosynthetic reactions by an order of magnitude, [25][26][27] and significantly decrease the enzymatic activity of cytochrome c oxidase. 28he maximal strength of macroscopic fields applied to a crystal during an X-ray crystallographic diffraction is around 10 7 V m −1 . 29The field-induced changes in diffraction patterns have been observed for sometime, [30][31][32][33][34][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. 29n 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.

B. Electronic structure calculations
Only the ground electronic states of the diatomics were considered: 1 + g for closed-shell homo-nuclear dimers, 1 + for closed-shell hetero-nuclear dimers, 2 r in the case of NO, and 3 − g in the case of O 2 .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), 37 and (ii) density functional theory (DFT) Becke's hybrid exchange functional 38 and Lee-Yang-Parr correlation functional 39 (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 −6 and 1 × 10 −6 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. 40All reported frequencies are those obtained directly from the calculation (unscaled). 41oth 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. 42here 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][44][45] with single and double excitations (CCSD), Møller-Plesset second order perturbation theory 46 with full excitation including core electrons, Hartree-Fock (HF) self consistent field (SCF) level, 47 and the DFT functional mPW1PW91. 48

C. Field strengths and directions
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πε 0 a 2 0 ≈ 5.14 × 10 11 V m −1 .In this work, electric fields from 5.14 × 10 8 to 5.14 × 10 10 V m −1 (=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. 49Any 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 |E| = E is the electric field strength, IP is the vertical ionization potential, both quantities in a.u., and ω 0 is the atomic unit of frequency (4 × 10 16 s −1 ).Tunneling ionization rates calculated from Eq. (1) using the experimental IPs of the nine diatomics are collected in Table S1 of the supplementary material. 42hemical reactions generally happen on a time scale of ∼10-10 4 fs. 52One can take the reciprocal of 10 −4 fs −1 as the maximal acceptable average tunneling ionization rate for a molecule in the field.An arbitrary order of magnitude threshold of around ω ∼ 10 11 s −1 can thus be considered reasonable.Given these considerations, the values listed in Table S1 suggest a field strength of 5.14 × 10 9 V m −1 (0.01 a.u.) as an upper ionization threshold with the exception of alkali metal dimers Li 2 and Na 2 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) 53 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, E + or − → E + mean a field of magnitude 1.03 × 10 10 V m −1 = 2.0 × 10 −2 a.u.oriented to point in the positive direction of the z-axis shown in Figure 1 while E − or ← − 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 +← − μ − 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 parallelplate 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 54 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 destabi-lizing (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 (fieldfree) 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.75 H ← F −0.75 or +0.26 H ← 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.22 C → O −1.22 and +0.44 N → 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][58][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. 60

III. RESULTS AND DISCUSSION
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 10 V m −1 ), 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 (μ), 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 homonuclear diatomics the two field directions are equivalent, the fields are assigned a positive sign.

A. Energy and dipole moment of a diatomic molecule in external homogenous electric fields 1. The energy expression
The expansion of the energy [E(E, r) ≡ E E ] of a diatomic molecule in an external homogenous electric field (E) as a power series takes the form 61   where E 0 , μ 0 , α //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; 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, μ 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 μ 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 μ 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. 62rther, 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.

Dipole moment and polarizability
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 2 for which α ≈ 0.8 Å 3 and up to the highly polarizable Cl 2 with α ≈ 4.6 Å 3 .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 2 ; CO and NO have small permanent dipole moments (∼0.1 debye) and sizeable polarizabilities (∼1.8 Å 3 ); and finally HCl has both a sizable dipole moment (∼1.1 debye) and a sizable polarizability (∼2.5 Å 3 ).The listing in the table also V m −1 .The convention of assigning directions to the field is given by the arrows parallel to the abscissa of the plot (bottom) in which the field changes direction halfway through the abscissa.In the inset of the bottom plot, each molecule is drawn in the orientation it has with respect to the external fields with a small arrow between the atomic symbols depicting the orientation of the permanent molecular dipole moment with respect to the external field (also see Figure 1).Except when stated otherwise, all plotted results were obtained at the (U)QCISD/6-311++G(3df,2pd) level of theory.] 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 2 , the one with the smallest polarizability among the homo-nuclear diatomics, while the molecule with the largest polarizability, Cl 2 , 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 10 V m −1 ) a sorting of the homonuclear molecules X 2 in order of increasing polarizability parallels a corresponding sorting in the magnitude of the stabilization energy.Thus, (α // in Å 3 (2.2167/−0.0800)< Cl 2 (6.1443/−0.2215).The mild discrepancy in the observed trend in the case of O 2 (slightly smaller polarizability than N 2 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 E = aE 2 , where a is the constant of proportionality, yields, when E is expressed in eV and E in multiples of 10 9 V m −1 , The signs of the permanent dipole moments indicate their relative orientation with respect to the coordinate system displayed in Figure 1.d The arrows between the atomic symbols indicate the direction of the dipole moment according to the physicist convention.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 E = bE, where b = 0.0383 (in multiples of 1.602 × 10 −28 C.m = 48.02debye, when E is expressed in eV and E in multiples of 10 9 V m −1 ) is the constant of proportionality with r 2 = 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.8358debye] and low polar-izability [α // = 0.8658 Å 3 ] (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.0846debye and μ NO (expt./calc.)= 0.1587/0.1247debye] but significant polarizabilities [α //,NO = 2.2085 Å 3 and α //,CO = 2.2911 Å 3 ].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 9 and 3.60 × 10 9 V m −1 , 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, E = aE 2 , yields for CO r 2 = 0.9258 and a = −8.158× 10 −4 (when E is in eV, and E is in multiples of 10 9 V m −1 ); while for NO r 2 = 0.8308 and a = −7.703× 10 −4 (in the same units).The strength of the linear correlation as judged from r 2 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 Å 3 ) and simultaneously has a strong permanent dipole moment (μ = 1.0949 debye) (Table II).Therefore, for HCl the two terms in Eq. ( 4) have important contributions.The fieldinduced 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 μ 0 .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 E = μ 0 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.

Description of trends
Figure 4 and Table I show the trends in R and in the actual equilibrium bond length under the external fields R E , respectively, for the nine molecules.Following a convention similar to the one in Eq. ( 4), R is defined: where R 0 is the field-free bond length.The plot for homonuclear 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), 42 the equilibrium bond length can be fitted closely to an expo-FIG.4. Plots of the change in the bond lengths ( R), at the optimized geometry, as a function of the electric field (E) strength for the homo-nuclear diatomics (top), and as a function of the field strength and direction for the hetero-nuclear diatiomics (bottom); all in mÅ.The change in the BL is defined: R = R E -R 0 , where R E is the bond length in the field and R 0 the bond length in absence of external fields (bond lengths of the ground state molecules are listed in Table I) (see end statement of the caption of Fig. 2 in square brackets).nential dependence on E: where γ is a constant of dimensions reciprocal of 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 2 ) ≤ r 2 ≤ 0.881 (H 2 ), one can see that the exponential relation captures reasonably well the bond stretching under the field for homodiatomics, with r 2 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 E − stretches HF and HCl by 0.0066 and 0.0092 Å, respectively.On the other hand, the antiparallel 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 Å.

Field-perturbed Morse potential
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. 63This author proposed V (r) = D e (1 − ζ ) 2 + F r as a generalized Morse potential 64 that includes a fieldperturbation term (Fr) where F = F(E) has the dimensions of force.The form of the equations that followed from this proposal 63 are correct and hence achieve excellent statistical fits to direct numerical results. 63Delley'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 field-free equilibrium distance R 0 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 μ 0 (r → ∞) = 0 and α 0 (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 = 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 1 vanishes since ∂μ 0 /∂r = 0, the reason for their IR inactivity.Thus for homo-nuclear diatomics Eq. ( 11) has the simpler ) 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 65 and the Raman intensity to the square of the derivative of the polarizability and that of its anisotropy. 66ince 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 2 includes, in addition to the field-free value, six values of the bond length up to a maximum field intensity of 1.03 × 10 10 V m −1 , 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 1 and c 2 are non-zero and the external fields range from  -1.03 × 10 10 V m −1 to +1.03 × 10 10 V m −1 (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.

C. Vibrational Stark shift 1. General considerations
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 67 (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.The direction of the field with respect to the permanent molecular dipole moment is indicated by the arrow, while the permanent dipole moment is indicated by the arrow between the atomic symbols.c For every field there are three rows: The top row is the value calculated directly from the Gaussian 09 program, the middle row is the value obtained from using the regression parameters c 1 and c 2 listed in Table III and obtained from a fitting of bond lengths using Eq. ( 10), and the third row is the signed % error of the fitted value.
effect which is maximal in the case of O 2 (k 0 = 26.26mdynes Å −1 , k E = 25.62 mdynes Å −1 ).The reduction in the force constant of homo-nuclear diatomics is not minimal for N 2 , a molecule bound by one of the strongest chemical bonds known [k 0 (N≡N) = 47.70 mdyne Å −1 ].Instead, bonds that are only 1/4th as stiff (H-H and F-F) exhibit the smallest field-induced softening: In mdynes Å −1 , the force constants of these two bonds drop from 11.52 to 11.35 for H 2 and from 10.94 to 10.78 for F 2 .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 2 and F 2 show the smallest Stark-shift, while H 2 -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 Å −1 : This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 140.184.72.92 On: Wed, 15 Apr 2015 12:57:51 FIG. 5. Plots of the change in the force constants in mdynes Å −1 ( k, left), and of the change in the harmonic frequencies in cm −1 ( ν, right plots) as functions of the electric field (E) strength for the homo-nuclear diatomics (top), and as functions of the field strength and direction for the hetero-nuclear diatomics (bottom).The change is defined: P = P E -P 0 , where P E and P 0 are either ν or k with and without E, respectively (the values of ν and k can be found in Table I).All data are obtained at the (U)QCISD/6-311++G(3df,2pd) level of theory except those for NO which were obtained at the UCCSD/6-311++G(3df,2pd) level for reasons described in Sec.III C 2. of the text (see end statement of the caption of Figure 2 in square brackets and also Figure 6).

Inability of UQCISD to reproduce the vibrational Stark-shift of nitric oxide
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). 42he 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.
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  and the corresponding squared correlation coefficients (r 2 ) 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 fieldfree 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.

A simple model accounting for the vibrational Stark-shift in diatomics
Following similar step as those of Delley 63 but starting from the field-perturbed potential [Eq. ( 8)], the model is elaborated further to rationalize the observed trends in the fieldresponse 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 0 gives which is the field free force constant.But since the field alters the equilibrium bond length, i.e., R 0 E − −− →R E , we keep ζ and solve for k E 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 ν E is given by the field-free value ν 0 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) ≡ κ 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 ν 0 weighted by κ.The latter, κ, can be evaluated, at least in principle, from field-free IR 65 and Raman 66 through its dependence on c 1 and c 2 defined in Eq. (12).For homo-nuclear diatomics, c 1 = 0, and hence Eq. ( 20) reduces to Regression estimates of c 1 and c 2 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.The direction of the field with respect to the permanent molecular dipole moment is indicated by the arrow, while the permanent dipole moment is indicated by the arrow between the atomic symbols.c For every field there are three rows: The top row is the value calculated directly from the Gaussian 09 program, the middle row is the value obtained from using the regression parameters c 1 and c 2 listed in Table III [and obtained from a fitting of bond lengths, Eq. ( 10)] using Eq. ( 19), and the third row is the signed % error of the fitted value.d The c 1 and c 2 parameters for NO were obtained as all others from a fitting of bond lengths at the UQCISD/6-311++G(3df,2pd), but the directly calculated vibrational frequencies were obtained from UCCSD/6-311++G(3df,2pd) as explained in the text.
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 1 and c 2 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.

IV. CONCLUSIONS
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 2 , N 2 , O 2 , F 2 , and Cl 2 and the hetero-nuclear diatomics include This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 140.184.72.92 On: Wed, 15 Apr 2015 12:57:51 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 2 (0.8 Å 3 ) to Cl 2 (4.6 Å 3 ); (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 Å 3 ), those with a negligible dipole moment but sizable polarizability such as CO (μ = 0.1 debye, α = 2.0 Å 3 ), and those with a large dipole moment and a large polarizability such as HCl (μ = 1.1 debye, α = 2.6 Å 3 ).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. 63The 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 2 , 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 IRactive, an additional parameter proportional to the rate of change of the dipole moment with the bond length (hence related to IR intensities), c 1 , 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.

FIG. 1 .
FIG. 1. Orientation of hetero-nuclear molecules along axis in the coordinate system where A and B refer to the less and more electronegative atom, respectively.Atom A is placed at the origin and the internuclear axis aligned along the z-axis.The electronegativity ranking is taken according to Pauling to be 53 H (2.1) < C (2.5) < N (3.0) ≈ Cl (3.0) < O (3.5) < F (4.0).The electric field E is positive when oriented as in this figure, pointing to the positive z-direction, and negative when pointing to the negative z-direction.The arrows between the atomic symbols of the diatomics point at the direction of the permanent (field-free) molecular dipole moment with the value of the dipole moment calculated at the QCISD level to two decimal places in debye (a ± sign indicates a parallel/antiparallel orientation with respect to both the z-axis and the E-field).

144101- 5 S
FIG.2.Plots of the molecular dipole moment (in debyes) as a function of the electric field (E) strength for the homo-nuclear diatomics (top), and as a function of the field strength and direction for the hetero-nuclear diatiomics (bottom).The following statement in square brackets applies to this figure and to Figures 3-5 as well: [The electric field magnitude is given in 10 9 V m −1 .The convention of assigning directions to the field is given by the arrows parallel to the abscissa of the plot (bottom) in which the field changes direction halfway through the abscissa.In the inset of the bottom plot, each molecule is drawn in the orientation it has with respect to the external fields with a small arrow between the atomic symbols depicting the orientation of the permanent molecular dipole moment with respect to the external field (also see Figure1).Except when stated otherwise, all plotted results were obtained at the (U)QCISD/6-311++G(3df,2pd) level of theory.]

FIG. 3 .
FIG.3.The change in the total energy ( E) without ZPE or vibrational corrections as a function of the electric field (E) strength for the homo-nuclear diatomics in milli-electon volts (meV) (top), and as a function of the field strength and direction for the hetero-nuclear diatomics in electron volts (eV) (bottom).The change in the total energy is defined: E = E E -E 0 , where E E is the total energy in a finite non-zero field E and E 0 the energy in absence of external fields (energies of the field-free ground state molecules and their perturbed values under the strongest field considered in this study are listed in TableI) (see end statement of the caption of Fig.2in square brackets).
) and where r is the internuclear separation, a = √ k/2D e is a constant of dimensions of [length] −1 related to the curvature (or force constant k) of the potential energy curve at the This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 140.184.72.92 On: Wed, 15 Apr 2015 12:57:51

FIG. 6 .
FIG.6.Plots of the force constant, k, of the bond in N-O (in mdyne Å −1 ) as a function of the electric field (E) strength and direction (in 10 9 V m −1 ) obtained from six different underlying electronic structure methods (levels of theory) all using the same basis set [6-311++G(3df, 2pd)]: Hartree Fock (H-F), Møller-Plesset second order perturbation theory (MP2), density functional theory using the two functionals B3LYP and mPW1PW91, coupled clusters with single and double excitations (CCSD), and quadratic configuration interaction with single and double excitations (QCISD).The relative field-molecule orientation is the same as in all previous figures.

TABLE I .
Molecular, bond, and atomic properties of diatomics with and without an electric field (E ± ) of ±1.03 × 10 10 17365 − 109.35884− 150.11470 − 199.27451 − 919.39999 − 100.32095 − 460.31704 − 113.14531 − 129.70479 0 − 1.17235 − 109.35590− 150.11162 − 199.27218 − 919.39185 − 100.33424 − 460.32223 − 113.14149 − 129.70083 ← − E − − 100.34986 − 460.33429 − 113.14401 − 129.70285 a Data based on calculations at the (U)QCISD/6-311++G(3df,2pd) level of theory for all properties and molecules [except for the force constants and frequencies of NO which were obtained at the UCCSD/6-311++G(3df,2pd) level of theory].bThearrow between the atomic symbols depicts the direction of the field free (permanent) dipole moment.Note that this direction may flip sides under a strong external field in the opposite direction [see also footnote (c)].c A negative dipole moment points to the left (−μ ≡ ← − μ ) and one that is positive to the right with respect to the other vectors indicated by arrows in this table [see also footnote (b)].

TABLE II .
Experimental and calculated polarizabilities (in cubic angstroms, Å 3 ) and permanent dipole moments (in debye) of the ground states of the diatomic molecules considered in this study.

TABLE III .
Regression parameters, correlation coefficients, and root mean square deviations of bond lengths obtained from a fitting to Eq. (10).a

TABLE IV .
Directly calculated equilibrium internuclear separations (bond lengths) in Å under different fields and those obtained from the regression parameters listed in Table III using Eq.(10).a a All data are based on calculations at the (U)QCISD/6-311++G(3df,2pd) level of theory.b This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 140.184.72.92 On: Wed, 15 Apr 2015 12:57:51

TABLE V .
Directly calculated harmonic vibrational frequencies (in cm −1 ) under different fields and those obtained from the regression parameters listed in TableIIIusing Eq. (19).a a All data except for NO are based on calculations at the (U)QCISD/6-311++G(3df,2pd) level of theory.b