Exploring the relationship between vibrational mode locality and coupling using constrained optimization Free

Vibrational spectroscopy is one of the most powerful means available for characterizing molecular structures and reactions. Techniques like Raman and multi-dimensional IR are capable not only of identifying functional groups and describing their environment but also of capturing dynamical behavior and solvent interactions.^1,2 Assigning these spectra, however, is challenging due to the number and overlapping nature of vibrational peaks within relatively narrow energy regions. Detailed interpretation of vibrational spectra therefore relies on simulation to identify the peaks and derive atomistic meaning.

In order to accurately model vibrational spectra, a number of powerful techniques have been developed. Vibrational self-consistent field (VSCF),^3–8 configuration interaction (VCI),^9,10 coupled cluster (VCC),¹¹ and perturbation theory (VPT)^12–15 are designed to systematically account for mode correlations and provide improved descriptions of vibrational properties. Regardless of which of these methods is chosen, a representation for the potential energy surface (PES) must be constructed.^16,17 This often requires inclusion of the effects of many-mode interactions by high-order polynomials, which is computationally expensive. For all but the simplest systems, it is usually necessary to use a truncated “n-mode” representation, where n is the highest order coupling term included in the model.^18,19 Even then, the poor scaling of even a 2-mode representation means it is desirable to omit mode-couplings wherever possible. The accuracy of these models is dependent on the magnitude of the omitted coupling terms, which is in turn sensitive to the vibrational mode coordinate system^20,21 used to construct the model PES.

Ideally, vibrational modes should be chosen for their ability to capture the underlying physics and the ease with which they can be interpreted. The standard basis vectors for describing nuclear vibrations are the normal modes (NMs), which are the eigenvectors (with nonzero eigenvalues) of the mass-weighted Hessian. Normal modes are frequently generated using Cartesian coordinates, but may also employ internal coordinates²² to better capture torsions, bends, and rotations that are not accurately described by a rectilinear coordinate system.^23–25 Because they are based on second-derivatives of the potential energy, normal modes provide a reasonable zero-order approximation for vibrational motion and are ideal when the potential is nearly or completely harmonic. In anharmonic systems, however, normal modes produce off-diagonal Hessian terms that vary rapidly as the molecule samples the vibrationally accessible region of phase space. Normal modes become an especially poor choice when investigating properties like dissociation²⁶ and Fermi resonance frequency shifts.²⁷ 

As system size increases, normal modes become more difficult to interpret. Where chemists frequently think of molecular and vibrational features in terms of small fragments, normal modes are generally non-local and can involve the motion of many nuclei, even in a sizeable system.²⁸ Due to their near-degeneracy in vibrational frequency, similar functional groups on a large molecule are not described independently by normal modes, even if significant spatial separation produces low degrees of interaction between the fragments.

Several alternatives to normal modes have been proposed that overcome problems associated with intrinsic delocalization and vibrational coupling away from the geometry of the initial Hessian calculation. One approach is to exploit the short-range nature of inter-nuclear interaction. These procedures usually start from the normal modes and perform unitary transformations to increase the locality of the vibrations.^28,29 Another popular approach is to choose a set of modes to localize based on heuristic considerations, as explored by Child,^30,31 Mortensen,³² and others.^29,33,34

Local modes (LMs) provide interesting advantages over normal modes for several reasons. The short-range interactions in a molecule lead to a rapid decay of mode-mode coupling terms as the distance between the modes increases.²⁹ Panek and Jacob³⁵ found this property of local modes to have a definite benefit to the convergence behavior of both VSCF and VCI calculations performed on an alanine hexamer. Also, because local modes represent motion of only a few atoms, interpretation of their motion is much simpler when compared with the highly delocalized normal modes. Local modes tend to itemize the molecular system into recognizable fragments (e.g., of amide, carbonyl, and alcohol groups), which are not only familiar to chemical intuition but also typically correspond to transferrable descriptions of local chemical units. Similar principles have been successfully shown to apply to localized orbitals,^36–39 which have been used extensively and successfully in electronic structure theory.^40–45 

