# How is the bond length calculated from the total electronic energy?

In the quantum chemistry course that I currently attend, it was said several times that one of the key quantities derived from molecules by means of computation is the total electronic energy $E_\text{el}$, because from this many derived quantities can be calculated, such as bond length or $J$-coupling values for NMR experiments.

However, I've never gotten a good explanation as to how to proceed when I've calculated a value for $E_\text{el}$ and want to find out the bond lengths in the molecule in question.

What are the formulae and/or computation steps necessary for this?

The traditional way of getting the right geometry for your system is by taking a fixed atomic structure, perform an electronic structure minimization, calculate the forces acting on the atoms from that, use them to alter the atomic structure according to the forces and then use the altered structure to start over again. This is repeated until you converge to the equilibrium geometry.

However there is an alternative (quite elegant) way for DFT calculations using pseudopotentials or the PAW method with plane-wave basis sets: the ab-initio molecular dynamics method introduced by Roberto Car and Michele Parrinello. The "killer feature" of the Car-Parrinello method is that it allows the simultaneous solution of the electronic structure problem and the equations of motion for the nuclei which allows for the relaxation of the nuclei to find stable structures, i.e. the equilibrium geometry, as well as for thermal simulations of solids and liquids - so this method does offer much more than a simple geometry optimization. How does it achieve that? In the Car-Parrinello approach, the total Kohn-Sham (KS) energy is the potential energy as a function of the positions of the nuclei. The forces from this energy then determine the molecular dynamics (MD) for the nuclei. The special feature of the Car-Parrinello algorithm is that it also solves the quantum electronic problem using MD. This is accomplished by deriving the equations of motion from a fictitious Lagrangian which contains a fictitious kinetic energy of the electronic (KS) wave functions and some constaints to ensure that the wave functions remain orthonormal:

\begin{align} \mathcal{L} &= \sum_{n} f_{n} m_{\psi} \langle \dot{\psi}_{n} | \dot{\psi}_{n} \rangle + \frac{1}{2} \sum_{i} M_{i} \dot{\vec{R}}_{i} \! {}^{2} + E_{\mathrm{DFT}}[\psi_{n}, \vec{R}_{i}] - \sum_{m, n} \Lambda_{n,m} \left( \langle \psi_{m} | \psi_{n} \rangle - \delta_{n,m} \right) \end{align}

where $f_{n}$ is the occupation of the $n^{\mathrm{th}}$ Kohn-Sham eigenstate $\psi_{n}$, $M_{i}$ and $\vec{R}_{i}$ are the mass and the position of the $i^{\mathrm{th}}$ nucleus, $\delta_{n,m} = \begin{cases}1 & n=m \\ 0 & n \neq m \end{cases}$ is the Kronecker delta, $\dot{\psi}_{n} = \frac{\mathrm{d} \psi_{n}}{\mathrm{d}t}$ and $\dot{\vec{R}}_{i} = \frac{\mathrm{d} \vec{R}_{i}}{\mathrm{d}t}$ are the first derivatives with respect to time and $\Lambda_{n,m}$ is a Lagrange multiplier ensuring the constraint of wave function orthonormality. The first term in this Lagrangian represents the fictitious kinetic energy of the wave functions. There is no physical correspondence for this term. Ideally one would choose the fictitious mass $m_{\psi}$ of the wave functions equal to zero. This introduction of this unphysical quantity also lead to naming the Lagrangian fictitious. The second term describes the classical kinetic energy of the nuclei. The third term $E_{\mathrm{DFT}}[\psi_{n}, \vec{R}_{i}]$ is the density functional total energy, which is a functional of the electronic wave functions and the atomic positions. The last term is the constraint of orthonormal wave functions $| \psi_{n} \rangle$. Applying the principle of least action to this Lagrangian leads to the following Euler-Lagrange equations (equations of motion for the ab initio MD simulation):

\begin{align} M_{i} \ddot{\vec{R}}_{i} &= \underbrace{- \nabla_{\vec{R}_{i}} E_{\mathrm{DFT}}}_{= \, \vec{F}_{i}} \\ m_{\psi} | \ddot{\psi}_{n} \rangle &= - \underbrace{ \frac{1}{f_{n}} \frac{\delta E_{\mathrm{DFT}}}{\delta \langle \psi_{n} |} }_{= \, \hat{H} | \psi_{n} \rangle } + \underbrace{ \sum_{m} | \psi_{m} \rangle \frac{\Lambda_{n,m}}{f_{n}} }_{\tilde \, | \psi_{n} \rangle \epsilon_{n}} \end{align}

