Equivalence of Electron-Vibration Interaction and Charge-Induced Force Variations : A New O ( 1 ) Approach to an Old Problem

Calculating electron-vibration (vibronic) interaction constants is computationally expensive. For molecules containing N nuclei it involves solving the Schrödinger equation for O(3N) nuclear configurations in addition to the cost of determining the vibrational modes. We show that quantum vibronic interactions are proportional to the classical atomic forces induced when the total charge of the system is varied. This enables the calculation of vibronic interaction constants from O(1) solutions of the Schrödinger equation. We demonstrate that the O(1) approach produces numerically accurate results by calculating the vibronic interaction constants for several molecules. We investigate the role of molecular vibrations in the Mott transition in κ-(BEDT-TTF)2Cu[N(CN)2]Br.


Introduction
Electrons in molecules and solids feel two forces: a direct Coulomb interaction between themselves and a Coulombic interaction with the atomic nuclei.As the vibrations of the nuclei are quantized the interaction with the nuclei takes the form of the vibronic interaction.Electron-vibration interactions play important roles across science.In physics vibronic interactions can give rise to superconductivity [1], spin and charge density waves [2], polaron formation [3] and piezoelectricity [1].In chemistry they impact electron-transfer processes [4], Jahn-Teller effects [5], spectroscopy [6,7], stereochemistry [5], activation of chemical reactions [5] and catalysis [5].In biology the vibronic interactions play important roles in photoprotection [8,9], photosynthesis [10] and vision [11].In engineering vibronic interactions play important roles in determining the efficiency of organic photovoltaic cells and organic light emitting diodes [12].It is therefore clear that one of the central tasks for condensed matter theory and theoretical chemistry is to accurately and efficiently calculate vibronic interaction constants.
The Holstein Hamiltonian [3], which describes the vibronic coupling, is ) where H 0 e is the Hamiltonian of the electrons in the absence of vibrations, ω i is the frequency of the mode, â( †) i annihilates (creates) a vibron in the i th mode and ĉ( †) µ annihilates (creates) an electron in the state µ.The coupling between the i th vibrational mode and the µ th electronic state is given by where λ µ is the energy of the µ th electronic state.Therefore, to calculate the vibronic interaction constants via the frozen phonon method, one begins by calculating the electronic structure and the normal modes of the system; one then makes a small displacement of the system along each normal coordinate, and thus calculates ∂λ µ /∂Q i .Therefore, if there are N normal modes, N such calculations must be performed.The frozen phonon method can be speeded up by the use of density functional perturbation theory (DFPT) or linear-response methods [13].However, such techniques still scale as O(N ), albeit with a smaller coefficient than the frozen phonon calculation.An alternative is to consider Janak's theorem [14] which states that where E is the total energy of a fictitious system consisting of the system of interest with a small extra charge in the µ th electronic state and n µ is the electronic occupancy of the state µ.It is important to realize that, within the density functional formalism n µ is not required to be an integer [14].It follows from Equations 2 and 3 that where (∂E/∂Q i ) x indicates that the derivative is taken after the charge is changed by x relative to the charge of the initially optimized geometry.We have used the fact that for a (meta)stable geometry Conceptually it is useful to observe that in the frozen phonon and DFPT approaches one considers the problem from the point of view of the electrons.That is, one calculates the vibronic interaction constants by making a small perturbation to nuclei and considering the resultant change on the electron.However, in the O(1) approach one asks the question what happens to the nuclei when one puts a small additional charge in the system.In many cases one is only interested in the coupling of the vibrations to the electronic states closest to the Fermi level, however, if one were interested in interactions with other states, an extension of Janak's theorem to other electronic states allows for the determination of these parameters straightforwardly within the O(1) method.
The primary observation reported here is that using Equation 4 the g i can be calculated by solving a single nuclear geometry.However in addition to pointing toward an efficient means for the calculation of the vibronic interaction, Equation 6 demonstrates a fundamental relationship between the electron-vibrational interaction and the charge induced geometrical relaxations.The former is generally viewed as a quantum-mechanical characteristic while the latter is viewed as a classical phenomena.Equation 6 simply states that the changes in the atomic Hellmann-Feynman forces due to charge addition or depletion depend on the same fundamental quantities as the electron-vibrational interaction.Below, we present results from this new O(1) method for calculation of the spin-vibron interaction (∂λ µ /∂Q i ).