As suggested by Reiher²⁸ and Steele,²⁹ localized vibrations can provide an improved basis for describing weakly coupled, anharmonic oscillators. For instance, X–H stretches (X = N, C, O, …) are known to benefit from localization as the effective anharmonic coupling decreases when the modes are localized.^30,31,46 For example, Figure 1(a) shows a representative stretching normal mode of an HF tetramer. Despite the system being composed of four nearly independent anharmonic oscillators, the vibrational motion is still delocalized. The right side of the figure shows the analogous local mode. Besides the interpretive benefit that comes from single-molecule motion, Yagi et al.²¹ have shown that the magnitude of anharmonic coupling is reduced as well by localizing these stretches.

Not all vibrational modes, however, benefit from total localization. Figure 1(b) shows the result of localizing a normal mode consisting of C–H bends and ring distortion in stilbene. The maximally localized vibration compresses this mode into motion of just two nuclei. The apparent frequency of the local mode (1701 cm⁻¹) deviates significantly from the harmonic frequency (1247 cm⁻¹), and such deviations in the diagonal elements of the Hessian only occur with the introduction of significant coupling to neighboring modes. As a result, the fully localized mode is a worse representation of the underlying physics and is not easy to interpret using chemical intuition.

These results suggest that an optimal degree of localization may exist somewhere between the delocalized normal modes and the fully localized modes. Examples of insight into this area include work by the Sibert group, where localized vibrations involving methylene fragments are used to precisely identify a number of structural properties in organic compounds.^27,33 By employing localized internal coordinate representations of the alkyl and aromatic C–H stretches within a molecule, the effects of Fermi resonances are captured, and rotational-conformers arising from “soft” alkyl torsions can be distinguished.^34,47,48 While these heuristic approaches must be formed on a case-by-case basis, the general principle appears sound and deserves further attention.

In this light, Yagi and co-workers²¹ proposed a general algorithm for optimizing vibrational basis sets. Since the VSCF energy is not invariant to unitary transformation of the coordinates, the coordinate system can be included as variational degrees of freedom (DOF) to minimize the energy of the target vibrational state. The resulting vibrational basis is typically neither completely localized nor completely delocalized. VCI computations using the optimized coordinates require lower-order excitations to achieve comparable accuracy to when normal modes are used as the coordinate system, resulting in faster calculations. Further analysis showed that reduced anharmonic coupling in the new coordinate system is responsible for the improved convergence. Unfortunately, the iterative VSCF energy minimization is a costly procedure,⁴⁹ restricting its applicability to smaller systems.

These considerations indicate that while localization does reduce mode-mode coupling at long distances, not all modes benefit from localization when coupling is considered. Herein, we propose a systematic means of traversing between the fully delocalized normal modes and fully local modes by explicitly limiting the intermode coupling. This differs from recent use of local modes as seen in Steele²⁹ and Reiher²⁸ in that we create bases of intermediate, rather than complete, localization. Such partial localization is achieved by including a constraint⁵⁰ that places a maximum threshold on the coupling in the mass-weighted Hessian during localization. This strategy allows mode pairs that are weakly coupled to localize while preserving the delocalized description of the remaining, more coupled mode pairs. Because this procedure requires only a single Hessian computation, this procedure is readily applicable to sizable polyatomic molecules. The theory motivating the partial localization scheme is discussed below, followed by algorithmic details and examples of its application to several polyatomic molecules.

We seek to determine when vibrational coordinate localization will minimize anharmonic coupling in the Hessian. To do so, we consider a one dimensional system consisting of three masses (Figure 2) as previously discussed by Steele and Cheng.²⁹ We define the potential to be anharmonic and let the masses be m₂ = 10m₁, m₁ = m₃. The total potential as a function of the positions x_i is therefore,

where the parameter ζ represents the anharmonicity in V(x). After mass-weighting, the coordinates $xi=ui/mi$⁠, the Hessian, A, becomes

