Visualization of electron density changes along chemical reaction pathways

We propose a simple procedure for visualising the electron density changes (EDC) during a chemical reaction, which is based on a mapping of rectangular grid points for a stationary structure into (distorted) positions around atoms of another stationary structure. Specifically, during a small step along the minimum energy pathway (MEP), the displacement of each grid point is obtained as a linear combination of the motion of all atoms, with the contribution from each atom scaled by the corresponding Hirshfeld weight. For several reactions (identity S 2, Claisen rearrangement, Diels-Alder reaction, cycloaddition, and phenylethyl mercaptan attack on pericosine A), our EDC plots showed an expected reduction of electron densities around severed bonds (or those with the bond-order lowered), with the opposite observed for newly-formed or enhanced chemical bonds. The EDC plots were also shown for copper triflate catalyzed N O fragmentation, where the N–O bond weakening initially occurred on a singlet surface, but continued on a triplet surface after reaching the minimum-energy crossing point (MECP) between the two potential energy surfaces. GRAPHICAL ABSTRACT


Introduction
In chemical reactions, the electron density of a molecular system undergoes significant changes. For simple reactions, such as identity S N 2 and Menshutkin reactions, electron density decreases around broken bonds, whereas the density increases at newly-formed chemical bonds. For other reactions, such as Claisen rearrangement and Diels-Alder reactions, density changes are also expected around double bonds converted into single bonds (or vice versa). It would be highly desirable to be able to directly visualise such electron density changes (EDC) during a chemical reaction. However, the EDC visualization is not as straightforward as it seems. In computational chemistry, the electron density of each stationary structure (reactant, transition state, and product) of a reaction can be readily computed on a rectangular grid, and then visualised via an isosurface representation. Naively, one can compute the EDC from one stationary structure to another by subtracting the density values on the same rectangular grid. But, the corresponding isosurface plot tells us more about the displacement of each atom (which we already know) than the actual change in the density around key atoms and bonds.
In this work, we aim to develop a simple procedure for enabling the EDC visualization. Before proceeding to introduce our procedure, however, we would like to point out that numerous electron structure-based tools have been developed for analysing chemical reactions. Below we briefly mention some of these tools, which we would loosely classify into atomic-population-based, bond-order/energy-based, and topology-based ones. The readers are referred to reviews for many other available tools [1][2][3][4].
Atomic charge populations are the intuitive coarsegrained representation of the electronic density [5]. But they are not uniquely defined, practically leading to numerous procedures for computing these charges [6]. These procedures include the decomposition of the oneparticle density matrix (actually its product with the overlap matrix, leading to Mulliken and Löwdin charges [7,8]), the partitioning of the electronic density (such as the Bader Atom-in-Molecule charges [9] and the Hirshfeld and related charges [10,11]), the fitting of electrostatic potentials (such as Merz-Kollman ESP charges [12][13][14] and CHELPG charges [15]), and the natural bond orbital analysis [16,17]. Besides examining the atomic charge differences between stationary structures, one can also analyse the Fukui functions, [18][19][20] which are often approximated as the electron density difference between a neutral molecule and its anion/cation. The corresponding 'atomic charges' are widely used to identify electrophilic and nucleophilic sites of reactants. Related to these conceptual density functional theory indices, the reaction electronic flux (REF) method decomposes the chemical potential calculated along the reaction coordinate into polarisation and charge transfer components [21][22][23][24]. This allows different regions along the reaction coordinate to be characterised by their dominance by either polarisation or charge transfer effects, thus providing additional chemical insights into the reaction mechanisms.
The bond order, especially the one proposed by Mayer, [25] has been widely used to characterise covalent bonds in molecules. The corresponding bond energy, which can be obtained using various schemes to decompose the total energy of a molecule, [26][27][28][29] can in principle also be employed to assess the strength change of chemical bonds in a reaction. Separately, the bond order can also be gauged via the local vibrational analysis [30], and a chemical reaction can be analysed using the related unified reaction valley approach (URVA) [31][32][33].
Topological analyses of electron density, density gradient, and wavefunction based quantities [such as the electron localisation function (ELF), [34,35] non-covalent interactions (NCI) analysis, [36] and dynamic Voronoi metropolis sampling (DVMS) [37,38]] have been proposed by several research groups. These tools originated from the idea of partitioning the molecular space into different regions, as in the aforementioned Bader's Atoms-in-Molecule (AIM) theory [9]. The ELF function from Becke and Edgecombe, [34] for instance, can be used to divide the electron density into the core basins connecting two nuclei and valence basins representing non-bonding regions and lone pairs [35] and to study the bond evolution in chemical reactions [39,40]. Similarly, the NCI plot can be used to identify the key noncovalent interaction regions in molecules at the equilibrium geometry and transition states [36]. More recently, the DVMS method [37,38] utilised a metropolis walk in the electron configuration space to acquire wavefunction tiles -signed regions of the wavefunction-that can be followed to visualise the motion of the electrons along a reaction pathway.
Going back to the EDC, which is the main focus of this manuscript, we will propose a simple method for directly visualizing EDC (without using wavefunction tiles) during a chemical reaction, which will complement the aforementioned analysis tools. The key step for our EDC visualization, as outlined in Section 2, is to map a rectangular grid for the first stationary structure (such as reactant) onto the space of the second stationary structure (such as transition state). Through evaluating the electron density on the distorted grid for the second structure, we can assess/visualise the changes to the electron density on the original rectangular grid for the first structure. Such EDCs will be shown in Section 4 for several model reactions, which conform well with our chemical intuition.
We would like to dedicate this article to Prof. Peter Gill, whose elegant and masterly applications of simple and advanced mathematical tools in the field of quantum chemistry have been a great inspiration to Y.S. and other authors. To a large extent, our work is an attempt to follow Peter's example and develop simple mathematical tools to advance the chemical understanding. However, as will become obvious, our grid mapping has a serious shortcoming in that it does not conserve the number of electrons. Hopefully, more rigorous procedures will be developed in the future by readers of this special issue.

