# How should I think about potential energy surfaces for reactions of complex systems?

I've recently started doing some computational chemistry, studying the reaction dynamics of some moderately complex organic and inorganic reactions. I'm doing things like finding maxima, minima, TS's, saddle points, geom. optimizations, and finding conformational isomers of all of the above. I'm really enjoying myself.

One thing I would like to better understand, though, is how to think of the potential energy surface (PES) of these sorts of systems. The way I see it, a PES is just a potential energy profile (PEP) with one or more dimensions beyond energy and reaction coordinate. I understand that you cannot really compute a complete PES for a complex system. Most of the PES's I've seen are for simple systems like this solvent study of an SN2 reaction:

(R. Otto et al. 2011, Nature Chemistry 4, 534-538)

This is simple enough to understand. They have added the dihedral angle as another dimension of their picture of the reaction dynamics. But what about a large system, especially a non-rigid one? Since there are so many more possible ways the system can move, does that mean its 'true' PES is polydimensional? I like to think of these mechanistic studies of big reactions as us 'exploring the PES', even though we cannot see the whole thing. How should I think about this surface? Am I understanding this concept properly?

• Yes, in general the PES has many dimensions and only some of them can be displayed in any one diagram. May 28 '16 at 11:55
• Is the dimensionality of the PES related to the number of motions, i.e.; 3N-6 for a given molecule? How should I think of the PES in the context of my computations? For example, sometimes I will look for a saddle point. Obviously this saddle point must exist on a surface of some sort. How should I visualize a saddle point on a polydimensional surface? May 28 '16 at 12:09

In order to understand what the PES really is, it is useful to see where this concept comes from.

A system of $M$ nuclei and $N$ electrons should be described by a multidimensional wavefunction $\Psi(\vec{r}_1,\dots,\vec{r}_N,\vec{R}_1,\dots,\vec{R}_M)$, where $\vec{r}$ denotes electronic coordinates and $\vec{R}$ denotes nuclear coordinates (we neglect spin for simplicity). Within the Born-Oppenheimer approximation, the total wavefunction is split in an electronic part $\phi$ and a nuclear part $\chi$: $$\Psi(\vec{r}_1,\dots,\vec{r}_N,\vec{R}_1,\dots,\vec{R}_M) = \phi(\vec{r}_1,\dots,\vec{r}_N;\vec{R}_1,\dots,\vec{R}_M)\chi(\vec{R}_1,\dots,\vec{R}_M)$$ where the electronic wavefunction depends parametrically on nuclear positions (as it will be clear in a moment).

If you inject the Born-Oppenheimer ansatz in Schrödinger equation for a system of interacting nuclei and electrons and you integrate out electronic degrees of freedom you get an equation for the nuclear wavefunction only $$\left[ -\sum_I \frac{1}{2M_I} \nabla^2 + \sum_{I<J} \frac{Z_IZ_J}{|\vec{R}_I-\vec{R}_J|} + E(\vec{R}_1,\dots,\vec{R}_M) \right]\chi = E_\text{tot} \chi$$ where $E(\vec{R}_1,\dots,\vec{R}_M)$ is the solution of the electronic problem $$\left[ -\sum_i \frac{1}{2}\nabla_{i}^2 + \sum_{i<j} \frac{1}{|\vec{r}_i-\vec{r}_j|} -\sum_{i,I} \frac{Z_I}{|\vec{r}_i-\vec{R}_I|} \right]\phi = E(\vec{R}_1,\dots,\vec{R}_M)\phi$$ As you can easily see the solution of the electronic problem (both the energy and the wavefunctions) depend parametrically on nuclear coordinates. The potential energy surface is the potential energy appearing in the equation of the nuclear wavefunction, i.e. $$U(\vec{R}_1,\dots,\vec{R}_M)=E(\vec{R}_1,\dots,\vec{R}_M) + \sum_{I<J} \frac{Z_IZ_J}{|\vec{R}_I-\vec{R}_J|}.$$ It is clear that finding the PES is a formidable task, since you have to solve the electronic problem for each possible set of nuclear coordinates. So, yes the PES is a multi-dimensional function of nuclear coordinates. If you have $M$ nuclei, your PES will be a function of at least $3M-6$ coordinates (once you remove translation and rotation of the whole system).

Since a multi-dimensional function is impossible to visualize, usually we study the PES by considering few interesting degrees of freedom (as the O-C distance and O-C-I angle in your example). Therefore, as you say, you can think as if you are exploring the PES locally and in a reduced space.

Since the nuclear wavefunction is sharply picked, (heavy) nucley can usually be considered as classical particles. Thus instead of the nuclear quantum Hamiltonian $$H_n = -\sum_I \frac{1}{2M_I} \nabla^2 + \sum_{I<J} \frac{Z_IZ_J}{|\vec{R}_I-\vec{R}_J|} + E(\vec{R}_1,\dots,\vec{R}_M)$$ you can consider its classical counterpart $$H_n^\text{cl} = \sum_I \frac{\vec{P}_I^2}{2M_I} + \sum_{I<J} \frac{Z_IZ_J}{|\vec{R}_I-\vec{R}_J|} + E(\vec{R}_1,\dots,\vec{R}_M) = \sum_I \frac{\vec{P}_I^2}{2M_I} + U(\vec{R}_1,\dots,\vec{R}_M)$$ where $\vec{P}_I$ is the momentum of the nucleus $I$. From the laws of classical mechanics you can easily verify that the force acting on nucleus I is simply given by $$\vec{F}_I = -\nabla_{\vec{R}_I} U(\vec{R}_1,\dots,\vec{R}_M).$$ So even if you consider the nuclei as classical particles, the PES has the same meaning, i.e. is the potential felt by the nuclei. Here you can easily apply a classical picture. A transition state is simply a state along your reaction coordinate where you have a local maxima of the PES for the reduced coordinates you are considering(as the C-O bond in you example).