If V_tot is harmonic (i.e., ζ = 0), or the Hessian is evaluated at the relaxed spring positions (u₁ = u₂ = u₃ = 0) then the eigenvalues of A are 0, 1, and 6/5 corresponding to a translational mode, an asymmetric mode, and a symmetric mode, respectively. We note that the two nonzero eigenvalues are quite close in frequency, indicating that localization may be appropriate for this system.⁵⁰ Using a unitary rotation matrix, U,

the two internal vibrations V can be transformed to create a new basis R, and Hessian, B,

We would like a measure of the overall coupling for the two vibrations, so we introduce the function $Fθ$ that integrates the total squared off-diagonal coupling over the span of vibrational motion in the harmonic oscillator ground state (i.e., from −1 to 1),

Evaluating the integrals gives,

where $mR=m2m1$ is the ratio of the masses of the center and terminal atoms. Solving for θ that minimizes F leads to two cases

Where the critical value for the anharmonicity is inversely proportional to the ratio of the masses $ζc=6mR2+4mR+3−1/2$⁠. Thus, in the case of nearly harmonic modes, the normal modes are preferred, whereas if strong anharmonicities are present, we obtain complete localization of the two modes. These results remain true even if all the modes of the system are included, though the form of $Fθ$ becomes more cumbersome.

This analysis shows that strongly anharmonic modes benefit from localization because the mixing of anharmonic vibrations in the normal modes introduces larger coupling over the region spanned by $Fθ$⁠. Further, the mass ratio is an important measure of harmonic coupling that arises in near-resonance vibrations. As noted by Child,⁵¹ two degenerate vibrations that are coupled will mix to form two decoupled normal modes with an energy difference proportional to the strength of the coupling. Thus, vibrations involving light atoms that are separated by heavier nuclei experience weaker coupling in the harmonic approximation. When considering strongly anharmonic vibrations like X–H stretches, the harmonic coupling introduced by localization is therefore less important than the decrease in anharmonic coupling. While Child’s metric justifies using local modes to describe overtone vibrations in X–H systems, our result shows this heuristic is generally applicable to any vibrations with significant mass ratios, m_R.

Partially localized vibrational modes with reduced intermode coupling can be constructed by combining a vibrational localization metric with a constraint on the coupling. First, the well-known Boys⁵² and Pipek-Mezey (PM)³⁶ orbital localization metrics are adapted to localize molecular vibrations.^28,29 PM orbitals³⁶ (vibrations) minimize the number of nuclei involved in a given orbital (mode), while the related Boys^38,39,52 metric maximizes the distance between the centers of the orbitals (modes). Consistent with the results of Steele and Reiher,²⁹ preliminary investigations found that using the atomic criterion (PM-like) to localize leads to vibrations that contain contributions from spatially separated nuclei, especially for larger systems (see Figure S1 of the supplementary material⁵³). Regardless of the proximity of those atoms, such projective metrics give equal locality to modes with the same number of atoms involved, so the distance (Boys-like) criterion, which avoids this issue, was chosen for subsequent use.

In the distance localization scheme, the center of each vibration, C_v, is calculated using

Where Q_iα,v is the magnitude of motion in the vth vibration for atom i in Cartesian coordinate α, and R_i is the position of the ith nucleus in the system. The localization metric, L, is then defined as the sum of the distances between pairs of all M modes,

Maximizing L is equivalent to maximizing the sum of the C_v terms.³⁶ In previous studies, L was maximized by performing pairwise rotations of the vibrational modes using a Jacobi sweep algorithm.^{28,29,36,39,52} During the sweeps, the normal mode vibrational basis, V₀, is rotated to yield the new modes, V, according to the unitary transformation over a vector of rotation angles, θ, giving

By transforming the basis, we also obtain a transformed Hessian matrix. This matrix, B, captures the harmonic bilinear coupling between the modified vibrational modes

Clearly, if V = V₀, then B is the diagonal matrix consisting of the eigenvalues. We measure the magnitude of off-diagonal coupling in the Hessian, using a scalar coupling constraint,