Method
Let us consider two stationary structures, 0 and 1. We will denote the corresponding atomic coordinates as R (0) A and R (1) A . The key step is to map a rectangular grid around the first structure, G In other words, each grid point in G k , which results from the molecular structure change.
Let us use the transformation of three acetylene molecules to benzene as an example to illustrate what a mapped grid should look like. The top panel of Figure 1 shows a rectangular grid (with a 0.1 Å grid size) around three acetylene molecules arranged in a hexagonal complex. Clearly, there are fewer grid points between the two carbon atoms within each acetylene molecule (separated by 1.215 Å) than between the closest carbon atoms from neighbouring molecules (distance: 2.410 Å). As the distance between the acetylene molecules shortens to 1.398 Å (which is the equilibrium C-C bond length of benzene), we expect the grid points between these molecules to be squeezed closer to each other, as shown in the bottom panel of Figure 1. Meanwhile, the C-C bond within each acetylene molecule gets stretched, so one would expect a larger spacing between the grid points between the two carbon atoms of each acetylene molecule.
There might be different ways for carrying out such mapping of a rectangular grid into a distorted one. In this work, we will utilise a continuous reaction pathway R A (s) that connects the two stationary structures, R (0) A (s = 0) and R (1) A (s = 1). This allows us to implement the mapping in Equation (1) as where w A,k (s) refer to a set of atomic weights for each grid point G k (s) at a given structure along the reaction pathway. For a discretised pathway, this can be carried out with which is illustrated in Figure 2. Furthermore, in this work, we will adopt the Hirshfeld function for the atomic weights [10] In this expression, the denominator is the so-called promolecule density, which is a sum of atomic densities ρ (free) B (r) from each free atom. After the mapping, the density changes to each rectangular grid point for structure 0 can be computed as

