Bridging electron and nuclear motions in chemical reactions through electrostatic forces from reactive orbitals

Electrostatic forces exerted by electrons on nuclei

To understand the forces acting on nuclei during chemical reactions, we begin by examining the electrostatic forces exerted by all electrons on the nuclei within a molecule. The Hamiltonian, \(\hat \), governing the system of electrons and nuclei²⁴, is given by

where \({\nabla }_{i}^{2}\) denotes the Laplacian with respect to electron i, r_iA is the distance from electron i to nucleus A, r_ij is the inter-electronic distance between electrons i and j, R_AB represents the internuclear distance between nuclei A and B, Z_A is the charge of nucleus A, and n_elec and n_nuc represent the numbers of electrons and nuclei, respectively. Atomic units are used (ℏ  =  e² =  m =  1, energies are in hartree, and distances are in bohr). Building upon electrostatic theory¹⁵, the Hellmann-Feynman force exerted by the electrons and nuclei on nucleus A is expressed as

where ρ(r) represents the electron density at position r, and R_A and R_AB are position vector of nucleus A and vector from nucleus A to nucleus B, respectively. According to the Hellmann-Feynman theorem, these forces represent classical electrostatic forces when the electron distribution is determined variationally¹⁶. The approximation in the right-hand side of Eq. (2) neglects physical effects that depend on nuclear coordinates, such as dispersion forces, which scale with the inverse sixth power of the interatomic distance. It is known that while dispersion forces can affect the structures of weakly bound molecules, they have little influence on molecular orbitals. Therefore, their contribution can generally be neglected, except in cases where the Hellmann-Feynman force is extremely weak.

For a wavefunction describing all electrons, the electron density is given by

which is integrated over all spatial and spin coordinates, x_n = (r_n, s_n) for n  = 1 to N, except r₁. Within independent electron approximations²⁵, such as the Kohn-Sham method²⁶, the electron density simplifies to

where ϕ_i is the i-th spin orbital wavefunction. The total electrostatic force on nucleus A is then expressed as the sum of contributions from electrons (\({{{{{\bf{F}}}}}}_{A}^{{{{{\rm{elec}}}}}}\)) and other nuclei (\({{{{{\bf{F}}}}}}_{A}^{{{{{\rm{nuc}}}}}}\)),

where the force contribution from the i-th orbital on nucleus A is

Using this framework, the influence of reactive orbitals on nuclear forces can be isolated. The force contribution from an occupied reactive orbital (ORO) is given by

where ϕ^ORO represents the wavefunction of the ORO. This formulation allows for the direct assessment of how variations in reactive orbitals influence nuclear forces during chemical reactions.

In this study, we investigate the primary electrostatic forces driving nuclear motions in chemical reactions, focusing specifically on the forces generated by variations in the ORO. The progression of electron transfer throughout a reaction is characterized by changes in OROs¹³. The electrostatic force vector resulting from ORO variations can be expressed as:

where we refer to \({{{{{\bf{F}}}}}}_{A}^{{{{{\rm{ROEF}}}}}}\) as the reactive orbital-based electrostatic force (ROEF) vector. This formulation isolates the contribution of ORO variations along the reaction pathway, which may involve interactions with unoccupied orbitals, while neglecting the contributions from other orbitals.

Relationship between electrostatic forces and orbital energies

Next, we examine the relationship between ROEF and the Kohn-Sham orbital energies²⁴. The Kohn-Sham orbital energy, ϵ_i, is defined as

where v_xc is an exchange-correlation potential functional, and h_i represents the one-electron Hamiltonian, given by

The derivative of the orbital energy with respect to the coordinates of nucleus A is then expressed as

where contributions from the v_xc term are negligible due to its limited explicit dependence on nuclear positions.

Reactive-orbital-based electrostatic force (ROEF)

The total electrostatic force vector on nucleus A, derived from the Kohn-Sham electronic energy, E_KS, can be written as

indicating that the force contribution from each orbital is proportional to its energy gradient with respect to nuclear coordinates. Substituting this relationship, the ROEF vector can be expressed as

where ϵ^ORO is the orbital energy of the ORO. This equation directly relates the orbital energy gradients to the forces exerted on the nuclei, providing a quantitative link between electronic structure and nuclear motion. By focusing on the ORO, this approach establishes a direct connection between orbital energy variations and nuclear forces, enabling a more detailed understanding of the forces driving chemical reactions at the atomic level.

Accuracy and interpretation of ROEF