Results
We now show that Equation 4 reproduces the values of the vibronic interaction constants calculated by the frozen phonon method for a number of small molecules.As discussed in the methods section, we have implemented our scheme for calculating electron-vibrational interaction constants within a densityfunctional-based computational code written by  using the Perdew-Burke-Ernzerhof [19] exchange correlation functional.
It is worthwhile to briefly comment on some of the approximations used during these numerics, which are not, formally, required by the method itself.All of the calculations below are performed within the Born-Oppenheimer approximation (BOA).It may not be immediately clear that g iµ is finite in the BOA.The BOA is to assume that the many-body wavefunction Ψ({R}, {r}), which in general depends on the set of nuclear co-ordinates, {R}, and the set of electronic co-ordinates, {r}, is separable into a nuclear wavefunction ϕ({R}) and an electronic wavefunction ψ({R}, {r}).That is, Ψ({R}, {r}) = ϕ({R}) ⊗ ψ({R}, {r}).Physically this corresponds to two assumptions: (i) the electronic wavefunction depends on the nuclear positions, R j , but not the momenta of the nuclei and (ii) the nuclear motion sees a smeared out potential from the speedy electrons.We then treat the nuclei as classical particles and so we only consider ψ({R}, {r}), and hence we calculate only the electronic density ρ({R}, {r}) within the DFT framework.It is the dependence of ρ({R}, {r}) on the nuclear coordinates (but not the nuclear momenta) that leads to finite g iµ .Finally we note that the BOA does not give the exact solution to Equation 1 unless all of the g iµ vanish.However for small g iµ the BOA gives a good approximation to, for example, the exact solution of the two site Holstein model [3].
In Figure 1 we plot the calculated vibronic interactions for tetrathiafulvalene (TTF), tetracyanoquinodimethane (TCNQ), bis(ethylenedithio)-tetrathiafulvalene (BEDT-TTF), C 60 and a number of diatomic and other small molecules calculated by the frozen phonon method against the vibronic interaction constants calculated by our method.It can be seen that the values calculated by either method are in good agrement with one another.

Figure 1.
Comparison of the electron-vibrational interaction constants for a variety of molecules calculated by the order-one method, g1 , and the frozen phonon method gfph (both in Ha/Bohr; giµ = g iµ √ 2m i hω 3 i ).It can be seen that the agreement between the frozen phonon and order-one methods is excellent.In the inset we plot the same data on a log-log scale which displays the excellent agrement achieved even for small values of g.For each molecule we plot the electron-vibrational interactions between the HOMO and LUMO levels and every vibrational mode.

~~Ñ
ote that, as the normal coordinate is the eigenvector of the dynamical matrix, the sign of Q i is not well defined.Therefore, although ∂E is negative definite, the sign of g iµ is also not well defined.However, the product g iµ ∂Q i is well defined and describes the distortion caused by the addition of a charge.For example, for the stretching mode in a simple diatomic molecule, the sign of the product g iµ ∂Q i indicates whether the bond is lengthened or shortened when a charge is added.
Because the calculation of vibronic interaction constants is now computationally inexpensive, it is possible to use the vibronic interaction to estimate other quantities.For example, the ionization energy (or electron affinity) of a molecule, as opposed to the vertical ionization energy (electron affinity), may be extracted from the information already calculated.To do this one determines the Hubbard U (a second derivative) from either the energy as a function of ∓δn µ and ∓2δn µ or from the appropriate eigenvalue as a function of ∓δn µ .As the vibronic interaction constants are known, the adiabatic ionization energy (electron affinity) can be calculated directly from Marcus-Hush theory [4].In Figure 2 we plot the difference between the adiabatic ionization energies (electron affinities) calculated in this way and the correct result (within the PBE framework) for the molecules considered in Figure 1 as a function of the the number of the atoms in the molecules.While for small molecules this approximation is not accurate, the results for the larger molecules are excellent.Clearly, it is for larger molecules that this approximation is needed, as the computational power required to relax the geometry grows with the size of the molecule.
Figure 2. The ionization potentials and electron affinities of the molecules considered in Figure 1.We plot the difference in energy, ∆E between the correct result (within the PBE approximation) and the results from approximating the ionization potential (electron affinity) by estimating the Hubbard U and relaxation energy from the O(1) method.It can be seen that the approximation improves dramatically as the number of atoms in the molecule, N , increases.The inset show the same data on a log-log scale, this suggests that the error in the approximation decreases as ∆E ∼ N −α with α ∼ 1.5-2.