Computational and implementation details
In this work, all the stationary structures (reactants, transition states, and products), the minimum energy crossing point (MECP), and the reactant pathways are optimised using the density functional theory (DFT) [41][42][43]. Specifically, the ωB97X-D functional [44] and the 6 − 31 + G * basis set [45,46], as implemented in Q-Chem version 5.2, [47] were employed. The transition state (TS) structures were obtained in three steps: (a) build a freezing string between the reactant and product geometries, [48,49] (b) take the highest-energy point along the string, and (c) refine the structure using Baker's partitioned rational-function optimisation (P-RFO) method [50]. Subsequent harmonic vibrational frequency analysis was performed to confirm a local minima (no imaginary frequency) for the reactant and product or a saddle point (one imaginary frequency) for the transition state. Intrinsic reaction coordinate (IRC) calculations were carried out to connect the TS to the reactant and product structures. The MECP structure for the spin-crossing reaction was optimised using the algorithm proposed by Harvey and coworkers [51] as implemented in the sobMECP software by Lu (http://sobereva.com/286). The EDC analysis was carried out using QMAnalysis, a locally developed Python program. Each EDC plot was created using a 0.1 Å grid, while each movie provided in the supplemental information was created with a 0.2 Å grid. For each analysis, it only requires: (a) the Q-Chem checkpoint files for the two end structures (reactant/TS, TS/product, reactant/MECP, or MECP/product); and (b) the pathway file (in xyz format) in between. The QMAnalysis program can be accessed at https://github.com/cc-ats/edc, which is being further developed. The checkpoint and pathway files can also be found in this GitHub repository.

Results
To test our EDC analysis tool, we considered six representative unimolecular and bimolecular reactions: (A) the identity S N 2 reaction; (B) the Claisen rearrangement of allyl vinyl aether; (C) the Diels-Alder cycloaddition between cyclopentadiene and acrylonitrile; (D) the [3 + 2] cycloaddition of the azomethine ylide pseudoradical with ethylene; (E) the S N 2 reaction between pericosine A and phenylethyl mercaptane; and (F) the splitting of N 2 O bound to copper triflate. When multiple bond forming and breaking are involved, some of the above reactions are asynchronous (such as the example Claisen rearrangement, Diels-Alder, and pericosine reactions), whereas others are concerted (such as the identity S N 2 and [3 + 2] cycloaddition). While most of the reactions are spin-conserving, the last reaction involves a crossing from a singlet surface of the N 2 O-Cu-Triflate complex to the triplet one. Besides the EDC plots shown in the following main text, the EDC movies of the studied reactions and the reaction energy profiles can be found in the supporting information (SI).

Identity S N 2 reaction
Bimolecular nucleophilic substitution S N 2 reactions are widely studied both experimentally [52,53] and theoretically [54][55][56]. Here, we shall use the identity S N 2 reaction between methyl chloride and chloride anion to test our new EDC method. The reactant structure (Figure 3(a)) was constructed from two fragments, namely chloride anion (Cl − A ) and methyl chloride (CH 3 Cl B ), which were separated by a Cl A -C distance of 3.22 Å. The nucleophilic attack of the chloride anion led to the transition state in Figure 3(b). At the transition state, the carbon atom of the methyl chloride is equidistant from two chlorine atoms with both C-Cl A and C-Cl B distances being 2.34 Å. From the TS structure, the product was generated with the C-Cl A bond distance of 1.82 Å, and chloride anion (Cl − B ) as the leaving group. Figure 3(c) showed the electron density change from the reactant to the transition state. As expected, it showed a decrease in the electron density from the C-Cl B σ bond (shown in magenta) of methyl chloride and the emergence of a new bond (shown in blue) between Cl A and the carbon atom of the methyl chloride. In this stage, the three C-H bonds of methyl chloride also gained a minor increase in their electron densities. A complete formation of the Cl A -C bond (with a large blue cloud) and the further weakening of the C-Cl B bond could be seen during the subsequent transition from TS to the product in Figure 3(d). Movies showing the EDC for entire pathway can be found in the supporting information for all reactions studied in our work, with sn2_iso02.mp4 for the identity S N 2 reaction.

Claisen rearrangement
The [3,3] sigmatropic rearrangement of the allyl vinyl aethers to the γ ,δ-unsaturated carbonyl compounds, [57] also known as the Claisen rearrangement or analogous Cope rearrangement, [58,59] were shown in Figure 4. In the reactant structure, the aether bond C 4 -O 3 has a normal bond length of 1.42 Å. As the reaction progressed toward the TS, one observed a lengthening of the C 4 -O 3 aether bond from 1.42 Å to 1.94 Å. This was accompanied by a conformational change, namely a change of the O 3 -C 4 -C 5 -C 6 dihedral angle, which reduced the C 1 -C 6 distance (between two terminal allyl carbon atoms) from 3.37 to 2.33 Å and resulted in a chair-like TS (Figure 4(e)). In the SI Figures S1 and S2, the chair form of the TS was shown to be 6.0 kcal/mol more stable than the boat form, consistent with previous computational analysis [60].
The EDC plots in Figure 4(d,e) showed that the weakening of the C 4 -O 3 σ bond (shown in magenta) and the strengthening of the C 2 -O 3 bond occurred in both reactant→TS and TS→product stages. On the other hand, the formation of the C 1 -C 6 σ bond and the weakening of C 1 -C 2 and C 5 -C 6 bonds occurred mostly after the TS.