The form of S is closely related to the second-order energy correction in perturbation theory, namely, $E(2)=∑i<jHˆij2Ei(0)−Ej(0)−1$⁠. S is therefore an approximate measure of the energetic effect of intermode coupling. By choosing a value of ε², L can be maximized with controlled intermode coupling. The constrained optimization is achieved by evolving the following differential equation:

Where λ is the Lagrange multiplier whose value is determined by

To ensure smooth behavior around the constraint, we include a sigmoidal function, $gS$⁠, chosen such that g′(S) > 0 and $gε2=1$⁠. This allows L to increase so long as S ≤ ε². Furthermore, S is guaranteed to increase (or remain constant) while S ≤ ε², since it is observed that ∇ L ⋅ ∇ S ≥ 0. The sigmoidal function was chosen to be,

such that the values range from 0 to 2. This has the added effect of providing a “restorative” force if the numerical step causes S to exceed the constraint. The width of the sigmoid is determined by the parameter α. In practice, α ≈ 0.1 ⋅ ε² yields a sufficiently smooth derivative near S = ε² while allowing rapid convergence. To solve Equation (13), we employ an explicit first order method, namely,

is integrated until L is sufficiently converged.

All Cartesian Hessians were generated using the density functional ωB97X-D^54–57 and the double-zeta, polarized 6-31G** basis^58,59 in the Q-Chem⁶⁰ software package, on structures optimized at the same level of theory. Vibration optimization was performed in a stand-alone program written in C+ +. The convergence criteria for the localization metric, L, were specified as Δ L ≤ 10⁻⁵ for smaller systems, and for the larger systems (diphenylethane and stilbene), Δ L ≤ 10⁻⁴.

To being examining the relationship between vibrational locality and coupling, a tetramer of HF provides a simple model of weakly coupled anharmonic oscillators. The four monomers are aligned in a plane in parallel orientation with intermolecular spacing of 3.0 Å. Although the calculations include all 3N − 6 internal DOF, we restrict our discussion to the four stretching modes because they make the greatest contributions to the changes in L.

Figure 3 shows the four highest frequency normal modes (H–F stretches) and their localized counterparts. In the normal mode basis, these vibrations are delocalized over the outer (ν₁, ν₂) and inner (ν₃, ν₄) pairs of HF monomers. Delocalization occurs because of the weak harmonic coupling caused by long-range interactions. However, since HF has strong anharmonic character while the coupling between the resonant vibrations is small, this small decrease in harmonic coupling is outweighed by the increased higher-order coupling terms.

By systematically increasing the off-diagonal coupling constraint, all of the vibrational modes progressively become localized, as shown in Figure 3. The graph also shows the evolution of the diagonal elements $Hii1/2$ as S increases. At S = 0, the values of $Hii1/2$ equal the frequencies of normal mode vibrations. The first step of partial localization corresponds to a significant increase in locality, which results from unmixing the two lower-energy stretching modes. Because the frequencies of these modes are nearly degenerate (Δ ω < 0.2 cm⁻¹), this introduces negligible change in the $Hii1/2$ terms associated with ν₁ and ν₂. When S = 0.15 cm⁻¹, the stretches on the inner HF pair also localize. For the fully local modes, the deviation in $Hii1/2$ for these vibrations remains small (∼3 cm⁻¹).

In the HF tetramer, the low coupling between the four H–F stretches means that delocalized normal modes are a poor representation of vibrational motion. This intuition is quantified by monitoring the mode-mode coupling while the modes are being progressively localized, where the stretches can become highly localized with minimal coupling penalty. To gain a deeper understanding of the relationship between locality and coupling, more complicated systems are now examined.