∆E (eV) ionisation potential electron affinity
It is important to stress that the calculation of the vibronic interaction constants is, in fact, intrinsically an O(1) problem as the dynamical matrix and the forces calculated on the geometries used to calculate the dynamical matrix together with the eigenvalue changes determined during the course of these calculations contain all the information required to calculate the g iµ .However, significant computational complexities must be overcome to retrieve this information.
In the course of the above calculations, a particular practical benefit of the O(1) method became clear.It is well known that the calculated electron-vibrational interaction is rather sensitive to the size of the displacement used in a frozen phonon calculation.In general, a simple criterion for the size of the displacement (e.g., x = √ 2∆/k, where ∆ is a small energy and k is the spring constant associated with the vibrational mode in question) does not produce uniformly reliable results for all possible vibrational frequencies in the frozen phonon calculations of the type discussed above.However, we have found that the results for the O(1) method are remarkably insensitive to the value of the "small" charge, δn µ , used in the calculation.
The calculation of isotope effects is an O(0) problem in our method.That is, only trivial matrix manipulation and no further solutions of the Schrödinger equation are required to calculate the vibronic interaction constants as the isotopic masses are varied.This contrasts to the frozen phonon method where, as the dynamical matrix and thus the normal modes change upon isotopic substitution, O(N ) additional calculations are required.
A particularly interesting isotope effect is observed in the superconductor κ-(BEDT-TTF) 2 Cu[N(CN) 2 ]Br.When the eight hydrogen atoms in the BEDT-TTF molecule are replaced by deuterium, the ground state is found to be a Mott insulator [20,21].Recent efforts to calculate the Hubbard U [22][23][24] and intramolecular hopping integrals [23][24][25] for this material do not provide an explanation of this.This material is part of a larger family: κ-(BEDT-TTF) 2 X.Other well known members of the family included κ-(BEDT-TTF) 2 Cu[N(CN) 2 ]Cl, which is a Mott insulator at ambient pressure and is driven metallic/superconducting by hydrostatic pressure, and κ-(BEDT-TTF) 2 Cu(SCN) 2 which is metallic/superconducting at ambient pressure [21].This has led to the idea of a unified phase diagram where "chemical pressure" (or hydrostatic pressure) can be used to tune the system across the Mott metal-insulator transition.Thus it is believed that κ-(BEDT-TTF) 2 Cu[N(CN) 2 ]Br is very close to the Mott transition [20,21,26] and that deuteration produces a small negative chemical pressure [20,21,27,28].However, the microscopic mechanism for this small negative chemical pressure (i.e., a small decrease in the ratio W/U , where W is the bandwidth and U is the effective repulsion between two electrons on the same dimer, (BEDT-TTF) 2 ) has not previously been explained.
In the crystal the highest occupied molecular orbital (HOMO) of the BEDT-TTF molecule is doped with holes from the anion layer, therefore one expects [3] that, in the small polaron limit, on deuteration, the bandwidth will change by a factor where W D (W H ) is the bandwidth of the deuterated (hydrogenated) system and g ih(D) (g ih(H) ) is the coupling between the i th mode and the HOMO in the deuterated (hydrogenated) molecule.The calculated frequencies and vibronic coupling constants are tabulated in Table 1.Thus polaronic band narrowing provides precisely the small negative chemical pressure needed to explain the Mott transition by dueteration in κ-(BEDT-TTF) 2 Cu[N(CN) 2 ]Br.
Table 1.Calculated vibrational frequencies, ω i , masses, M i , and vibronic coupling to the HOMO, g ih , for protonated and dueterated BEDT-TTF.