where $\ddot{\psi}_{n} = \frac{\mathrm{d}^{2} \psi_{n}}{\mathrm{d}t^{2}}$ and $\ddot{\vec{R}}_{i} = \frac{\mathrm{d}^{2} \vec{R}_{i}}{\mathrm{d}t^{2}}$ are the second derivatives with respect to time, $\hat{H}$ is the Kohn-Sham hamiltonian, $\epsilon_{n}$ is $n^{\mathrm{th}}$ Kohn-Sham eigenvalue and $\vec{F}_{i}$ is the force acting on the $i^{\mathrm{th}}$ nucleus. The first set of equations is just Newton's equations of motion for the nuclei moving under the forces derived from $E_{\mathrm{DFT}}$.

The stationary solution of the second set of equations is equivalent to the Kohn-Sham equations, since for a steady state all time derivatives vanish so that:

\begin{align} \hat{H} | \psi_{n} \rangle &= \sum_{m} | \psi_{m} \rangle \frac{\Lambda_{n,m}}{f_{n}} \end{align}

Constructing a matrix $\mathbf{\Lambda}$ from the Lagrangian multipliers $\Lambda_{n,m}$ shows that $\mathbf{\Lambda}$ is proportional to the transpose of the matrix representation of $\hat{H}$, i.e. $\Lambda_{n,m} \propto H_{m,n}$. Diagonalizing $\mathbf{\Lambda}$ leads to the eigenvalues $\epsilon_{n}$ of the Kohn-Sham equations. So, this Lagrangian does not lead to the time-dependent Schroedinger equation for the electronic wave functions. Rather it creates a dynamic that (in the ideal case) keeps the electrons always in their ground state.

• Under which conditions does the Car-Parrinello method break down? Excited state analysis already seems to be out of the question... – tschoppi Oct 11 '14 at 21:57
• @tschoppi I've only seen it being used with plane-wave DFT methods. But it is more generally applicable than this. I think you can even use it with a CI calculation but the CP-method is rather expensive since it makes extensive use of forces which are usually quite costly to calculate. As for excited states: There are papers indicating it can be done (see for example here and here). – Philipp Oct 12 '14 at 0:47

The geometry of a system is what is known (or given) and a corresponding energy is computed. If you have an energy then you should have the geometry (or in your case, bond length).

The wave function, $\Psi$, in the Schrodiner equation is simply a function of all the positions of the particles within the system (nuclei and electrons). Typically chemists invoke the Born-Oppenheimer approximation which assumes that the nuclei are stationary (frozen in space) with respect to the very 'fast-moving' electrons. Therefore, we freeze the nuclei in space and compute the energy of this stationary state. Move the nuclei in your molecule a little bit and you will get a new energy.

As Wildcat was explaining, optimization procedures essentially moves the nuclei in your molecule around in space, computing an energy for each configuration. The goal of this type of routine is to find an energy that is lower than all the rest. Take, for example, the $\ce{H2}$ molecule. It contains one bond length and therefore the only coordinate we need to consider in an 'optimization procedure' is the internuclear separation of the two $\ce{H}$ atoms. It is the chemist which will choose the initial positions of these atoms. Say we initially set these atoms to be 200 pm apart. We compute an energy. Now, can we do better? (i.e. can we perturb this geometry in such a way that gives rise to a lower energy?) Well we can certainly try other bond lengths so we now set our bond length to 150 pm and compute a new energy. Lo and behold we see that this new energy is lower than the first! We keep trying different bond lengths until we reach a point where the energy increases no matter what we do with the separation (whether we increase it or decrease it). Once we are at this point, we have found the lowest energy configuration of $\ce{H2}$. Of course optimization routines are smarter than doing a brute-force approach to sampling all the possible quantities each geometrical parameter can take on. Furthermore, second-derivatives are essential for characterizing the nature of the stationary state that was determined (i.e. is the stationary state a true minimum, a transition state or some higher order saddle point?).

• So if I got this correctly: It is not from the calculated quantity, but by way of calculating it, that we find the bond length? – tschoppi Oct 8 '14 at 15:14
• The bond length is what 'you' make it. From that you compute an energy. Of course optimization procedures will take what you give it and manipulate the geometry as it tries to find the lowest energy. The geometry information is printed to the output file for every iteration the procedure takes. – LordStryker Oct 8 '14 at 19:09

There is no formula. In short, for an electronic state of interest you do what is known as the geometry optimisation, i.e. for the chosen electronic state you solve the electronic Schrödinger equation for a number of fixed (or clamped) nuclear configuration chosen in a certain way by the employed minimisation procedure with the aim to find the minimum of the electronic energy. And the nuclear configuration that corresponds to minimum of the electronic energy for the chosen electronic state describes the equilibrium geometry of a molecule being in the chosen electronic state. Then from such equilibrium nuclear configuration you can easily extract structural parameters, such as, for instance, bond lengths.

• NMR shielding is a bit more complicated story. You need to become familiar with the basics of the so-called response theory, to understand how some electric and magnetic properties, such as NMR shielding, are calculated. – Wildcat Oct 8 '14 at 11:57