An H-bonded dimer of water molecules is a commonly studied vibrational model with moderately coupled modes. Unlike the HF tetramer, each of the X–H stretches has a different frequency as a result of the local electronic environment. However, the mixing of stretches in the normal modes complicates the elucidation of these differences. For instance, the H-bonded O–H stretch requires two normal modes (Figures 4(B) and 4(D)) to describe the motion of that hydrogen atom. However, by localizing the stretching vibrations (Figure 4), the properties of each O–H stretch are isolated. After S = 15cm⁻¹, the four stretches are nearly identical to the fully localized representations (Figure 4, right) and the relative energy of vibration D′ is clearly stabilized by the H-bonding interaction compared to the other three O–H vibrations.

In the local mode representation, the diagonal terms are in agreement with expectations about the relative strength of each O–H bond, with the hydrogen bond donor (D′) being the weakest of the four and its neighboring hydrogen having the strongest bond (B′) when compared to the stretches on the H-bond-accepting water molecule. Additionally, since the total changes in $Hii1/2$ remain small, localizing these modes introduces minimal coupling. Minimal coupling upon localization, however, is not necessarily assured for the remaining vibrational modes in this system.

Two bending modes and six unimolecular rotational modes comprise the remaining internal degrees of freedom in the water dimer (Figure 5). In the normal modes (left), the molecular rotations are delocalized and the monomers simultaneously rotate along different axes (O, P). Alternately, only one molecule rotates while the other remains stationary in the corresponding local mode (O′). As in the X–H stretches, localizing these motions to single monomers occurs without introducing large amounts of harmonic coupling. This type of vibrational localization, where vibrations are isolated to individual monomers, should be highly beneficial for treating clusters of molecules.

In contrast to rotations in O, full localization of vibrations P and K yields modes with significant coupling and large shifts in $Hii1/2$⁠. Although the bending modes (J and K) near 1700 cm⁻¹ are already localized to single monomers in the normal mode representation, K begins to mix with the low-energy mode P as the coupling metric S is increased. At full localization, modes P and K have transformed into modes P″ and K′, which each consist of motion of a single hydrogen atom. The large change in $Hii1/2$ for modes K and $PΔHii1/2=−650cm−1and700cm−1,respectively$ requires the introduction of significant off-diagonal coupling (S = 50 cm⁻¹). Full localization of these modes also leads to a loss of the chemically intuitive bending motion, where in its place is a coordinate system corresponding to a nearly arbitrary vector localized on a single atom.

In contrast, partial localization allows sufficient locality to resolve the single-molecule rotations while preventing the formation of strongly coupled modes like K′ and P″. By S = 20 cm⁻¹, all of the rotations are localized to single molecules (resembling O′ and P′), while the bends (J and K) remain intact as localized intramolecular bends. At this point in the localization process, deviations in $Hii1/2$ remain small over all of the modes, indicating that this basis remains a reasonable description of the harmonic potential. Thus, the partially localized modes around S = 20 cm⁻¹ are easily interpreted modes with relatively low overall coupling. Importantly, by this stage, 96% of the maximum locality has been reached, indicating that the remaining 4% of locality is responsible for 60% of the coupling!

The 12-dimensional, closely interacting vibrations of ethene serve as the next case for careful examination. Figure 6 shows the partial localization process for the twelve vibrational modes of ethene. At S ≈ 1.3cm⁻¹, the four C–H stretching normal modes (far left) collapse into degenerate pairs of symmetric and asymmetric vibrations on the left and right methylene fragments, where these C–H stretches resemble the X–H stretching normal modes of the water dimer. Unlike the water dimer stretches, all of the C–H stretches experience identical environments, leading to totally degeneracy in fully localized ethene C–H stretches by S ≈ 7cm⁻¹.

While localization of the four C–H stretches onto the two CH₂ fragments occurs quite rapidly, only after S ≈ 32cm⁻¹ do the C–H bends fully merge (Figure 7). Because bending motions affect the forces on neighboring nuclei more strongly than stretches, they are expected to localize at larger values of S than the stretching modes. However, even for these more strongly coupled modes, a clear point of intermediate localization occurs at S = 5 cm⁻¹ where the two symmetric H–C–H bends become degenerate. The remaining three modes between 980 cm⁻¹ and 1070 cm⁻¹ correspond to out of plane hydrogen motion, which require a fourth degree of freedom to localize completely.