From the relationship established in Eq. (13), we infer that orbital energy variations influence electrostatic forces during chemical reactions. Specifically, for OROs that exhibit decreasing orbital energies as the reaction progresses, the resulting electrostatic forces are expected to align with the reaction direction, as described in Eq. (8). These reaction-aligned electrostatic forces create trenches on the PES, delineating the intrinsic reaction coordinate (IRC) and defining the reaction pathway. OROs act as the driving force for a reaction when their electrostatic forces are sufficiently strong to establish the IRC. Since orbital energy remains constant during idealized electron transfer²⁷, orbital energy variations are minimal in reaction stages dominated by electron transfer. This behavior is reflected in the early stages of many reactions, as shown in Figs. S1 and S2 of the Supplementary Information. According to Eq. (13), when the orbital energy gradient approaches zero, the corresponding electrostatic force vector vanishes. However, even small gradients in orbital energy, arising from structural transformations along the IRC, determine the magnitude of the electrostatic force. Accurate calculation of these forces, therefore, hinges on precise determination of orbital energy variations induced by structural changes. This underscores the necessity of a robust theoretical framework capable of reliably predicting orbital energy gradients across different molecular structures. LC-DFT provides a critical tool in this regard, as it quantitatively reproduces orbital energies with high fidelity, making it indispensable for such analyses.

The Pulay force, an artificial force arising from the use of basis sets, is a factor potentially related to ROEFs^28,29,30. Previous studies report that the magnitude of Pulay forces for all electrons is less than 10 kJ mol⁻¹ bohr⁻¹, with even smaller contributions for individual molecular orbitals. By contrast, ROEF derived from orbital energy gradients are on the order of several hundred kJ mol⁻¹ bohr⁻¹. Given the negligible magnitude of Pulay forces compared to ROEFs, they are excluded from this analysis for simplicity.

An important insight from the interplay between ROEFs and orbital energy variations relates to the virial theorem³¹, which states that the kinetic energy is half the potential energy but with the opposite sign for nuclear-electron interactions. Under independent electron approximations, this relationship holds for individual electrons³². When orbital energy variations during orbital mixing are small, nuclear-electron potential changes are similarly minimized, suppressing fluctuations in ROEFs, as shown in Eqs. (10) and (12). While electron transfer induces polarization in the electron distribution, ROEF remains stable. However, structural deformations near the transition state can generate significant ROEFs, influencing the reaction pathway.

This study highlights the critical role of electron transfer as the driving force for structural transformations during reactions. By calculating the electrostatic forces associated with ORO energy gradients along the IRC, we trace the forces exerted on atomic nuclei, revealing how changes in electron distribution induced by electron transfer shape the reaction pathway.

Three-body model for evaluating ROEFs

Using the framework of ROEFs, we analyze a diverse set of atom-transfer reactions, comprising 16 hydrogen transfer reactions, one heavy atom transfer, two nucleophilic substitutions, and five unimolecular reactions, for the forward and backward processes, based on their calculated IRCs. For the ROEF analysis, atom-transfer reactions were modeled using a three-body approach.

Figure 1 illustrates this framework, applied to the NH + CH₄ → NH₂ + CH₃ reaction (a) and the ROEF vector calculation method (b). These ROEF values provide insight into the forces driving atomic movement caused by electron redistribution during ORO variations. This three-body model provides a simple and intuitive way to visualize the correlation between the reaction direction and electrostatic forces. However, it should be noted that this model is not applicable to more complex reactions such as catalytic or enzymatic processes, which are beyond the scope of the present study. To visualize such correlations in complex reactions, a new model will be required—one that projects the electrostatic force vectors onto the nuclear displacement vectors along the reaction coordinate. ROEF variations along the IRCs for all calculated reactions, along with corresponding molecular structures and atomic coordinates, are presented in Figs. S1 and S2 of the Supplementary Information.

a Example calculation model for the NH + CH₄  → NH₂ + CH₃ reaction at the transition state, where atoms C1, N2, and H6 are labeled as A, B, and C, and the movement of atom C from A to B is considered. b Method for calculating ROEFs in atom-transfer reactions: F_A, F_B, and F_C are the ROEF vectors calculated for atoms A, B, and C, respectively, \({{{{{\bf{F}}}}}}_{{{{{\rm{A}}}}}}^{{\prime} }\) and \({{{{{\bf{F}}}}}}_{{{{{\rm{C}}}}}}^{{\prime} }\) are the projection vectors of F_A and F_C in the direction of the A-C bond, F_C” and F_B” are the projection vectors of F_C and F_B in the direction of the C-B bond. These projections enable the calculation of ROEF vectors, F_BC and F_CA, crucial to the reaction. For unimolecular reactions, only one ROEF vector set corresponds to the reaction direction, and the remaining vectors are set to zero.