Looking at local classical and quantum forces in stable crystals using multipole-model refined electron densities and orbital-free DFT approximations

Two distinct approaches, that of energy and that of force, are adopted in quantum mechanics to get insights on chemical processes. In the second one, the net forces acting on the electrons and nuclei in a system (Ehrnefest and Hellmann–Feynman forces, respectively) are determined and a local version of the approach, in terms of force density fields rather than forces, has also been proposed for electrons. This is the path followed by Tsirelson & Stash (2020) in this issue of Acta Crystallographica Section B, to study for the first time the spatial distribution of the electronic forces of different nature acting in stable crystals. Interestingly, by relying on approximations taken from orbital-free DFT, all components of the inner-crystal force can be easily retrieved from multipole-model refined experimental electron densities and their derivatives. No less important is that these calculations are becoming easily doable for any X-ray density crystallographer thanks to a new version of the computer program WinXPRO, purposely developed in the study which is discussed in this commentary.

In their article, Tsirelson and Stash (T&S) describe a method to obtain a local mapping of electrostatic and quantum forces acting on electrons in stable crystals. To achieve their goal, T&S have carefully selected and ingeniously combined several ingredients that might not be immediately familiar to crystallographers. The ingredients include the quantum mechanical force approach in chemistry (Deb, 1973), the theory of gradient dynamical systems (Abraham & Shaw, 1992), Orbital-Free Density Functional Theory (Wesolowski & Wang, 2013) and the old, but recently broadened and revitalized field of Quantum Crystallography (Genoni et al., 2018). Each of these bits and pieces of science deserves a brief introduction on its own, before analysing the role they play when combined together (Fig. 1).

Figure 1 Ingredients and concepts from several disciplines are ingeniously combined in a Quantum Crystallography approach to analyse the balance of local forces acting on electrons in stable crystals (see text). Forces are derived as gradients of their associated classical and orbital-free DFT quantum mechanical potentials, evaluated using experimental electron densities and derivatives from high-resolution X-ray diffraction experiments. The central plot [taken from Fig. 5(b) in Tsirelson & Stash (2020)] shows the kinetic potential iso-contour lines and the kinetic force field arrows in the (110) plane of diamond. The space of the crystal turns out to be partitioned in atomic-like basins.

In quantum chemistry, rather than using the more customary energy approach, chemical bonding and molecular processes may be studied using the force approach, that is in terms of the forces of various nature acting on the electrons and nuclei in a system. In the local force form adopted by T&S, forces acting on electrons are expressed in terms of force density fields (Deb, 1979; Bader et al., 2007; Pendás & Hernández-Trujillo, 2012; Pendás et al., 2016). The latter may be obtained either from the expectation values of the gradients of the corresponding potential operators or, as performed by T&S, as gradients of the associated scalar potential fields.

These scalar potentials and their gradient fields (the forces) are analysed by T&S in terms of their topology via the theory of gradient dynamical systems, using its typical entities (attractors, repellers, basins, gradient paths, critical points) and concepts (homeomorphism). It is a tool common to most of the Quantum Chemical Topology methods, such as Bader's Quantum Theory of Atoms in Molecules (QTAIM) (Bader, 1990). Attractors and repellers are those local potential maxima and local potential minima where gradient paths (the force lines) converge to and are expelled from, leading to a partitioning of chemical systems into space-exhaustive pieces, named basins. Each basin is traversed by all ascending gradient paths (GPs) terminating at its own 3D attractor (an alternative partitioning in terms of descending GPs and 3D repellers is equally possible). In most cases the basins bear a chemical meaning. Those with 3D attractors are often atom-like in shape and in QTAIM they are identified as the (quantum) atoms in molecules. The critical points (CPs) are those points where the gradient of the potential, the force, vanishes and are characterized by their signature. The signature is given by the dimensionality of the space traversed by the ascending GPs expelled from the CP minus that of the space traversed by the ascending GPs with terminus at the CP. Signatures −3 and +3 correspond to attractors (local maxima) and repellers (local minima), while the −1 and +1 signatures denote first- and second-order saddle points, respectively. If two potential fields exhibit the same number of CPs for each signature, the fields are said to be homeo­morphic within the system where this occurs. Basins sharing a portion of their bounding surface imply a CP of −1 signature lying on their common surface and being the potential maximum for it. Two unique ascending GPs originate at this CP. One GP ends up at the attractor of one of the two basins sharing the surface and the other at the attractor of the remaining basin. The union of such paths is a line joining these two attractors and with maximum potential value relative to any lateral displacement from the line. In QTAIM it is known as a bond path, a line of maximum electron density ρ linking the nuclei of a pair of atoms. Accumulation of ρ along the bond path has been associated with chemical bonding (Bader, 1990). T&S discuss the analogy and dissimilarity of the topology of the scalar potential and force fields they analyse with that typical of QTAIM (where ρ is taken as the potential and its gradient ∇ρ as the corresponding force).