Vibrational localization is typically performed on a subset of the 3N Cartesian degrees of freedom. At minimum, exclusion of the rotational and translational dimensions is necessary for localization, and often a much smaller set of vibrations (e.g., only X–H stretches) are selected. By examining localization of all 3N degrees of freedom in ethene, instead of 3N − 6 or fewer, it will be shown that excluding degrees of freedom is a crucial step for otherwise unconstrained localization. When all 3N modes are included in the localization and S is unrestricted, the vibrational modes completely collapse to single-nucleus displacements. In other words, full locality consists of modes that are simply rotations of the Cartesian basis vectors on each atom. This single-nucleus locality results in significantly greater coupling and large deviations in $Hii1/2$ when compared to the harmonic basis. As shown in Figure 8, localization continues far beyond the S ≈ 40cm⁻¹ limit that maximized locality in the 3N − 6 vibrational basis.

Despite the major differences at maximum localization when compared to the 3N − 6 degree of freedom localization, the early behavior (S < 45 cm⁻¹) of the internal modes is consistent between the two models. To demonstrate this effect, we compare the localization of ethene vibrations in the internal, 3N − 6, basis and the full, 3N, basis (Figure 9). Total localization of the 3N basis leads to deviations of hundreds of wavenumbers in $Hii1/2$⁠, showing clearly that local modes are highly dependent on the chosen subspace upon which localization is performed. In the limit of a large system where 3N ≈ 3N − 6, localizing the full set of vibrational modes becomes similar to choosing an arbitrarily aligned atom-centered Cartesian coordinate system. Such an arbitrary coordinate system is essentially guaranteed to have large values of coupling, and little useful vibrational frequency information.

While localization procedures often begin from the 3N − 6 normal mode eigenvectors, the neither atomic nor distance localization metrics guarantee unique solutions. To demonstrate this, six new sets of modes for ethene were generated by rotating all pairs of normal modes using random values of theta $0≤θ<π4$⁠. By localizing with the distance criterion, these randomized initial conditions yield sets of modes where S varies by almost 60 cm⁻¹ (Figure 10). The normal-mode frequencies are presented on the y-axis (red bars) for comparison. Interestingly, the locality metric L varies by less than 0.01% between all calculations. Thus, while there are many solutions to the vibrational localization problem, the maximum attainable locality is intrinsic to the system. When comparing two representations with the same L, the one with a lower value of S will be the preferred solution.

Consistent with the behavior of molecular orbital localization procedures,^{37,43,61–64} the final set of localized modes are dependent on the path taken during optimization. Existing localization procedures directly maximize L, whereas the partial localization procedure follows a progressively relaxed constraint on S. For ethene, progressively relaxing the coupling constraint leads to a representation where S is 8 cm⁻¹ smaller than when localization is maximized directly from the eigenvectors. However, monitoring S only informs the extent of harmonic coupling, so we must develop a new metric to assess the anharmonic terms that localization seeks to minimize.

While normal modes are formed via minimization of off-diagonal coupling in the Hessian, this is only true for the geometry at which they are generated. Unless the PES is truly harmonic, even small changes in the geometry lead to finite coupling terms. Thus, a more meaningful measure of the quality of a vibrational basis will account for coupling in a region of space, akin to the use of $Fθ$ in the theory section.^65–67 Ideally, this spatial region will have physical relevance to the process being examined, which in this case is the local region of space accessed by the vibrations. As a specific metric, we consider a set of geometries that span the zero-point vibrational motion of the molecule. Since generating a full grid is prohibitively expensive, we limit our calculations to linear scans along the normal mode coordinates. Thus, N_H = (3N − 6)^∗n Hessian calculations are used to generate a grid of n Hessian calculations per mode (n = 11 for this investigation). For a pair of modes p and q in a basis localized for some value of S, the pairwise coupling strength, $CpqS$⁠, is defined as the root mean square of all of the elements $HpqS$

where $HS=VSTHVS$⁠.