Diels-Alder cycloaddition
The [4 + 2] Diels-Alder cycloaddition between cyclopentadiene and acrylonitrile is well-known for its asynchronous reaction pathway [61][62][63][64][65]. In the reactant   Figure 5(a)), the cyclopentadiene and acrylonitrile molecules are separated by 3.17-3.24 Å. As the reaction progressed, the two fragments moved closer to each other. At the TS structure ( Figure 5(b)), the distance between the C 4 carbon of cyclopentadiene and the C 6 carbon of acrylonitrile was reduced to be 2.09 Å. Meanwhile, the C 1 -C 7 distance was reduced only to 2.43 Å, suggesting an asynchronous cycloaddition caused by the presence of the cyano group during the distance reduction between the double bond and (gradually nonplanar) cyclopentadiene. After the transition state, the formation of C 4 -C 6 bond and C 1 -C 7 bonds were both completed, with the final bond distances being 1.56 Å.

structure (
As shown in Figure 5(d), the reaction started with a transfer of electrons of three C-C double bonds of acrylonitrile and cyclopentadiene (shown in magenta) to the formation of the C 4 -C 6 and C 1 -C 7 bonds (shown in blue), respectively. Noticeably, a more substantial electron density change occurred between the C 4 and C 6 atoms, consistent with a short distance at the transition state. At the same time, a relatively smaller electron density change was observed between C 1 -C 7 atoms, indicating a slower bond formation. In addition to the partial formation of these new σ bonds, the C 2 -C 3 bond of the cyclopentadiene also gained electron density, bringing it a partial double-bond character. The EDC from the TS to the product showed the completion of the reaction ( Figure 5(e)), where the π bond of acrylonitrile (C 6 -C 7 ) is the primary site for losing electrons while C 1 -C 2 and C 3 -C 4 bonds of cyclopentadiene continued to lose electrons (in magenta). On the other hand, the electron density on the cyano group underwent only marginal changes throughout the reaction.

Cycloaddition between azomethine ylide and ethylene
The pseudoradical (pr)-type azomethine ylide and ethylene can undergo the [3 + 2] cycloaddition to form regioisomeric pyrrolidines [66,67]. Domingo and coworkers suggested that the pseudoradical nature of azomethine ylide is the primary reason for the low activation energy barrier (1.2 kcal/mol) of the above cycloadditon reaction [68,69]. Here, the reactant structure consisted of two fragments: the azomethine ylide and ethylene ( Figure 6(a)). The initial distances are 3.10 Å for both C 1 -C 5 and C 3 -C 4 . During the reaction, these fragments came closer, reducing the C 1 -C 5 and C 3 -C 4 distances from 3.10 Å to 2.49 Å at TS. Subsequently, during the transition from the TS to the product, both bond lengths were further decreased to 1.55 Å. In contrast to the Diels-Alder reaction in Section 4.3, this [3 + 2] cycloaddition reaction proceeded in a concerted fashion.
The EDC plots in panels (d) and (e) of Figure 6 confirmed the concerted mechanism of this cycloaddition reaction. During the transition from the reactant to TS, the π bond of ethylene and π orbitals of azomethine ylide lost some electron density (shown in magenta). At the same time, two nascent bonds with an equal amount of electron density could be seen (marked blue) between two fragments. During the subsequent transformation from the TS to product, two single bonds between ethylene and azomethine ylide fragments were completely formed (shown in blue), whereas the π orbitals of both ethylene and azomethine ylide orbitals were significant donors of the electron density (shown in magenta).

S N 2 reaction between pericosine A and phenylethyl mercaptan
To neutralise the pungent odour of skunk spray, which arises mainly from mercaptans, Du and coworkers reported a chemoreactive fungal metabolite of Tolypocladium sp. termed pericosine A [70,71]. With the help of a Michael acceptor -methoxycarbonyl group -on C 4 and chlorine atom on C 6 position, the pericosine A scaffold can act as an electron sink during its reaction with the nucleophilic mercaptan. The reaction between phenylethyl mercaptan and pericosine A (shown in Figure 7) was reported to follow two potential mechanisms: syn-S N 2 and anti-S N 2 , [70] depending on the stereochemistry of the attack of thiolate on the pericosine molecule. In our calculations, we focused on the syn-S N 2 reaction pathway shown in Figure S6(a). At the TS of the syn-S N 2 reaction pathway, the S-C 5 distance was found to be 2.45 Å, and the C 1 -Cl bond is 1.91 Å long. During the transition from TS to the products, the S-C 5 bond length decreased from 2.45 to 1.86 Å, whereas the C 1 -Cl bond length increased from 1.91 to 2.92 Å.
The EDC plot in Figure 7 showed the main density changes for the syn-S N 2 reaction mechanism. During the transition from the reactant to TS, the C 5 -C 6 π bond lost some electron density (shown in magenta in Figure 7(d)), and the new S-C 5 bond was partially formed (shown in blue). From the TS to product, the S-C 5 bond was completely formed, along with a π bond migration from C 5 -C 6 to C 6 -C 1 in Figure 7(e). Notably, the C 1 -Cl bond was broken (shown in magenta) only after the TS, leading to the release of a chloride ion. This is also consistent with the change of bond lengths shown in Figure S6

Spin-crossing splitting of N 2 O bound to copper triflate
A reaction involving the splitting of N 2 O through the cleavage of the N-O bond with a copper(I) triflate has been recently studied using DFT calculations [72]. As shown in Figure 8, the reactant structure was a singlet state with the copper triflate bound to N 2 O. However, the dissociation of N 2 O resulted in the production of singlet N 2 and a triplet copper-oxo species. Thus, the overall reaction exhibited a spin-crossing transition from the singlet reactant state to the triplet product state. This resulted in an MECP as the highest energy point along the reaction path shown in Figure S7. At MECP, the N-O bond stretched to 1.79 from 1.21 Å, and the N 2 O molecule adopted a bent conformation (Figure 8(b,e)). After the MECP, the reaction proceeded on the triplet surface. When the reaction moved to the product, the N-O bond was fully broken, and the Cu atom was predicted to transition from Cu(I) to Cu(II) oxidation state.
The EDC plot in Figure 8(d) showed a cleavage of the N-O bond and strengthening of the Cu-O bond as the reaction moved towards the MECP. After the MECP, the plot in Figure 8(e) showed a further decrease in the electron density around the N-O bond and a small increase in the electron density of the N-N bond. Furthermore, the EDC around the copper atom reflected changes in the electron density of d orbitals as the copper transitions between oxidation states.

Conclusions
In this work, we proposed a simple scheme for mapping grid points from one molecular structure to another. Such a mapping allows for a direct visualization of electron density changes (EDC) along the reaction pathway. For six representative reactions, the acquired EDC plots/movies provided intuitive pictures of the (concerted or asynchronous) motion of electrons during the reaction process. Our EDC approach thus provides a handy tool for chemical reaction analysis, which complements many other existing ones in the toolbox (including the recent DVMS method).
While this work was entirely focussed on the electron density, we expect that our grid mapping is applicable to the study of changes to other scalar or vector quantities from quantum chemistry calculations, such as the energy densities [73,74] and the electrostatic potential/field, which can potentially enhance our understanding of reactive and binding processes as well as conformational changes.
We note, however, our grid mapping does not conserve the number of electrons. This is especially obvious in the TOC figure, which shows the EDC during the transformation from three acetylene molecules into benzene. There, the net increase in the electron density between acetylene molecules (marked blue) appears to be more substantial than the decrease of the electron density (shown magenta) along the original C-C triple bond.
We believe that, to overcome this, different weights need to be assigned to the points of the original rectangular grid as well as the distorted grid. One can potentially use an auxiliary grid, such as the atom-centred Lebedev grid, [75] and project the corresponding weights onto the rectangular and distorted grids. However, it remains to be tested.