Orbital-free density functional theory (OF-DFT) and the methods derived thereof represent a fascinating challenge to go beyond the established routes in computational DFT. In principle, Hohenberg–Kohn theorems make it possible to describe the ground state of a multi-electron system exactly without constructing the wavefunction explicitly. Yet, the extreme difficulty of building sufficiently accurate approximations to the electron density (ED) functional for the kinetic energy of a many-electron system progressively led to the adoption of orbital-based formal frameworks in practical DFT (the Kohn–Sham DFT method being the most well known solution). For this reason and despite its conceptual superiority, OF-DFT was put aside for a long time since in its early Thomas–Fermi implementation all stable molecules were found to dissociate. However, in the last two decades the situation has dramatically changed and research on OF-DFT revitalized. Since their computational cost grows quasi-linearly with the size of the system, OF-DFT methods have a great appeal for treating extremely large systems, while the conventional orbital-based methods have a well known highly non-linear dependence. This pleasing OF-DFT feature has fuelled the studies on the exact properties of the kinetic energy density functional and the hunt for its increasingly better approximations, making current OF-DFT computations more reliable by far (Wesolowski & Wang, 2013). Besides, the kinetic energy and the electrostatic, exchange and correlation energy densities, all entering with their corresponding functionals in OF-DFT approaches, were demonstrated to be of great importance for interpreting key chemical concepts (Bader, 1990; Gatti, 2005; Wesolowski & Wang, 2013; Pendás et al., 2016). Using information from experimental EDs and their derivatives, Tsirelson (2007) has in the past pioneered the study of bonding in crystals in terms of these various local OF-DFT energies. T&S now take a significant step forward. Rather than the energies, they study the OF-DFT potentials and forces derived thereof, using again experimentally derived densities and examining which forces behave as heterotropic and which as homotropic (Revici, 1961). The latter tend to preserve homogenous disorganization of the matter, while the former favour its local organization, such as the electronic shell structure around nuclei.

Next to `orbital-free', the little known discipline of Quantum Crystallography (QCr) appears into the title of T&S's paper before the aim of `viewing forces in crystals' is declared. QCr was originally defined as the operation of using `crystallographic information to enhance quantum mechanical calculations and the information derived from them' (Massa et al., 1995) or, conversely, of that of performing `ad hoc quantum mechanical calculations to enhance the accuracy of the crystallographic refinement' (Huang et al., 1999). None of these definitions strictly applies to the work performed by T&S, yet it surely does to the more recent, highly broadened QCr view (Genoni et al., 2018), embracing many other active and emerging research areas involving quantum mechanics and scattering experiments. Significantly, in 2017 the IUCr Commission on Charge and Spin Momentum Densities was renamed as the Commission on Quantum Crystallography. Orbital-free QCr can then be seen as a successful attempt of using precious information obtained from crystallographic methods for evaluating potential and forces in crystals within a fundamental quantum mechanical approach.

After illustrating their multi-faceted scientific framework, the results obtained by T&S deserve a brief overview and concluding comments. Within DFT, and considering a number of simplifications (Tsirelson & Stash, 2020), the local one-electron potential μ(r) is given by the sum of a kinetic potential, ν_kin(r), and of the so-called potential acting on an electron in a molecule (PAEM; Zhao & Yang, 2014), ν_PAEM(r). For a stable system, μ(r) vanishes everywhere. Thereby, the two electron potentials balance themselves everywhere, ν_kin(r) = −ν_PAEM(r), and so do their corresponding forces −f_kin(r) = f_PAEM(r) [with ∇ν_kin(r) = −f_kin(r) and −∇ν_PAEM(r) = f_PAEM(r)]. Since ν_kin and ν_PAEM differ in sign, they have the same total number of CPs and placed at the same locations, but each with reversed signature. For instance, the local maxima of one potential become the local minima of the other and vice versa. The kinetic potential includes two terms, one due to the electron fluctuations arising from the uncertainty principle and the other to the corrections arising from the wavefunction antisymmetry requirements. Both terms are clearly quantum mechanical in nature. Instead, PAEM is composed of the electrostatic, the exchange and the correlation potential contributions, the first one being purely classical.

T&S have analysed all ν_kin and ν_PAEM components for two simple stable crystals, diamond and rock salt, using ED and its derivatives obtained from high-resolution X-ray diffraction experiments, followed by an accurate multipole-model refinement of the electronic structure factors. Many are the intriguing results of their study. First, it is found that the forces of kinetic origin draw electrons from the nuclei and are heterotropic in nature, while all f_PAEM components behave as homotropic. The latter push electrons to the nuclei that therefore act as sinks for the PAEM field, rather than as local maxima for ν_kin. Hence, only the kinetic potential and its associated forces are able to create the electronic shells yielding atomic-like structures and basins resembling those of ρ in the QTAIM, yet with noticeable differences (Tsirelson & Stash, 2020). Both ν_kin and ν_PAEM fields concur to the accumulation of electronic charge along curved lines (GPs) joining nearest atoms. However, not all atomic pairs linked by bond paths in the ρ field (QTAIM) are found to be joined by the corresponding GPs in the ν_kin field. This fact has an impact on the topologies of the two fields that are therefore not always homeomorphic. Even more importantly, it opens the question on the physical meaning of both the lack of a shared portion of bounding surface and of a common CP of signature −1 in some circumstances (in the present case, the Cl⋯Cl interaction in the NaCl crystal, which is missing in the ν_kin field and, instead, present in the ρ field). Interestingly, the missing CPs in the ν_kin field topologically correspond to the QTAIM bond critical points (bcps), whose meaning is and has been widely debated in the literature [see for instance, Pendás et al. (2007)]. Understanding the reasons behind the possible dissimilarities of the ρ and the ν_kin fields could shed further light on the physical interpretation of bcps. It will certainly be a challenge for forthcoming contributions in the area pioneered by T&S in their paper.