The distribution of pairwise coupling values for ethene is presented in Figure 11 for a number of vibrational mode representations. The seven data sets include the NMs, three partially localized modes (S) generated in the 3N − 6 localization (Figure 6), LMs generated directly from the normal modes, and two of the random initialization distributions (R) with S = 107 and 116 cm⁻¹, respectively.

The first four distributions span the partial localization scan over S (Figure 6) and show that large changes to the locality of the vibrations can have small effects on the overall inter-mode coupling. Indeed, a great deal of localization is attainable without significant coupling cost, but ethene is a small enough system that the full benefits of locality are not realized. Indeed, randomized local modes have greatly increased coupling values (Figure 11 is on a log scale), suggesting local modes can even be detrimental in some cases. Thus, we investigate diphenylethane and stilbene as examples of large systems where localization of modes can lead to more significant spatial and energetic decoupling.

Thus far, we have examined the effect of partial localization for several small model systems, but the full benefits of this procedure will only be apparent with larger molecules. Diphenylethane and stilbene are useful examples because they are large enough that there is reduced coupling between nuclear motions on opposite sides of the molecule. Furthermore, their similar structure allows a close examination of the effect of increasing conjugation between the two phenyls.

Diphenylethane consists of two covalently connected phenyl groups (Figure 12, top) that appear to be ideal fragments for partial localization. Because of the relatively weak electronic coupling from the ethane linker, partial localization should show clear benefit early in the optimization process. This is indeed the case, as shown by the anharmonic coupling distributions for several values of S in Figure 13. As the vibrations of diphenylethane localize, the coupling distribution broadens, but aside from two mode pairs whose coupling remains close to the normal mode maximum coupling, at S = 75 cm⁻¹ the couplings are greatly decreased. A minimum for the interquartile range of anharmonic coupling is reached around S = 345 cm⁻¹, after which the coupling increases once again. As with the small molecules discussed before, most of the benefit of localization is gained well before the maximum of S is reached. In this case, the basis with the smallest anharmonic coupling terms occurs at 20% of the maximum value of S (L = 65%).

Partial localization also allows mode pairs to be distinguished by spatial separation. While the distance between mode centers in the normal modes is typically small (Figure 14(a)) due to their high degree of delocality, partial localization significantly increases the average distance between mode centers (Figures 14(b) and 14(c)). Interestingly, while the distance between mode pairs increases rapidly at first, there is little change in the distribution after S = 345 cm⁻¹ (see supplementary material Figure S2⁵³). This is consistent with the relationship between L and S for the smaller molecular examples. Furthermore, the short-range nature of vibrational interactions leads to an inverse relationship between mode pair distance and inter-mode coupling strength. Figures 14(d)-14(f) partition the coupling by pair distances, where the strength of coupling within the 2.0 Å and 4.0 Å distances emphasizes the relationship between locality and coupling in the various vibrational mode representations. In the delocalized normal modes, the majority of the pairs (42%) are near neighbors and as a result, experience moderate coupling. The coupling decay is immediately apparent in Figure 14(e) where the majority of the mode pairs (41%) are simultaneously weakly coupled and well separated. This distribution is ideal for anharmonic calculations because many of the weakly coupled terms can be identified by distance alone, without having to explicitly calculate anharmonic coupling terms. While the totally localized modes (S = 1706 cm⁻¹) have a larger number of mode pairs separated by more than 4.0 Å, the coupling strength for mode pairs closer than 2.0 Å shifts up by a factor of ∼5 compared to the minimum coupling at S = 345 cm⁻¹.

Stilbene presents a unique challenge to vibrational localization since the full conjugation makes it difficult to assign decoupled fragments. While the spatial separation (as with diphenylethane) would imply that the vibrations on the two phenyl rings are only moderately coupled, the extended π system makes it unclear exactly how the vibrations can be partitioned. Despite this challenge, partial localization analysis shows that intermediate localization is indeed beneficial in this system (Figure 15).