I never looked to this into details but I think computational programs climb the PES uphill from the initial and final state of your system (i.e. they follows -$\vec{F}_I$) in order to find the TS.

• This is the sort of answer I was looking for. Based on your explanation, I think my conceptual understanding was alright, but seeing some of the math and seeing how Born-Oppenheimer fits into things is very helpful! One thing, though; if I'm exploring a reaction mechanism by modelling species and TS's using DFT, am I really exploring the PES at all, since I am just looking at the system as a sort of 'fluid of electrons', rather than solving nuclear wavefunctions? May 28 '16 at 16:13
• @gannex Maybe the fact that I considered the nuclei as quantum particles was a bit confusing... I added a section to my answer to (hopefully) clear this confusion.
– user23061
May 28 '16 at 17:12
• @gennax Basically you use DFT to solve the electronic problem and find E and then you add to E the classical nucleus-nucleus repulsion in order to obtain the total PES (which I called U).
– user23061
May 28 '16 at 17:44
• It can also be noted that often the coordinates used for a system are determined by using principle component analysis and making these reaction coordinates simply the coordinate transformations which are most likely to occur. May 28 '16 at 18:30
• Just to add that for TS search usually some eigenvector following based algorithm is used. May 28 '16 at 19:53

In this context, the PES is not proper of a reaction but the nuclei involved.

You should take in mind that, in absence of external fields (and taking apart spin), it is just a functional of the nuclei positions, and as there are many ways to describe nuclei positions there are many ways of visualize the PES.

For a "large" number of nuclei, the dimension of the PES will be very large. Ideally, 3$N$-6 with the considerations above.

How it looks like? It is hard to imagine and represent in an 2D screen. By the way, I never found it useful, no more than to imagine that the reaction is a ball rolling down a hill.

To turn it into something that can be useful, the PES is projected into a subspace spanned by the vectors linked to the coordinates of interest for the case (for example, two bond distances).

Edit:

- Why make a distinction between a 'reaction' and the movement of the system's nuclei

It could be done, but with "the PES is not proper of a reaction but the nuclei involved" I meant that the whole PES is independent of the trajectory followed by the nuclei during the reaction.

- Can you show me an example of a projection of a PES? (Or do you just mean the PEP, as it is shown in most computational journals?)

I mean the projections of the PES usually plotted. With PES is projected into a subspace spanned by the vectors linked to the coordinates of interest for the case I meant the following. The whole functional cannot be plotted in high dimensional cases. So, the trick is chose a set of variables that allow the description of the nuclei positions BUT also allow just use a few to describe the important part (which depends on your needs).

A simple example would be the PES corresponding to the model of a benzene molecule (built for teaching purposes) in a laboratory when you are interested in teaching parabolic movement. Of course you can describe it with Cartesian coordinates of each nuclei in R^3 (12*3=36 coordinates). But they are not very useful, you need to make appear the $x^2$ somehow for your young students. So you make a change of basis and get a new set of coordinates. Three translational of the geometrical centre and internal coordinates. Then, you can make a matrix representation of this vector: [int1, int2, ..., int33, z, x, y]$^T$. It is horrible to write here (you'll see why soon) so fortunately found a little ball. And you reduce it to [z, x, y]$^T$ But your student have to plot the famous parabolic case in a paper, and you only need $x$ and $y$, so you can just discard $z$ or make a projection by applying the matrix representation of an operator $P= P^2$.

$$\begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \times \left[ \begin{array}{c} z \\ x \\ y \end{array} \right] = \left[ \begin{array}{c} 0 \\ x \\ y \end{array} \right]$$

now is simpler. The change is subtle (maybe pedantic), one case is an approximated plot of something exact, the second one is an exact plot of something approximated. By the way, the plot is the same...

My advise is: Just look at the PES as a simple function $R^N \rightarrow R$ with $N$ at least as 3$n$-6. There is nothing special about it that makes it different from what you studied in your algebra/calculus classes. The only one particularity is what the variables represents (nuclear positions), it was well described by the @R.M. answer.

I really never find it too useful more than a fast inspection of simple energy curve. I found the mathematical machinery of the optimisation algorithms much more useful. But of course, I looks very nice in papers, I would recommend a good colour scheme.

• It has sort of helped. Why make a distinction between a 'reaction' and the movement of the system's nuclei. That seems a bit pedantic to me, but I get what you mean. Can you show me an example of a projection of a PES? (Or do you just mean the PEP, as it is shown in most computational journals?) What I really want to do is get a better understanding of how I am interacting with the PES of a system when I am looking for certain minima/maxima in my calculations. I'm just looking for a better qualitative understanding of what I am doing when I apply this process. May 28 '16 at 12:59
• Your welcome @gannex May 28 '16 at 17:02
• Why down-voted ? May 28 '16 at 17:02
• I didn't downvote, but why you say that the PES is an $R^N\mapsto R^M$ function?
– user23061
May 28 '16 at 19:08
• @R.M. Because I was still thinking in the the PES simplification. I should put $R^N \rightarrow R$. Thank you for notice the mistake. May 28 '16 at 19:47