Methodological Details
To implement our method one simply recalculates the electronic structure with a small change in the charge in the orbital of interest.This allows the electronic orbitals to relax to account for the additional charge.However, in this process, the nuclear structure remains fixed as that of the equilibrium geometry of the charge neutral system.The forces on the molecule can be calculated using the Hellmann-Feynman theorem.As the dynamical matrix and vibrational eigenvectors are already known from the calculation of vibrational spectrum, the individual vibronic interaction constants can readily be calculated by projecting onto the basis of normal modes.To perform the electronic structure calculations we have followed the ansatz that the Kohn-Sham wavefunctions can be expanded in terms of a set of atom-centered gaussiantype orbitals.The methods for performing these calculations are discussed in Refs.[16][17][18][29][30][31][32].
The directional vibrational derivative (Equation 4) is determined in a standard way by taking the scalar product the energy gradient with the vibrational eigenvector.Provided a variationally optimized quantum-mechanical wavefunction is determined, Hellmann and Feynman have shown that this energy gradient associated with an atomic displacement is exactly equal to the electric field due to all other electronic and nuclear charges in the system.When finite-basis sets are used the energy gradient has an additional term, commonly referred to as the Pulay correction, which is required to determine the correct energy gradient.We emphasize that this additional term should be and is contained in our calculations.See Reference [16] for discussion and references regarding Hellmann-Feynman-Pulay forces.

Conclusion
Prior to concluding, a few additional remarks are appropriate with regard to assessing the simplifications garnered from the O(1) approach.First, it is worthwhile noting that if one determines the electron-vibrational interaction by performing 3N electronic structure calculations, it is in fact possible to determine the vibronic interactions associated with all the electronic states rather than only the LUMO and HOMO levels.If vibronic interactions from levels other than the LUMO or HOMO are desired, the approach discussed here can be modified by varying the charge of the the particular electronic state of interest.Since one is generally interested in vibronic interactions associated with frontier orbitals, the charge-variation approach discussed here should be more efficient for most applications.A more general statement is that the scaling associated with the charge variation method scales linearly with the number of electronic states needed for the model Hamiltonian.In the unlikely event that an electronic system has O(3N ) interesting electronic states, the traditional approach may then be as efficient as the O(1) approach.
A second issue is related to the total computational cost since the O(1) method requires that the vibrational modes be determined.In our work we have determined the dynamical matrix by numerically determining derivatives of the interatomic forces as a function of atomic displacement.This method, while reasonably standard, requires O(3N ) electronic structure calculations in the most general case.Therefore the most pessimistic assessment is that, given the equilibrium geometry, the determination of the vibronic interaction still requires O(3N ) calculations since one needs this many calculations to determine the vibrational eigenvectors.However, this should then be contrasted with the O(6N ) calculations required to obtain the same information by tradition approaches.One counterpoint to this assessment is that the dynamical matrix is required for other phenomena such as infrared and Raman spectra and one general goal should be to determine efficient post processing algorithms to determine the interaction between a bath of vibrons and the interesting aspects of their environment.
Another issue is the widely held hypothesis that vibrational energies for large systems can generally be determined faster than O(3N ) by using molecular dynamics followed by a Fourier transform of the velocity-velocity autocorrelation function.Once the vibrational frequencies are determined, a Fourier transform of the 3N -dimensional trajectory also allows for the determination of, at least, the non-degenerate vibrational eigenvectors.Coupled with such a strategy, it may indeed be possible find the vibronic interactions of an exceedingly large system with scaling below O(3N ) even if the cost of determining the vibrational degrees of freedom are included.However, the most rigourous and relevant statement is that the vibronic interaction may now be determined as quickly as the vibrational spectrum.
In summary we have demonstrated that by considering the forces on the nuclei due to the addition or subtraction of an arbitrarily small electronic charge one may calculate the vibronic interaction constants as an O(1) problem.This method was shown to be numerically accurate for a large number of small molecules.