Especially at small values of S (∼125 cm⁻¹), the overall decrease in coupling for stilbene is even more favorable than in diphenylethane. Because of the size of the molecule, Figure 16 only presents a few of the vibrational modes in stilbene, but the full $Hii1/2$ plot is shown in the supplementary material (Figure S3⁵³). The localization of the C–H stretches behaves identically to the previously discussed systems, with immediate localization followed by very little change in the modes until the maximal coupling is allowed.

The effects of partial localization on the distribution of coupling terms (Figure 17) are quite marked. Even in this strongly coupled system, partial localization provides a significant ability to lower coupling values compared to normal and local modes. The partially localized (50% of L_max) representation at S = 128 cm⁻¹ (8% of S_max) has the lowest overall anharmonic coupling distribution. As with diphenylethane, the formation of long distance tail provides a predictive metric for inter-mode couplings to be excluded from subsequent calculations.

This paper presents a detailed analysis of the transition between delocalized normal modes and fully localized modes by sequentially maximizing mode locality with constrained coupling. Particularly for large systems, the suggested algorithm provides a useful metric for determining how much locality is possible without introducing large couplings into the vibrational Hamiltonian.

As seen in Figure 18, the majority of localization (>75%) is achieved before 30% of the maximal coupling for each system, which provides a useful metric to determine the optimal degree of localization. Alternatively, if the main concern is eliminating coupling terms from the potential energy expansion, identifying the basis where the distance distribution stops evolving is also a viable choice (see Fig. S2⁵³). Partially localized modes are shown to reduce the square total Hessian off-diagonal coupling terms in the vicinity of the minimum energy geometry, making them highly useful as a basis in vibrational structure simulations.³⁵ The method is particularly advantageous in that it requires the selection of a single parameter to perform localization of all vibrational modes rather than a heuristically chosen set of modes, as has been frequently employed in previous methods.

Additional analysis of fully localized modes found that they are non-unique and can depend on the initial randomized mixing of normal modes. The maximization of the locality objective function, L, therefore depends on the path taken. By taking the path that maximizes L with constrained coupling, the presented algorithm yields an improved fully local mode basis with significantly reduced off-diagonal coupling.

The analysis of this manuscript suggests that the benefits of localization come from decreasing anharmonic coupling at the expense of increasing harmonic coupling. If a potential is truly harmonic, delocalization would be desirable because it produces the simplest possible description of the vibrational motion. Essentially, all molecular systems, however, have some degree of anharmonicity in the potential, so it is necessary to counterbalance the quality of the harmonic approximation against the coupling produced by higher-order terms. Visualizing the relationship between Hessian diagonal elements and coupling for the span ranging from normal modes to a fully localized basis, such as shown in Figures 3–6, demonstrates the intrinsic difficulty in providing compact local mode representations while maintaining reduced values of coupling.

Partially localized modes provide a straightforward connection between the various bases used in anharmonic structural analysis and produce results similar to those from direct minimization of the intermode coupling as seen in a recent paper from our group.⁶⁸ By bridging the normal modes and the fully localized modes, partial localization provides new insight into why the local modes often provide a better basis for treating anharmonic systems.^28,29 Similar to the results of Hirata et al.,²¹ we note that the quality of a vibrational basis is dependent on balancing the magnitude of harmonic and anharmonic coupling terms inherent in the representation. The tendency of certain vibrational bands to localize before others also helps to explain the success of certain heuristic-based partial localizations such as those employed by Sibert and co-workers on C–H^47,48 and carbonyl⁶⁹ stretches.

Regardless of the possible benefits of partially localized modes, certain vibrations, such as internal rotations, will not be well described in Cartesian coordinates and significant mode-mode coupling will exist regardless of the degree of locality introduced. A partial localization algorithm that uses a curvilinear coordinate system would further reduce coupling and may improve the overall performance of this method. The optimization of a local mode basis that includes rotational or floppy vibrations is a natural area of further study.

The authors would like David Braun for his continued support in operating our computational resources. We also thank the National Science Foundation for the support given under Award Nos. DMS 1436075 and MPS 1551994.