# Chemical reactivity from linear response eigenfunctions and eigenvalues

Rémi Grincourt, Guillaume Hoffmann, Frédéric Guégan, Vincent Tognetti, Laurent Joubert, Henry Chermette, Alejandro Toro-Labbé, Christophe Morell

► To cite this version:

Rémi Grincourt, Guillaume Hoffmann, Frédéric Guégan, Vincent Tognetti, Laurent Joubert, et al.  
Chemical reactivity from linear response eigenfunctions and eigenvalues. The Journal of Chemical  
Physics, 2025, 163 (13), pp.134118. 10.1063/5.0277034. hal-05300071

HAL Id: hal-05300071

https://hal.science/hal-05300071v1

Submitted on 6 Oct 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
AUTHOR QUERY FORM

<figure><img src="image_2.png" /></figure>

AIP
Publishing

Journal: J. Chem. Phys.

Article Number: JCP25-AR-01519

Please provide your responses and any corrections by annotating this
PDF and uploading it to AIP's eProof website as detailed in the
Welcome email.

Dear Author,

Below are the queries associated with your article. Please answer all of these queries before sending the proof back to AIP.

**Article checklist:** In order to ensure greater accuracy, please check the following and make all necessary corrections before returning your proof.

1. Is the title of your article accurate and spelled correctly?

2. Please check affiliations including spelling, completeness, and correct linking to authors.

3. Did you remember to include acknowledgment of funding, if required, and is it accurate?

<table><thead><tr><th>Location in article</th><th>Query/Remark: click on the Q link to navigate to the appropriate spot in the proof. There, insert your comments as a PDF annotation.</th></tr></thead><tbody><tr><td>Q1</td><td>Please check that the author names are in the proper order and spelled correctly. Also, please ensure that each author's given and surnames have been correctly identified (given names are highlighted in red and surnames appear in blue).</td></tr><tr><td>Q2</td><td>For the benefit of the reader, please insert a brief explanation of the meaning of the boldface values that appear in Tables III-VII and X in the caption. For example "Boldface denotes..."</td></tr><tr><td>Q3</td><td>References 19 and 35 contain identical information. Please check and provide the correct reference or delete the duplicate reference. If the duplicate is deleted, renumber the reference list as needed and update all citations in the text.</td></tr><tr><td>Q4</td><td>Please confirm the content in Refs. 29, 33, 34, and 57, as we have inserted the required information.</td></tr><tr><td>Q5</td><td>Please confirm ORCIDs are accurate. If you wish to add an ORCID for any author that does not have one, you may do so now. For more information on ORCID, see https://orcid.org/.<br><br>Rémi Grincourt – 0009-0009-7799-9696<br>Guillaume Hoffmann – 0000-0003-1851-5451<br>Frédéric Guégan – 0000-0002-4932-8643<br>Vincent Tognetti – 0000-0002-0292-3509<br>Laurent Joubert – 0000-0003-2649-112X<br>Henry Chermette –<br>Alejandro Toro-Labbé – 0000-0001-9906-2153<br>Christophe Morell – 0000-0002-6321-8723</td></tr><tr><td>Q6</td><td>Please check and confirm the Funder(s) and Grant Reference Number(s) provided with your submission:<br>Agence Nationale de la Recherche, Award/Contract Number 23-CE29-0023, (ANR-11-LABX-0029, ANR-18- EURE-0020<br>XL CHEM<br>Conseil National des Universités<br>Please add any additional funding sources not stated above.</td></tr></tbody></table>

Thank you for your assistance.
# Chemical reactivity from linear response eigenfunctions and eigenvalues

Cite as: J. Chem. Phys. 163, 000000 (2025); doi: 10.1063/5.0277034

Submitted: 22 April 2025 • Accepted: 21 August 2025 •

Published Online: 9 99 9999

<figure><img src="image_3.png" /></figure>

View Online

<figure><img src="image_4.png" /></figure>

Export Citation

<figure><img src="image_5.png" /></figure>

CrossMark

Rémi Crincourt,<fnref n="1" /> ID
Guillaume Hoffmann,<fnref n="1" /> ID
Frédéric Cuégan,<fnref n="2" /> ID
Vincent Tognetti,<fnref n="2" /> ID
Laurent Joubert,<fnref n="3" /> ID
Henry Chermette,<fnref n="1" />
Alejandro Toro-Labbé,<fnref n="4" /> ID
and Christophe Morell<fnref n="1" />,ª ID

## AFFILIATIONS

## ABSTRACT

The atom-condensed linear response function matrix can be diagonalized, and the associated eigenvectors thus obtained form a complete basis set on which any density deformation can be projected. In this paper, one proposes to interpret these vectors as atom condensed Electron Density Deformation Modes. Each eigenvector represents a mode describing how electron density can deformed in response to an external potential perturbation. It is hypothesized that these modes can reveal the preferred direction or pattern of electron flow, helping identifying the reactive region of a molecule. A relationship between electron density polarization energy and hardness variation is also derived. Practical applications on organic reactions are provided.

Published under an exclusive license by AIP Publishing. https://doi.org/10.1063/5.0277034

## I. INTRODUCTION

Conceptual density functional theory (C-DFT) provides a physical framework to investigate the chemical reactivity and selectivity of chemical processes.<fnref n="1" />⁻<fnref n="3" /> Since the dawn of this theory, lots of concepts already developed by experimental chemists have found a sound mathematical definition and a chemical meaning. Electronegativity,<fnref n="4" /> chemical hardness,⁵ softness,⁶ and Fukui functions⁷ are prototypical examples of such rationalization. In addition, reactivity principles, such as Sanderson's principle of electronegativity equalization,⁸ have been proven within this framework. Other principles, for instance, the Pearson Maximum Hardness principle,⁹ have also emerged over the years.

Since the 1960s, organic chemistry textbooks¹⁰ have often rationalized chemical reactivity and selectivity through the famous Klopman-Salem model.¹¹,¹² In this model, the driving force of a chemical reaction is either ascribed to electrostatic control or frontier molecular orbitals (FMOs) overlap. Electron density and Fukui functions¹³⁻¹⁵ are the corresponding descriptors allowing local characterization of the molecule selectivity. During the late 1990s, Parr

and Yang¹⁶ and Chattaraj¹⁷ revisited the Klopman-Salem model within the framework of C-DFT. They showed that another driving force should also be included in the model, namely, electron density polarization.

Until the past decade, the linear response function (LRF, also named the linear response kernel), the fundamental descriptor of electron density polarization, was only barely investigated. This was surely due to its mathematical complexity. Indeed, for an electron system, LRF is defined as the derivative of the electron density at point $\vec{r}$ when the external $v$ potential is changed at point $\vec{r}'$. This results in a six-dimension function due to the two sets of spatial coordinates, which cannot be visualized without projection in a three-dimension space, an approach that was adopted, for instance, by Mussard and Ángyán¹⁸ and the Brussels group.¹⁹⁻²²

Seminal studies were also reported by Sebastiani and coworkers both in the static²³,²⁴ and frequency-dependent²⁵ regimes, but they might seem to have been a little bit overlooked within the chemistry community. Recently, Senet first,²⁶ then Fias, Geerlings, and De Proft,¹⁹⁻²² have proposed a method to work around this issue. Instead of looking at LRF wholly, they provide a way

<footnote n="1">¹ Université Claude Bernard Lyon 1, UMR 5280 CNRS, 5 rue de la Doua, 69100 Villeurbanne, France</footnote>

<footnote n="2">² IC2MP UMR 7285 Université de Poitiers, Poitiers, France</footnote>

<footnote n="3">³ Institut CARMEN UMR6064 (ex lab. COBRA), Université de Rouen, INSA Rouen, CNRS, 1 Rue Tesnière, Mont-Saint-Aignan, 76821 Cedex, France</footnote>

<footnote n="4">⁴ Department of Chemistry and Pharmacology, QTC, Pontifical University of Chile, Santiago, Chile</footnote>

<footnote>ª Author to whom correspondence should be addressed: christophe.morell@univ-lyon1.fr</footnote>
to compute the LRF eigenvectors and eigenvalues. Basically, they reduced the six dimension of LRF to an infinite linear combination of functions (eigenvectors) of three coordinates.

Starting from their work, we propose in this paper a chemical interpretation of these eigenvectors as electron density deformation modes in an atom-condensed perspective. Indeed, the electron density deformation modes of LRF form a basis set (BS) on which any density deformation can be projected on. Hence, any density deformation can be written as linear combination of LRF eigenvectors. In addition, this opens a way to quantify the energy required to reshuffle the electron density, a process that occurs in any chemical reaction.

This paper is organized as follows: Sec. II deals with the theoretical background and introduces the diagonalization of the kernel. In addition, in this section, new equations relating the polarization with the variations of hardness and softness, as defined in C-DFT. Section III then describes in detail all the computations performed to illustrate this work. In Sec. IV, several examples are given to show how the electron density deformation modes can be used to rationalize the chemical behavior of organic functional groups. This paper ends with concluding remarks and introduces the next related study, namely, the investigation of the kernel of softness.

## II. THEORETICAL BACKGROUND

### A. Polarization energy and distortion modes from the linear response function

A molecular quantum system is fully characterized by the external potential $v$ generated by the nuclei, considered here as pure point charges, and the total number of electrons, $N$, so that the molecular energy, $E_{mol}$, can be written within the Born-Oppenheimer approximation, as

$$E_{mol} = E_e[N, v(\vec{r})] + V_{NN}, \quad (1)$$

where $E_e$ stands for the electronic energy and $V_{NN}$ for the repulsion between nuclei. At constant $N$, the $E_{mol}$ evolution is thus due to the changes in $v$. This occurs either during a rearrangement of the molecular backbone (isomerism) or at the early stage of a chemical reaction. The changes in the electric field felt by the electrons induce a polarization of the electron density, and the energy variation up to the second order is

$$dE_{mol}[\delta v(\vec{r})] = \int \left( \frac{\delta E_e}{\delta v(\vec{r})} \right)_N \delta v(\vec{r}) d\vec{r} + \frac{1}{2} \iint \left( \frac{\delta^2 E_e}{\delta v(\vec{r}) \delta v(\vec{r}')} \right)_N \delta v(\vec{r}) \delta v(\vec{r}') d\vec{r} d\vec{r}' + dV_{NN}. \quad (2)$$

Rayleigh-Schrödinger perturbation theory (RSPT) allows for identifying the first-order functional derivative with electron density $\rho(\vec{r})$, while the second-order derivative is, by definition, the static LRF, which is symmetric [i.e., $\chi(\vec{r}, r') \equiv \chi(r', \vec{r})$] in virtue of the Schwarz theorem for second derivatives. Equation (2) can thus be rewritten,

$$dE_{mol}[\delta v(\vec{r})] = \int \rho(\vec{r}) \delta v(\vec{r}) d\vec{r} + \frac{1}{2} \iint \chi(\vec{r}, \vec{r}') \delta v(\vec{r}) \delta v(\vec{r}') d\vec{r} d\vec{r}' + dV_{NN}. \quad (3)$$

The first right-hand side (RHS) term in Eq. (3) is the interaction of the undistorted density with the external potential, while the third RHS term is obviously the variation of the repulsion between nuclei. As for the second RHS term, it describes the interaction of the reshuffled density with the external potential since $\chi(\vec{r}, r') = \frac{\delta \rho(\vec{r})}{\delta v(\vec{r}')}$ evaluates the variation in electron density at point $\vec{r}$ when the external potential is changed at point $r'$ (or vice versa). This electron density change can be directly related to the variation $|\delta\psi_0|$ of the ground state wavefunction $|\psi_0\rangle$ (real-valued functions are here considered throughout this paper) in terms of the Hermitian electron density operator $\hat{\rho}(\vec{r}) = \sum_{i=1}^N \hat{\delta}(\vec{r} - \vec{r}_i)$ (where $\hat{\delta}$ denotes the Dirac multiplicative operator),

$$\begin{aligned} \delta\rho_0(\vec{r}) &= \langle\psi_0 + \delta\psi_0|\hat{\rho}|\psi_0 + \delta\psi_0\rangle - \rho_0(\vec{r}) \\ &\approx \langle\psi_0|\hat{\rho}|\psi_0\rangle + 2\langle\psi_0|\hat{\rho}|\delta\psi_0\rangle - \rho_0(\vec{r}) \\ &= 2\langle\psi_0|\hat{\rho}|\delta\psi_0\rangle. \end{aligned} \quad (4)$$

RSPT gives
$$|\delta\psi\rangle = \sum_{k \neq 0} \frac{\langle\psi_k|\delta\hat{v}|\psi_0\rangle}{E_0 - E_k} |\psi_k\rangle, \quad (4)$$

where the infinite sum runs over all electronic excited states, so that
$$\delta\rho_0(\vec{r}) = 2 \sum_{k \neq 0} \frac{\langle\psi_k|\delta\hat{v}|\psi_0\rangle}{E_0 - E_k} \langle\psi_0|\hat{\rho}(\vec{r})|\psi_k\rangle. \quad (5)$$

We now introduce the transition densities²⁷ as the $\rho_0^k(\vec{r}) = \langle\psi_0|\hat{\rho}(\vec{r})|\psi_k\rangle$ matrix elements. Noticing that $\langle\psi_0|\delta\hat{v}|\psi_k\rangle = \int \rho_0^k(\vec{r}') \delta v(\vec{r}') d\vec{r}'$, 
$$\delta\rho(\vec{r}) = 2 \sum_k c_k \rho_0^k(\vec{r}), \quad (6)$$

with
$$c_k = \frac{\int \rho_0^k(\vec{r}) \delta v(\vec{r}) d\vec{r}}{E_0 - E_k}, \quad (7)$$

which holds for any pure local small perturbation $\delta v(\vec{r})$. Note that these equations have been also at the origin of the recent investigation on the temperature and external fields extensions of conceptual DFT.²⁸

Finally, remarking that $\delta c_k / \delta v(\vec{r}') = \rho_0^k(\vec{r}') / (E_0 - E_k)$, taking the functional derivative with respect to $v(\vec{r}')$ of Eq. (6) straightforwardly gives²⁹,³⁰

$$\chi(\vec{r}, \vec{r}') = 2 \sum_k \frac{\rho_0^k(\vec{r}) \rho_0^k(\vec{r}')}{E_0 - E_k}. \quad (8)$$

The second-order energy variation, $E^{(2)} = \frac{1}{2} \iint (\chi(\vec{r}, \vec{r}') \delta v(\vec{r}) \delta v(\vec{r}') d\vec{r} d\vec{r}'$, which can also be called polarization energy, then simply reads

$$E^{(2)} = \sum_k \frac{\int \rho_0^k(\vec{r}) \delta v(\vec{r}) d\vec{r}'}{E_0 - E_k} = - \sum_k c_k^2 (E_k - E_0), \quad (9)$$
the usual expression from RSPT. It clearly appears that $E^{(2)}$ is negative for any type of potential perturbation, corresponding to the stabilization due to the relaxation of the electron density. It is also worth noticing that, while the introduction rooted our approach within the C-DFT framework, all equations reported so far are actually NOT stemming from DFT. Indeed, all derivations are based on wavefunction perturbation theory and not on DFT equations (nowhere the founding Hohenberg-Kohn theorems were used).

Starting from the latter equation, some of the present authors have proposed several indices to characterize the global and local polarizability of a molecule,³¹,³² and even at an atoms-in-molecules scale.³³,³⁴ We hope that the following sections may open new avenues to this topic.

Regarding numerical approaches and applications, it could be remarked that the described approach cannot be directly applied within the popular Kohn-Sham (KS) DFT framework since KS-DFT and its time-dependent extension only provide wavefunctions for the fictitious non-interacting electronic system. One has thus to rely alternatively on the coupled perturbed KS formalism that gives, after a tough derivation,

$$\chi(\vec{r}, \vec{r}') = \sum_{i,b} \sum_{j,c} a_{i,b,j,c} \phi_i(\vec{r}) \phi_b(\vec{r}') \phi_j(\vec{r}') \phi_c(\vec{r}), \quad (10)$$

where $\phi$ generically denotes a (real-valued) KS molecular orbital (MO), with indices $i$ or $j$ referring to an occupied KS MO, and $b$ or $c$ to a virtual one. The intricate $a_{i,b,j,c}$ elements include MO energies, as well as matrix elements for the exchange-correlation kernel (see Ref. 35 for the explicit expressions). Equation (10) is hence the first one that really belongs to C-DFT (more specifically to KS C-DFT).

In summary, LRF can be evaluated in both wavefunction and DFT contexts, so that it can be seen as a “universal” kernel that is not rooted in a given theoretical framework, which thus emphasizes its chemical content. However, its numerical evaluation depends on the framework: RSPP scheme for wavefunction methods, coupled-perturbed KS within KS-DFT since, as already mentioned; KS does not provide, even if the exact exchange-correlation functional was used, the (ground or excited) exact wavefunctions.

In the following, only KS DFT will be considered. Within this approach, a popular approximation, often coined “frozen orbital” in the chemistry community and “independent particle approximation” (IPA) in the physics one, simply reads for fully doubly occupied MOs (with Kronecker's $\delta$ symbols),

$$a_{i,b,j,c}^{\text{IPA}} = -\frac{4\delta_{ij}\delta_{bc}}{\varepsilon_c - \varepsilon_j}. \quad (11)$$

LRF then takes the simpler form

$$\chi(\vec{r}, \vec{r}') = -\sum_{i,b} \frac{4}{\varepsilon_b - \varepsilon_i} \phi_i(\vec{r}) \phi_b(\vec{r}') \phi_i(\vec{r}') \phi_b(\vec{r}), \quad (12)$$

which is the only one implemented in the software used for the numerical applications in Sec. IV.

As LRF is a symmetric definite negative kernel, it can be shown²⁶,³⁶ (as the extension of Mercer's theorem for non-compact sets applied to the definite positive $-\chi$ kernel) to be expandable, according to

$$\chi(\vec{r}, \vec{r}') = \sum_k \lambda_k \chi_k(\vec{r}) \chi_k(\vec{r}'), \quad (13)$$

where the eigenvalues $\lambda_k$ are negative and $\chi_k$ belong to a countable set of orthonormal functions,

$$\int \chi_k(\vec{r}) \chi_l(\vec{r}) d\vec{r} = \delta_{kl}. \quad (14)$$

As a side note, and in order to avoid possible misconception, it is instructive to come back to Eq. (8). Indeed, one may think at first sight that $\rho_0^k$ are actually $\chi_k$. However, while the electronic wavefunctions are mutually orthogonal, the associated transition densities are in general not orthogonal (of course, they can be orthogonalized using, for instance, the well-known Gram-Schmidt process. This will be the topic of a forthcoming paper). Consider as a single example the hydrogen atom, for which (considering the spherically symmetric wavefunctions) $\rho_0^1(\vec{r}) = \phi_{1s}(\vec{r})\phi_{2s}(\vec{r})$ and $\rho_0^2(\vec{r}) = \phi_{1s}(\vec{r})\phi_{3s}(\vec{r})$. The overlap between $\rho_0^k$ and $\rho_0^0$ can be analytically calculated and is found close to 0.0021. This example also shows, for the very same reason, that neither of the $\phi_i\phi_b$ values in Eq. (12) are the targeted eigenvalues. Nevertheless, these eigenvalues and eigenfunction offer a direct decomposition of the electron density change,

$$\delta\rho(\vec{r}) = \int \chi(\vec{r}, \vec{r}') \delta v(\vec{r}') d\vec{r}', \quad (15)$$

$$= \sum_k \left( \lambda_k \int \chi_k(\vec{r}') \delta v(\vec{r}') d\vec{r}' \right) \chi_k(\vec{r}). \quad (16)$$

The second-order correction to the energy is then given by

$$E^{(2)} = \frac{1}{2} \sum_k \lambda_k \iint \chi_k(\vec{r}) \chi_k(\vec{r}') \delta v(\vec{r}) \delta v(\vec{r}') d\vec{r} d\vec{r}', \quad (17)$$

$$= \frac{1}{2} \sum_k \lambda_k \left( \int \chi_k(\vec{r}) \delta v(\vec{r}) d\vec{r} \right)^2. \quad (18)$$

It is valuable to notice that this polarization energy can also be conveniently mapped, as proposed by some authors of this article, affording a visual representation of the associated chemical effects. Indeed, if we choose the potential variation proportional to a Dirac distribution centered in $\vec{R}$ according to $\delta v_R(\vec{r}) = \epsilon\delta(\vec{r} - \vec{R})$ (with $\epsilon \ll 1$), Eq. (18) gives

$$\delta E^{(2)}(\vec{R}) = \frac{\epsilon^2}{2} \sum_k \lambda_k \chi_k(\vec{R})^2. \quad (19)$$

This equation clearly allows for pointing in real space the positions where the electron density is most modified by the variation of the external potential. These positions are those of maximum magnitude of the quantity $\lambda_k\chi_k(\vec{R})^2$. This function combines an eigenvalue whose dimension is the reciprocal of an energy and an eigenvector that is a density. It is, therefore, alike a local softness. Hence, one can be hopeful that this quantity might be a good descriptor of local polarizability. Several authors have already made the link between local polarizability and local softness, lending credibility to this hypothesis.³⁷⁻³⁹

More generally, the LRF eigenfunctions actually form a complete BS. As a consequence, any local function $f$ of this kind can be written as a linear combination of the LRF eigenfunctions,
$$f(\vec{r}) = \sum_k d_k \chi_k(\vec{r}), \qquad (20)$$

with coefficients obtained by scalar products,

$$d_k = \int f(\vec{r}) \chi_k(\vec{r}) d\vec{r}. \qquad (21)$$

For instance, $f(\vec{r})$ could be the dual descriptor (DD), $\Delta f(\vec{r})$,⁴⁰ or any state specific DD ($\Delta f_i(\vec{r}) = \rho_i(\vec{r}) - \rho_0(\vec{r})$).⁴¹ Let us recall that the DD is the second derivative of $\rho$ when $N$ is varied at fixed $v$,

$$\Delta f(\vec{r}) := \frac{\delta^3 E_e}{\delta N^2 \delta v(\vec{r})} = \left( \frac{\partial^2 \rho(\vec{r})}{\partial N^2} \right)_{v(\vec{r})}. \qquad (22)$$

Conversely, LRF represents the $\rho$ change when $v$ is varied but at a fixed $N$ value. In a way, Eq. (20) allows making a direct connection between these two processes.

$f$ can also be the electron density deformation itself, $\delta\rho$. In this case,

$$d_k = \int \delta\rho(\vec{r}) \chi_k(\vec{r}) d\vec{r} = \lambda_k \int \delta v(\vec{r}) \chi_k(\vec{r}) d\vec{r}, \qquad (23)$$

where the last equality comes from Eq. (16) and the fact that any expansion in a BS is unique.

Even if it should have been expected, the latter equation is quite remarkable. Indeed, the variation of the external potential appears to be encoded within the variation of the electron density. Hohenberg and Kohn, in their landmark paper,¹² showed that a unique relation exists between the ground state electron density and the external potential. Equations (18) and (23) show that a unique relation exists between the variation of the external potential and the variation of the electron density. These equations also show that these eigenfunctions, $\chi_k(\vec{r})$, can then be considered as normal distortion modes of the electron density or electron density deformation/distortion Modes shortened as EDDM.

It should be noted that these EDDMs are fundamentally different in their construction and properties from the normal modes of electron displacements proposed by Nalewajski et al.⁴³ while studying fluctuations in atomic populations under “generic” chemical events [that is, enabling both $N$ and $v(\vec{r})$ variations], as here $N$ is by construction fixed—hence no electron gain/loss mode are present within EDDMs.

Using Eq. (23),

$$E^{(2)} = \frac{1}{2} \sum_k \lambda_k \times \left( \frac{d_k}{\lambda_k} \right)^2 = - \frac{1}{2} \sum_k \frac{d_k^2}{|\lambda_k|}. \qquad (24)$$

Equation (24) can be seen as the stabilization energy when a set of fractional densities $d_k^2$ are stabilized by the energy $\frac{1}{\lambda_k} < 0$. In conclusion, the knowledge of the set of coefficients $d_k$ provides a way to evaluate the magnitude of the stabilization energy for a given distortion of the electron density. Obviously, this equation is mathematically valid only for non-zero eigenvalues. However, integrating each side of Eq. (8) over the $\vec{r}^7$ variable,

$$\forall \vec{r} \int \chi(\vec{r}, \vec{r}') d\vec{r}' \propto \int \rho_0^k(\vec{r}') d\vec{r}' = \langle \psi_0 | \psi_k \rangle = 0, \qquad (25)$$

which can be recast as

$$\forall \vec{r} \int \chi(\vec{r}, \vec{r}') \times Cst(\vec{r}) d\vec{r}' = 0 \times Cst(\vec{r}). \qquad (26)$$

This would suggest that 0 is an eigenvalue corresponding to an eigenfunction that is a pure constant function (Cst). However, this last one cannot be normalized over whole space, does not convey any chemical information, and also will obviously make a vanishing contribution in Eq. (13). We thus exclude it from Eq. (24).

## B. More on dual descriptor and distortion modes from the linear response function

As the second derivative of $\rho$ with respect to $N$ (at fixed $v$), $\Delta f$ exactly integrates to 0 on whole space, as does $\delta\rho$ induced by any variation of the external potential $\delta v$ (at fixed $N$). It is thus tempting to interpret $\Delta f$ as a particular density distortion, even though we cannot provide a formal derivation of the following statement. It must, however, be noted that it has been shown previously that the thermal fluctuation of electron density at finite temperature is directly proportional to the 0 K dual descriptor,⁴⁴ hence supporting the use of the DD as a density distortion.

We, therefore, assume that there exists, for any molecular system, a $\delta v$ value that can generate a $\delta\rho$ identical to the $\Delta f$ function of this system. Equation (23) is nonetheless a clue that such potential exists. For this particular $\delta v$, we can use the results of Sec. II A,

$$E^{(2)} = -\frac{1}{2} \sum_k \frac{d_{k,\Delta f}^2}{|\lambda_k|}, \qquad (27)$$

with $\Delta f$ being decomposed as a linear combination of the electron density distortion modes with the following coefficients:

$$\Delta f(\vec{r}) = \sum_k d_{k,\Delta f} \times \chi_k(\vec{r}). \qquad (28)$$

Equation (27) is the first working equation allowing associating an energy variation to the DD. One can go ever further by relating the variation in the chemical hardness $\eta$ with the dual descriptor. Indeed, as within the canonical ensemble $\Delta f(\vec{r}) = (\frac{\delta\eta}{\delta v(\vec{r})})_N$, one has for any external potential variation,

$$\delta\eta = \int \Delta f(\vec{r}) \delta v(\vec{r}) d\vec{r}. \qquad (29)$$

If we now apply this general formula to the potential variation that generates the DD as the electron density variation,

$$\delta\eta = \sum_k d_{k,\Delta f} \int \chi_k(\vec{r}) \delta v(\vec{r}) d\vec{r}. \qquad (30)$$

Each integral is equal [see Eq. (23)] to $d_{k,\Delta f}/\lambda_k$, so that

$$\delta\eta = -\sum_k \frac{d_{k,\Delta f}^2}{|\lambda_k|} = 2E^{(2)} < 0. \qquad (31)$$

This means that modifying the electron density distortion following the DD will always decrease the hardness of a system. As a consequence, it will increase its reciprocal, namely, molecular softness $S$,

$$\delta S = - \frac{\delta \eta}{\eta^2} = S^2 \sum_k \frac{d_{k, \Delta f}^2}{|\lambda_k|}. \qquad (32)$$
As one can note, this result provides us with a complementary element to the maximum hardness principle.⁴⁵ Indeed, acknowledging that the DD reflects the electron density deformation associated to chemical reactivity, one could propose that “to afford chemical reactivity, molecules will naturally tend to lower their chemical hardness.”

The above-mentioned equations provide metrics to monitor the energetic variations due to a chemical reaction. Nonetheless, one must keep in mind that when working in the canonical ensemble, the number of electrons has to be constant. This limitation may prevent the use of these equations for comparing different molecules within a family. In addition, as far as we know, no principle of hardness or softness variation that would help predicting the direction of a chemical reaction has ever been reported so far.

### C. Fundamental electron density distortion mode

At this stage, one can wonder if all EDDMs have a significant contribution to the polarization energy or if only one or few of the EDDM contribute to most of it. From Eq. (18), it is plain to see that the total (negative) polarization energy increases with the (negative) eigenvalue of a specific mode and with the square of the coupling between a specific mode and the external potential creating the polarization. In addition, the EDDMs are orthogonal between themselves. As a consequence, with a well-chosen external potential, one can select a specific distortion mode and quench all other modes.

Pushing farther this line of argument, from Eq. (18), one can surmise that the electron density deformation modes (EDDM) are ranked by their eigenvalue. For an equivalent coupling between external potential and eigenvector, the lowest eigenvalue corresponds to the easiest way to distort the electron density. Thus, the set of eigenvectors and eigenvalues indicates how the electron density is distorted and how much easy it is to trigger this deformation. However, it is very unlikely that the external potential selects only one of the EDDMs. Actually, the results presented in Sec. III give evidence that it is very rare, at least for the DD decomposition. On the other hand, it is also very rare that a decomposition that spans over all the EDDMs. Generally, the DD is accurately decomposed along only few EDDMs.

### D. Electrical analogy

The relationship between the variation in electron density and the variation in external potential may lead us to propose an electrical analogy. Such analogies in C-DFT have proven to be useful for the heuristic design of new descriptors.⁴⁶ Here, we will consider a capacitor.⁴⁷⁻⁵¹ Indeed, for a simple capacitor made up of two conductors under electrical influence, the relationship between stored electric charge Q and voltage U is given by the $Q = C \times U$ relationship. For complex cases, the same equation can be written using a matrix formulation,

$$[Q] = [C][U]. \qquad (33)$$

On the other hand, the molecular (electronic) charge polarization $\delta\rho$ due to an applied potential $v(\vec{r})$ is given by

$$\delta\rho(\vec{r}) = \int \chi(\vec{r}, \vec{r}') \delta v(\vec{r}') d\vec{r}'. \qquad (34)$$

<figure><img src="image_6.png" /><figcaption>FIG. 1. Schematic representation of the analogy between a charged capacitor and a molecular dipole.</figcaption></figure>

<table><caption>TABLE I. Analogy between the chemical descriptors and the electronic properties of a capacitor.</caption><thead><tr><th>Chemical descriptors</th><th>Electronic variables</th></tr></thead><tbody><tr><td>$-\delta\rho(\vec{r})$</td><td>$\leftrightarrow$ $\delta Q$</td></tr><tr><td>$-\chi(\vec{r}, \vec{r}') d\vec{r} d\vec{r}'$</td><td>$\leftrightarrow$ $d^2 C$</td></tr><tr><td>$\delta v(\vec{r})$</td><td>$\leftrightarrow$ $\delta V$</td></tr></tbody></table>

The latter equation can be seen as the scalar product of the linear response function and the external potential, or as their matrices product. Consequently, this prompts us to propose an analogy between a charged capacitor and two atoms in a molecule whose electron density is polarized, see Fig. 1. Correspondence between variables are given in Table I.

To compare the properties of a capacitor to those of LRF, one can look at the analogous for the energy stored within a capacitor, with the work being

$$d^2 W_i = \frac{1}{2} d^2 C_i V^2. \qquad (35)$$

The stored energy $E_{cap}$ can be calculated, and using the analogy, the polarization energy can be computed,

$$E_{cap} = \iint \sum_{i} d^2 W_i = \frac{1}{2} \sum_{i} \iint d^2 C_i V^2, \qquad (36)$$

$$\leftrightarrow \frac{1}{2} \sum_{i} \lambda_{i} \left( \int \chi_{i}(\vec{r}) \delta v(\vec{r}) d\vec{r} \right)^{2} = E^{(2)}. \qquad (37)$$

This shows that the energy stored within a capacitor and the stabilization energy are analogous. It is important to note that while the energy stored in a capacitor is positive, the stabilization energy is negative. This is due to the opposite behaviors of electron density compared to capacitors. A negative potential will repulse the electron density, therefore, increasing the positive charges.

Furthermore, due to the fact that the sum of the eigenvalues gives the sum of the capacitance, the analogy leads to a physical interpretation of the diagonalization. The diagonalization of the atomically condensed (see Sec. III for more details) LRF corresponds to finding the analogous system of capacitors that are connected in parallel in the system,

$$-\chi(\vec{r}, \vec{r}') d\vec{r} d\vec{r}' = -\sum_{i} \lambda_{i} \chi_{i}(\vec{r}) \chi_{i}(\vec{r}') d\vec{r} d\vec{r}' \leftrightarrow \sum_{i} d^2 C_i = d^2 C, \qquad (38)$$

### III. COMPUTATIONAL DETAILS

All calculations have been performed using the AMS2023 package (and more specifically, the ADF⁵² module). The molecules have
been fully geometrically optimized in gas phase using the PBE0 hybrid exchange-correlation functional⁵³ along with the Slater type triple-ζ TZP basis set for all atoms. Scalar relativity effects have been added for all atoms.⁵⁴ The frozen core setting has been set to none, so that all core, valence, and virtual orbitals are taken into account in the LRF computation. The nature of the stationary points on the potential energy surface has been determined by computing the (harmonic) vibrational frequencies. Conceptual DFT descriptors have been calculated using the standard implementation in ADF. This means that the DD was evaluated using the FMO approach,

$$ \Delta f(\vec{r}) \approx \rho_{\text{LUMO}}(\vec{r}) - \rho_{\text{HOMO}}(\vec{r}). \quad (39) $$

This function has then been condensed⁵⁵ at the atomic level within the framework of Bader's Quantum Theory of Atoms-In-Molecules (QTAIM)⁵⁶ framework, according to

$$ \Delta f(A) = \int_{\Omega_A} \Delta f(\vec{r}) d\vec{r}, \quad (40) $$

where $\Omega_A$ denotes the QTAIM atomic basin for atom A. All atomic values for a given molecule of $N_{at}$ atoms have then been collected into one (molecular) $\vec{\Delta}f$ vector, which offers a kind of (finite size) coarse-grained representation of the DD function.

Similarly, LRF has been QTAIM condensed in the IPA approximation [Eq. (12), consistent with the use of the FMO approximation for the DD]. Then, a full exchange-correlation response is not taken into account. However, some exchange-correlation effects are actually introduced in an indirect way since the KS orbitals that are used in Eq. (12) are computed using the PBE0 potential. We refer the interested reader to the insightful paper by Wang et al. on this topic.⁵⁷ The QTAIM-condensed values are thus evaluated, according to

$$ \chi_{AB} = -\sum_{i,b} \frac{4}{\epsilon_b - \epsilon_i} \int_{\Omega_A} \phi_i(\vec{r}) \phi_b(\vec{r}) d\vec{r} \int_{\Omega_B} \phi_i(\vec{r}') \phi_b(\vec{r}') d\vec{r}'. \quad (41) $$

This evaluation only requires atomic molecular orbital overlaps and affords (diagonal) atomic (for $A = B$) and (off-diagonal) diatomic ($A \neq B$) elements. The condensed LRF matrix [$\chi_{AB}$] is thus of dimension $N_{at} \times N_{at}$. Note that as Bader's atoms constitutes an exhaustive non-overlapping partition of $\mathbb{R}^3$ and as a direct consequence of Eq. (25),

$$ \forall A \quad \chi(A, A) + \sum_{B \neq A} \chi(A, B) = 0. \quad (42) $$

Note that for the studied systems, the LRF computation takes only few seconds. The current implementation involves all molecular orbitals obtained at the SCF step. Obviously, it can be adapted in case of very large systems, for instance, by considering only orbitals lying in specific energy windows. So far, we have not found examples for which this approach was required.

Obviously, it is also technically possible to extend such a condensation scheme to other atomic partitions, but currently, only QTAIM condensation has been implemented in ADF software. However, we have some evidence (see the Appendix) that the condensed-LRF values are generally nicely correlated from one condensation scheme to another one, so that the trends obtained at the QTAIM level will remain valid with other atomic partitions.

Then, for every molecule, the condensed LRF matrix [$\chi_{AB}$] has been then diagonalized, giving $N_{at}$ eigenvectors, $\vec{\chi}_k$, obtained as the linear combination of the canonical vectors $\vec{e}_k$ (each one representing one atom) of the $\mathbb{R}^{N_{at}}$ Euclidean vector space. Obviously, this approach can benefit from the usual computational techniques used in quantum chemistry. For instance, for bigger molecules, the total or partial diagonalization of the condensed-LRF matrix could be efficiently achieved by the Lanczos or the Davidson algorithms. The $\vec{\Delta}f$ molecular DD vector can then be exactly decomposed on the $\vec{\chi}_k$ basis set according to

$$ \vec{\Delta}f = \sum_{k=1}^{N_{at}} d_k \vec{\chi}_k, \quad (43) $$

where each $d_k$ coefficients is obtained by the $\vec{\Delta}f \cdot \vec{\chi}_k$ standard scalar product in $\mathbb{R}^{N_{at}}$. A home-made script finally allows for graphically representing the LRF eigenvectors, in which the coefficient on each atom is symbolized by a sphere with a proportional radius. All values (unless otherwise explicitly stated) are given in atomic units (a.u.).

## IV. RESULTS AND DISCUSSION

The aim of this section is to highlight the potential interest of EDDMs to understand the reactivity of atoms in molecules. The examples chosen are organic chemistry textbook examples that exhibit a marked reactivity. Thus, methane, fluoro, chloro, and bromo methanes have been processed. The results are presented in Subsection IV B, after the special case of diatomic molecules that can be treated in part analytically (as an illustrative toy model for our methodological approach). The reactivity evolution of carbon within a series of carbonyl compounds is then presented in Subsection IV C.

### A. Diatomic molecules

As the consequence of the integration of the DD function to 0 over the whole space, one has for any AB diatomic molecule, $\vec{\Delta}f = \begin{pmatrix} -f \\ f \end{pmatrix}$. In addition, the LRF matrix for any AB molecule takes the following form [as a direct consequence of Eq. (42)] [$\chi$] = $\begin{pmatrix} -a & a \\ a & -a \end{pmatrix}$, where $a > 0$. This matrix is straightforwardly diagonalized by the $\vec{\chi}_1 = \begin{pmatrix} -1 \\ 1 \end{pmatrix}$ and $\vec{\chi}_2 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$ eigenvectors, associated to the $\lambda_1 = -2a < 0$ and $\lambda_2 = 0$ eigenvalues. The $\vec{\Delta}f$ decomposition is then trivial with $\vec{\Delta}f = f\vec{\chi}_1$. The corresponding energy and the hardness variation are given by the following formula: $E^{(2)} = -\frac{f^2}{4a}$ and $\delta\eta = -\frac{f^2}{2a}$.

As shown by this diatomic model and by our subsequent numerical results, the off-diagonal elements are all strictly positive. This means that all of these matrices belong to the Metzler matrices (MMs) family. From the mathematical point of view, it is known that any MM has an eigenvector in the nonnegative orthant. Here, this is confirmed by our calculations with the eigenvector $\vec{\chi}_2$ corresponding to the null eigenvalue. More generally, indeed, Eq. (42) is actually equivalent to the fact that the $\vec{u} = (11...1)^T$ vector verifies [$\chi]\vec{u} = 0$.

It can be further mentioned that while the constant function involved in Eq. (26) does not belong to the relevant Hilbert space [$L^2(\mathbb{R})$], $\vec{u}$ actually belongs to the relevant condensation Euclidean
<table><caption>**TABLE II.** Dual descriptor magnitude, LRF magnitude, and hardness evolution for hydrogen halides (in a.u.)</caption><thead><tr><th>Molecule</th><th>f</th><th>a</th><th>S</th><th>η<sub>0</sub></th><th>δη</th><th>η<sub>f</sub></th></tr></thead><tbody><tr><td>HF</td><td>0.39</td><td>0.38</td><td>2.30</td><td>0.435</td><td>-0.200</td><td>0.235</td></tr><tr><td>HCl</td><td>0.50</td><td>1.06</td><td>2.94</td><td>0.339</td><td>-0.117</td><td>0.222</td></tr><tr><td>HBr</td><td>0.47</td><td>1.41</td><td>3.31</td><td>0.302</td><td>-0.078</td><td>0.224</td></tr><tr><td>HI</td><td>0.43</td><td>1.75</td><td>6.10</td><td>0.164</td><td>-0.052</td><td>0.112</td></tr></tbody></table>

space. Furthermore, the diagonal elements are largely higher in absolute value than the off-diagonal ones. As a consequence of the Gershgorin circle theorem, the eigenvalues are quite close to the diagonal elements, and thus are negative, which is fully consistent with our numerical results.

Finally, it should be remarked that, mathematically speaking, the two following sequences are certainly not equivalent: (i) LRF diagonalization followed by atomic condensation and (ii) atomic condensation and LRF matrix diagonalization. This issue, however important, will not be pursued in this paper.

In Table II, the numerical values for *f*, *a*, *η*, and hardness variation are reported for hydrogen halides. As expected, hydrogen halides hardness decreases when one goes down the periodic classification of elements, namely from fluoride to iodide. Softness, defined as the reciprocal of hardness, increases along the chemical series. The *a* parameter follows the softness trend, emphasizing once more the relationship between polarizability and softness. It is also quite interesting to notice that one can compute the evolution of hardness when the density is deformed by a $\delta\rho$ equal to the DD. Interestingly, the biggest variation is obtained for hydrogen fluoride. This is quite expected: because of its initial large hardness value, to reach a similar electron density fluctuation in HF, one would expect is hardness to vary significantly.

## B. Methane and halo-methanes

We now turn our attention to more complex molecules. More specifically, in order to highlight the evolution of the methyl radical when it is substituted by halogen atoms, the EDDMs of fluoro, chloro, bromo, and iodo methanes have been computed, with methane taken as reference. The results are displayed in Tables III–VII and Fig. 2.

It does not come as a surprise that all molecules exhibit a constant eigenvector (same value for all the atoms) whose eigenvalue is zero. As already mentioned, this translates an unchanged electron density when the external potential is shifted by a constant

<table><caption>**TABLE III.** Linear response function eigenvectors and eigenvalues for methane.</caption><thead><tr><th>Atoms (k)</th><th>χ<sub>1</sub> (k)</th><th>χ<sub>2</sub> (k)</th><th>χ<sub>3</sub> (k)</th><th>χ<sub>4</sub> (k)</th><th>χ<sub>5</sub> (k)</th></tr></thead><tbody><tr><td>C1</td><td><b>0.894</b></td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.447</td></tr><tr><td>H2</td><td>-0.224</td><td>0.544</td><td>0.000</td><td>0.674</td><td>0.447</td></tr><tr><td>H3</td><td>-0.224</td><td>0.454</td><td>0.000</td><td>-0.738</td><td>0.447</td></tr><tr><td>H4</td><td>-0.224</td><td>-0.499</td><td>-0.707</td><td>0.032</td><td>0.447</td></tr><tr><td>H5</td><td>-0.224</td><td>-0.499</td><td>0.707</td><td>0.032</td><td>0.447</td></tr><tr><td>λ (a.u.)</td><td>-4.142</td><td>-1.161</td><td>-1.161</td><td>-1.161</td><td>0.000</td></tr></tbody></table>

<table><caption>**TABLE IV.** Linear response function eigenvectors and eigenvalues for fluoro methane.</caption><thead><tr><th>Atoms (k)</th><th>χ<sub>1</sub> (k)</th><th>χ<sub>2</sub> (k)</th><th>χ<sub>3</sub> (k)</th><th>χ<sub>4</sub> (k)</th><th>χ<sub>5</sub> (k)</th></tr></thead><tbody><tr><td>C1</td><td><b>0.889</b></td><td>0.000</td><td>0.000</td><td>0.096</td><td>0.447</td></tr><tr><td>F2</td><td><b>-0.129</b></td><td>0.000</td><td>0.000</td><td><b>-0.885</b></td><td>0.447</td></tr><tr><td>H3</td><td>-0.253</td><td>-0.813</td><td>0.073</td><td>0.263</td><td>0.447</td></tr><tr><td>H4</td><td>-0.253</td><td>0.469</td><td>0.668</td><td>0.263</td><td>0.447</td></tr><tr><td>H5</td><td>-0.253</td><td>0.344</td><td>-0.741</td><td>0.263</td><td>0.447</td></tr><tr><td>λ (a.u.)</td><td>-3.567</td><td>-1.115</td><td>-1.115</td><td>-1.013</td><td>0.000</td></tr></tbody></table>

<table><caption>**TABLE V.** Linear response function eigenvectors and eigenvalues for chloro methane.</caption><thead><tr><th>Atoms (k)</th><th>χ<sub>1</sub> (k)</th><th>χ<sub>2</sub> (k)</th><th>χ<sub>3</sub> (k)</th><th>χ<sub>4</sub> (k)</th><th>χ<sub>5</sub> (k)</th></tr></thead><tbody><tr><td>C1</td><td><b>0.886</b></td><td>0.124</td><td>0.000</td><td>0.000</td><td>0.447</td></tr><tr><td>Cl2</td><td><b>-0.341</b></td><td><b>0.827</b></td><td>0.000</td><td>0.000</td><td>0.447</td></tr><tr><td>H3</td><td>-0.181</td><td>-0.317</td><td>0.816</td><td>0.000</td><td>0.447</td></tr><tr><td>H4</td><td>-0.181</td><td>-0.317</td><td>-0.408</td><td>-0.707</td><td>0.447</td></tr><tr><td>H5</td><td>-0.181</td><td>-0.317</td><td>-0.408</td><td>0.707</td><td>0.447</td></tr><tr><td>λ (a.u.)</td><td>-4.395</td><td>-1.573</td><td>-1.107</td><td>-1.107</td><td>0.000</td></tr></tbody></table>

<table><caption>**TABLE VI.** Linear response function eigenvectors and eigenvalues for bromo methane.</caption><thead><tr><th>Atoms (k)</th><th>χ<sub>1</sub> (k)</th><th>χ<sub>2</sub> (k)</th><th>χ<sub>3</sub> (k)</th><th>χ<sub>4</sub> (k)</th><th>χ<sub>5</sub> (k)</th></tr></thead><tbody><tr><td>C1</td><td><b>0.870</b></td><td>0.209</td><td>0.000</td><td>0.000</td><td>0.447</td></tr><tr><td>Br2</td><td><b>-0.420</b></td><td><b>0.790</b></td><td>0.000</td><td>0.000</td><td>0.447</td></tr><tr><td>H3</td><td>-0.150</td><td>-0.333</td><td>0.442</td><td>-0.687</td><td>0.447</td></tr><tr><td>H4</td><td>-0.150</td><td>-0.333</td><td>0.374</td><td>0.726</td><td>0.447</td></tr><tr><td>H5</td><td>-0.150</td><td>-0.333</td><td>-0.816</td><td>-0.040</td><td>0.447</td></tr><tr><td>λ (a.u.)</td><td>-4.899</td><td>-1.721</td><td>-1.102</td><td>-1.102</td><td>0.000</td></tr></tbody></table>

<table><caption>**TABLE VII.** Linear response function eigenvectors and eigenvalues for iodo methane.</caption><thead><tr><th>Atoms (k)</th><th>χ<sub>1</sub> (k)</th><th>χ<sub>2</sub> (k)</th><th>χ<sub>3</sub> (k)</th><th>χ<sub>4</sub> (k)</th><th>χ<sub>5</sub> (k)</th></tr></thead><tbody><tr><td>C1</td><td><b>0.853</b></td><td>0.268</td><td>0.000</td><td>0.000</td><td>0.447</td></tr><tr><td>I2</td><td><b>-0.473</b></td><td><b>0.759</b></td><td>0.000</td><td>0.000</td><td>0.447</td></tr><tr><td>H3</td><td>-0.127</td><td>-0.342</td><td>0.816</td><td>0.000</td><td>0.447</td></tr><tr><td>H4</td><td>-0.127</td><td>-0.342</td><td>-0.408</td><td>-0.707</td><td>0.447</td></tr><tr><td>H5</td><td>-0.127</td><td>-0.342</td><td>-0.408</td><td>0.707</td><td>0.447</td></tr><tr><td>λ (a.u.)</td><td>-5.449</td><td>-1.813</td><td>-1.106</td><td>-1.106</td><td>0.000</td></tr></tbody></table>

("translation of the electron density," by analogy with the normal mode analysis in theoretical vibrational spectroscopy).

It is also worth noticing that the eigenvalues display a pattern that looks uncannily familiar. Indeed, methane has one low-lying eigenvalue, followed by three degenerated eigenvalues, respectively, -4.14 and -1.61. This pattern is similar to the MO eigenvalues of methane. Of course, this pattern is a consequence of the $T_d$

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623
<figure><img src="image_7.png" /><figcaption>FIG. 2. Schematic representation of the EDDMs ordered by row for CH₄, CH₃F, CH₃Cl, and CH₃I. The atomic components of the vectors are represented by a sphere of proportional radius. The color code is as follows: positive values are in white and negative in red. λₖ represents the eigenvalue of the corresponding eigenvector. For CH₃Br, the fourth mode is located on the hydrogen's in front.</figcaption></figure>

symmetry point group, with the low-lying eigenvector being of A₁ symmetry, while the three others are of T₂ symmetry.

The same applies to all halogenated methanes. All molecules exhibit four different eigenvectors, with two of them being degenerate. It is obviously the same pattern as that of the valence MO eigenvalues of these molecules that belong to the C₃v symmetry group point. The degenerate eigenvectors belong to the E irreducible representation, while the other eigenvectors are of the A₁ symmetry (see Fig. 2). Due to this degeneracy, the EDDMs are invariant with respect to a rotation of π/3 or 2π/3. Another expected feature is the decrease in the eigenvalues when going down the periodic table column of the halogen family. Taking the lowest eigenvalue, for example, the magnitude drops from ≈ −3.56 for fluorine to ≈ −4.39 for chlorine, to ≈ −4.89 for bromine and to ≈ −5.44 for iodine. The eigenvalues magnitude can accordingly be used as a pointer of the polarizability of the atoms involved in the eigenvector.

Now, turning to the eigenvectors, it is assumed that these functions can be seen as EDDMs. The contribution of each mode to LRF is proportional to its eigenvalue. As such, the most contributing mode is the one with the lowest eigenvalue. For all halogeno

methane and for methane itself, the lowest mode points to carbon as the most polarizing atom in the molecule. When said differently, carbon seems to be the softest atom, the one on which it is the easiest to modify density.

A closer look at the EDDM values show that the neighboring halogen atom, or the lack thereof in the methane case, appears to influence the carbon polarizability. In Table VIII, the magnitude of |λ₁χ₁(C)²| seems to be a good marker of the evolution of carbon polarizability. The carbon in CH₃Cl is less polarizable than that of methane, the carbon in CH₄ itself being less polarizable than that in CH₃Cl, being less polarizable than that in CH₃Br, and so forth for CH₃I.

<table><caption>TABLE VIII. Lowest electron density distortion mode condensed on carbon for methane, fluoro, chloro, and bromo methanes.</caption><thead><tr><th>Molecule</th><th>CH<sub>4</sub></th><th>CH<sub>3</sub>F</th><th>CH<sub>3</sub>Cl</th><th>CH<sub>3</sub>Br</th><th>CH<sub>3</sub>I</th></tr></thead><tbody><tr><td>|λ<sub>1</sub> × χ<sub>1</sub>(C)<sup>2</sup>|</td><td>3.31</td><td>2.82</td><td>3.44</td><td>3.70</td><td>3.96</td></tr></tbody></table>
<table><caption>**TABLE IX.** Lowest electron density distortion mode condensed on the halogen X = H, F, Cl, Br, I for methane, fluoro, chloro, and bromo methanes.</caption><thead><tr><th>Molecule</th><th>CH<sub>4</sub></th><th>CH<sub>3</sub>F</th><th>CH<sub>3</sub>Cl</th><th>CH<sub>3</sub>Br</th><th>CH<sub>3</sub>I</th></tr></thead><tbody><tr><td>|λ<sub>1</sub> × χ<sub>1</sub>(X)<sup>2</sup>|</td><td>0.20</td><td>0.06</td><td>0.51</td><td>0.86</td><td>1.22</td></tr></tbody></table>

The second lowest mode points the halogen as the most polarizable atom for CH<sub>3</sub>Cl and CH<sub>3</sub>Br, but not for CH<sub>3</sub>F. For the latter molecule, the mode with fluorine as the most polarizable atom has the highest eigenvalue. Obviously, this translates the resistance of fluorine to change its electron density, in other words its hardness. Another way to look at it is to compare the values of the quantity |λ<sub>1</sub>χ<sub>1</sub>(X)<sup>2</sup>| condensed on the halogen atom from the first EDDM. The values are reported in Table IX. Its first entry corresponds to the value on one hydrogen atom in lieu of halogen. It is plain that hydrogen is more polarizable than fluorine, but less than the other halogens. As expected, the polarizability of the halogen increases if one goes down the periodic table.

The last part of our analysis focuses on the DD decomposition over the EDDMs. The values are reported in Table X. It is first interesting to notice that Δf projects mainly on the fourth EDDM for CH<sub>3</sub>F, while it projects on two EDDMs, the first and the second, for both CH<sub>3</sub>Cl and CH<sub>3</sub>Br. For CH<sub>3</sub>F, this means that the external potential at the origin of this particular deformation of the electron density is collinear to χ<sub>4</sub>(r̄). For CH<sub>3</sub>Cl, CH<sub>3</sub>Br, and CH<sub>3</sub>I,

<table><caption>**TABLE X.** Dual descriptor decomposition over the EDDMs.</caption><thead><tr><th>Δf</th><th>χ<sub>1</sub></th><th>χ<sub>2</sub></th><th>χ<sub>3</sub></th><th>χ<sub>4</sub></th><th>χ<sub>5</sub></th><th>E<sup>(2)</sup> (eV)</th></tr></thead><tbody><tr><td>f<sup>(2)</sup><sub>(CH<sub>3</sub>F)</sub></td><td>0.100</td><td>0.000</td><td>0.000</td><td><b>-0.809</b></td><td>0.000</td><td>-0.324</td></tr><tr><td>f<sup>(2)</sup><sub>(CH<sub>3</sub>Cl)</sub></td><td><b>-0.555</b></td><td><b>-1.051</b></td><td>0.000</td><td>0.000</td><td>0.000</td><td>-0.386</td></tr><tr><td>f<sup>(2)</sup><sub>(CH<sub>3</sub>Br)</sub></td><td><b>-0.804</b></td><td><b>-0.721</b></td><td>0.000</td><td>0.000</td><td>0.000</td><td>-0.217</td></tr><tr><td>f<sup>(2)</sup><sub>(CH<sub>3</sub>I)</sub></td><td><b>-0.917</b></td><td><b>-0.584</b></td><td>0.000</td><td>0.000</td><td>0.000</td><td>-0.171</td></tr></tbody></table>

the external potential is a combination of χ<sub>1</sub>(r̄) and χ<sub>2</sub>(r̄). One can also wonder whether the EDDM number indicates the kind of electronic control the molecule complies with. Indeed, the fact that Δf for CH<sub>3</sub>F projects mainly on the fourth EDDM, corresponding to the highest eigenvalue, could indicate that reactions on CH<sub>3</sub>F are neither controlled by charge transfer (MO overlap) nor by density polarization.

### C. Carbonyl compounds

In Fig. 3, the EDDMs for ethanamide, ethanoic acid, methyl methanoate, ethanoyl chloride, and ethanal are reported along with their eigenvalues.

Just as for halogeno methanes, a constant eigenvector with zero eigenvalue exists, but it has been removed from the picture as it does not provide any relevant chemical information. It

<figure><img src="image_8.png" /><figcaption>**FIG. 3.** Schematic representation of the EDDMs ordered by row for ethanamide, ethanoic acid, methyl methanoate, ethanoyl chloride, and ethanal. The atomic component of a vector is represented by a sphere of proportional radius. Color code: positive values are in white and the negative values are in red. λ<sub>k</sub> represents the eigenvalue of the corresponding eigenvector.</figcaption></figure>
<table><caption>TABLE XI. First electron density distortion mode condensed on carbon are used for ethanal and acyl chloride. Second electron density distortion mode condensed on carbon for ethanoic acid, formyl acetate, and ethanamide.</caption><thead><tr><th>Molecule</th><th>CH<sub>3</sub>CHO</th><th>CH<sub>3</sub>COCl</th><th>CH<sub>3</sub>COOH</th><th>HCOOCH<sub>3</sub></th><th>CH<sub>3</sub>CONH<sub>2</sub></th></tr></thead><tbody><tr><td>|λ<sub>1/2</sub> × χ<sub>1/2</sub>(C)<sup>2</sup>|</td><td>2.3</td><td>2.2</td><td>1.9</td><td>1.2</td><td>1.6</td></tr><tr><td>|λ<sub>1/2</sub> × χ<sub>1/2</sub>(O)<sup>2</sup>|</td><td>0.9</td><td>0.8</td><td>1.6</td><td>1.5</td><td>0.9</td></tr></tbody></table>

is also worth noticing that EDDMs involving the atoms of the methyl group in the carbonyl family have the same shape as the EDDMs in halogeno methanes. For instance, $χ_1, χ_5, χ_6$ of methyl methanoate (HCOOCH$_3$) or ethanamide (CH$_3$CONH$_2$) exhibit the same pattern as $χ_1, χ_2, χ_3$ of any halogeno methane (CH$_3$X). Interestingly, the magnitude of the eigenvalues appears to be partly due to the shape of the eigenvector and to the atoms involved. This is pretty obvious when one compares $λ_2(CH_3X) ≈ λ_3(CH_3X) ≈ -1.1$ to $λ_5(HCOOCH_3) ≈ -1.1$ and $λ_6(HCOOCH_3) ≈ -1.08$.

We can now inspect the reactivity of these functional groups. The carbon atoms of the carbonyl group are known to exhibit a marked affinity for nucleophiles. It is well-established that this affinity can be mitigated by the substituents bonded to the carbonyl groups. Organic chemists have established a rule of thumb based on the stabilization or the enhancement of the electrophilicity of the carbon atom, i.e., the more electron donating the substituent, the less reactive the carbon. Consequently, the order of reactivity toward a nucleophilic attack goes as follows: acyl chloride > aldehyde > carboxylic acid > ester > amide, with the last three compounds as almost unreactive. It is the reason why esterification is often used to protect a C=O bond in organic synthesis.

The EDDM associated with the carbonyl group is the first EDDM for ethanal and ethanoyl chloride, while for methyl methanoate, ethanoic acid, and ethanamide, it is the second. It is worth noticing the switching between the first and second EDDM from CH$_3$COOH going upward in the table. Thus, the carbonyl group is the most easily polarizable group for the first two molecules. Looking only on this EDDM gives the stabilization of a nucleophilic approach on the carbon and, therefore, an insight into reactivity.

The eigenvalues corresponding to the EDDMs of the C=O bond follow an almost similar order to the rule of thumb, with stabilization in the following order: ethanal -4.4 > ethanoyl chloride -4.3 > methyl methanoate -3.5 > ethanoic acid -3.3 > ethanamide -3.3. It is also noticeable that the associated eigenvalues can be separated into two different groups: ethanal and ethanoyl chloride have eigenvalues around -4.4, while the other have eigenvalues around -3.3.

Looking now at an atomic polarizability, the contribution of the carbon atom of the carbonyl group is evaluated. The magnitude of $|λ_i| × χ_i(X)^2$ is calculated for the first EDDM for CH$_3$CHO and for CH$_3$COCl, and for the second EDDM for the remaining molecules of the family. The figures are reported in Table XI. The chosen modes are the ones involving most the C=O double bond. It is plain to see that the EDDM that involves the carbonyl atoms (carbon and oxygen) is almost centered entirely on those atoms. No other atom notably contributes to this eigenvector.

Moreover, the EDDM exhibits an opposite sign on carbon and oxygen. This means that the electron density evolves in opposite

direction, one atom gaining density when the other undergoes a loss. This could be expected from a nucleophilic attack targeting the carbonyl carbon. As expected, this one is the most polarizable, at least for the first two molecules of the series. The methyl carbon seems to be the more polarizable atom for the latter three molecules. However, one can wonder whether this surprising mode is not biased by condensation. In addition, for methyl methanoate, one notices that oxygen is the most polarizable atom. Interestingly, the magnitude of $|λ_i| × χ_i(X)^2$ seems to provide a good scale of the carbonyl carbon reactivity along the chemical series.

Nonetheless, one must be cautious as the descriptor is defined within the canonical ensemble and, therefore, comparisons between molecules with different number of electrons are questionable. This is why in the following sections, only isomers are compared.

## D. Comparing the reactivity of isomers with EDDMs

### 1. C$_2$H$_4$O isomers

Three isomers of C$_2$H$_4$O exist, namely, ethanal, ethenol, and ethylene oxide (oxacyclopropane). They are ranked by decreasing stability or increasing order of reactivity in Fig. 4.

Even if ethanal EDDMs have been already discussed in Sec. III, they are nonetheless represented again along with EDDMs of ethenol and ethylene oxide for the sake of comparison in Fig. 5.

The first interesting feature to notice is the ranking of the lowest eigenvalue, and $λ_1(CH_3CHO) ≈ -6.5$ for ethenol is much lower than that for ethylene oxide, $λ_1(CH_2-O-CH_2) ≈ -4.7$, which is in turn lower than that for ethanal, $λ_1(CH_3CHO) ≈ -4.5$. This ordering seems to reflect the order of reactivity between isomers. The corresponding eigenvector also correctly displays the easiest way to distort the electron density. Thus, for this set of molecules, the lowest eigenvector points the C=O double bond, the carbons of the cycle, and the C=C double bond, respectively, for ethanal, ethylene oxide, and ethenol.

As carbonyl compounds have already been discussed, we are focusing on ethylene oxide. The order of reactivity points the carbons of the ring as the most reactive functional group within the molecule, which is what is experimentally observed. Indeed, the main reactivity consists in the opening of the ethylene oxide ring by nucleophilic attack targeting these carbons. The second and third

<figure><img src="image_9.png" /><figcaption>FIG. 4. Ethanal, ethylene oxide, and ethenol by decreasing order of stability.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 5. Schematic representation of the EDDMs of ethanal (top), ethylene oxide (middle), and ethenol (bottom). The atomic component for a vector is represented by a sphere of proportional radius. Color code: positive values are in white and negative values are in red. λₖ represents the eigenvalue of the corresponding eigenvector.</figcaption></figure>

vectors point the oxygen as the second most reactive atom. It is very likely that these vectors reflect on the complexation of oxygen of the ethylene oxide ring by Lewis acids.

As for the reactivity of ethenol, it is well known that the molecule has a very short lifetime as it transforms into ethanal by the called keto-enol prototropy process, as shown in Fig. 6.

This isomerism is almost always shifted toward the ketone isomer. Only few examples of enol are known. Phenol is one, in which the aromaticity favored the enol form. In the chosen example, the prototropy between ethanal and ethenol, the stable isomer is ethanal. Even though the molecule is not stable, one can expect C=C to be the most reactive functional group. This is what is represented by the first eigenvector. This eigenvector is anti-symmetric with respect to the C=C double bond. The EDDM exhibits opposite sign for the carbons. Therefore, when one carbon gains electron density, the other loses electron density. The magnitude of the eigenvector on the

terminal carbon χ(C₁) ≈ -0.76 is bigger than that on the substituted carbon, a prototypical application of the famous Markovnikov rule.

Looking now at the ethenol-ethanal prototropy, the mechanism likely requires the breaking of the oxygen hydrogen bond, freeing a proton that shifts toward the terminal carbon. The equilibrium is totally shifted toward the formation of ethanal. The mechanism as it is described in organic textbooks is schematized in Fig. 7. Regarding the electron reshuffling during the prototropy, one can surmise that the whole mechanism can be viewed as a combination of the first and second EDDMs. Indeed, from the electron view point, there is a shift of the electron pair of the oxygen hydrogen bond toward the substituted carbon, inducing a loss in density at the oxygen end and a gain at the substituted carbon end. This also induces a gain in electron at the terminal carbon by polarization.

This first part ends up in the formation of the C=O bond. The process also frees a proton. This part of the reshuffling would be well represented by the second EDDM. In synchronicity, the double bond shift toward the terminal carbon provides the electron density necessary to capture a proton. This second process would be best represented by the first EDDM. All in all, the evolution of the density could be qualitatively approximated as

$$|\Delta\rho(k)| \approx |\chi_1(k) - \chi_2(k)|. \qquad 851$$

In conclusion, the electron density variation in the prototropy between ethanal and ethenol can be, at least qualitatively, given by a combination of EDDMs. A quantitative relation would necessitate the computation of the changes in external potential, due to the

<figure><img src="image_11.png" /><figcaption>FIG. 6. Prototropy reaction interconverting ethenol and ethanal.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 7. Prototropy reaction interconverting ethenol and ethanal.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 8. Cyanamide (left) and carbodiimide (right) isomer.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 9. EDDMs of cyanamine (top) and carbodiimide (bottom).</figcaption></figure>

motion of the proton, and the application of Eq. (23) for each and every EDDM.

## 2. Carbodiimide cyanamide

The isomery between cyanamide and carbodiimide is represented in Fig. 8. It is well-established that cyanamide is the stablest isomer, but not by much since carbodiimide has been detected in deep space by its maser emission.⁵⁸

Reactivity of cyanamide can be guessed from its chemical formula. Basically, it is a CN nitrile functional group bonded to an amine. However, since the lone pair of the amino nitrogen is conjugated with the nitrile moiety, its nucleophilicity or basicity is likely weaker than a regular amine, just like an amide. On the other hand, as the nitrile group is connected to a nitrogen atom, more electronegative than a carbon atom, its reactivity is likely stronger than a regular nitrile. Using the same line of arguments, the reactivity of carbodiimide can also be predicted. It can be considered a heteroalene, i.e., a system with two vicinal double bonds in which the two terminal carbons have been replaced by nitrogen. The two π systems are orthogonal and, therefore, without noticeable overlap. Overall, carbodiimide should behave like an allene with a lone pair on each terminal nitrogen atom. The EDDMs of both systems are presented in Fig. 9.

Quite interestingly, the distribution of the first EDDM in both cases can be rather easily matched to mesomeric forms, as shown in Fig. 10.

Furthermore, they are in rather good agreement with the expected reactivity. It is interesting to notice first the gap between the first two eigenvalues of cyanamide, $\lambda_1(\text{NH}_2\text{CN})$

$\approx -6.9\lambda_2(\text{NH}_2\text{CN}) \approx -3.0$, as if the first reactivity (nitrile) is more pronounced than the second (amino group). Another striking feature is the asymmetry of the first EDDM in which, the coefficient on nitrogen (0.73) is higher than that on carbon ($-0.69$). This asymmetry likely translates the polarization of the $C\equiv N$ triple bond due to the lone pair carried by the imino nitrogen.

Reactivity of carbodiimide is actually well-documented. The heteroallene carbon is generally targeted by nucleophiles, and the two terminal nitrogens by electrophiles. Overall, the molecules behave as much as a nucleophile as an electrophile (comparable to the small difference in the first two eigenvalues seems to reflect this experimental fact. It is also plain to see that the first EDDM points carbon as the most reactive atom, while the second EDDM points nitrogen as the most reactive atom. Furthermore, the second EDDM is only located on N, in line with the lone-pair character of the associated electronic system, while the second is more delocalized, in line with the expected delocalizability of the associated $\pi$ electron cloud. Satisfyingly, the first two EDDMs report quite clearly the reactivity of the carbodiimide functional group.

## V. CONCLUSIONS

In the past decade, a lot of attention has been paid to the linear response function and especially to its eigenfunctions and eigenvalues. In this paper, a chemical meaning is given to these functions as Electron Density Distortion Modes (EDDMs), i.e., the different ways in which the electron density of a molecule can be deformed. Every distortion of the electron density of a molecule can be written as a linear combination of EDDMs. Thus, EDDMs provide a new mathematical framework to describe chemical reactivity and selectivity. More specifically, the EDDMs ascribed to the lowest eigenvalues correspond to the easiest way to distort the electron density, and therefore, the total spectrum of eigenvalues sorts by decreasing order the reactivity of a molecule.

Several examples were discussed to show how these eigenfunctions can open new avenues to predict the chemical behavior of common organic functional groups. Limitations are also highlighted. Using the condensed-to-atoms linear response function matrix instead of the actual six-dimensional function prevents form accessing stereo selectivity information. A step forward would be to easily compute the non-condensed eigenvectors, and will be the subject of a future investigation.

In addition, as the linear response function is defined within the canonical ensemble, only systems with the same number of

<caption>FIG. 10. Proposed mesomeric form accounting for the two first EDDM of cyanamide (top and center) and for the first EDDM of carbodiimide (bottom).</caption>
electrons could be safely compared and enter the scope of this descriptor. One can, therefore, suggest that the kernel of softness,⁵⁹ which is the linear response function counterpart within the C-DFT grand canonical ensemble, may provide a more complete view of chemical reactivity and selectivity of a molecule, by adding charge transfer to polarization. It will also relate all these concepts to Kohn's nearsightedness of matter principle.⁶⁰ By this way, all possible reactivity behaviors might be properly studied.

## ACKNOWLEDGMENTS

This work was funded by the ANR "QCTES" grant (Project No. 23-CE29-0023). This work was partially supported by the University of Rouen Normandy, INSA Rouen Normandy, the "Center National de la Recherche Scientifique" (CNRS), the European Regional Development Fund (ERDF), Labex SynOrg (Grant No. ANR-11-LABX-0029), Carnot Institut I2C, the graduate school for research XL-Chem (Grant No. ANR-18-EURE-0020 XL CHEM), and the "Région Normandie." F.G. and C.M. acknowledge the National Universities Commission (C.N.U) for one semester sabbatical in 2024. C.M. also thanks A.T.L. for the invitation to a research stay in Chile.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Rémi Grincourt:** Conceptualization (equal); Formal analysis (equal); Software (lead); Visualization (equal). **Guillaume Hoffmann:** Data curation (equal); Writing - review & editing (equal). **Frédéric Guégan:** Data curation (equal); Formal analysis (equal); Writing - review & editing (equal). **Vincent Tognetti:** Conceptualization (equal); Formal analysis (equal); Writing - original draft (equal); Writing - review & editing (equal). **Laurent Joubert:** Data curation (equal); Writing - original draft (equal); Writing - review & editing (equal). **Henry Chermette:** Conceptualization (equal); Data curation (equal); Writing - original draft (equal). **Alejandro Toro-Labbé:** Data curation (equal); Formal analysis (equal); Writing - original draft (equal). **Christophe Morell:** Conceptualization (lead); Formal analysis (lead); Writing - original draft (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## APPENDIX: INFLUENCE OF THE COMPUTATIONAL PROTOCOL

In this appendix, we compare our numerical values to those reported by Sablon et al.⁶¹ In this paper, the full matrices for six small molecules (CO, H₂O, HCN, NNO, NH₃, H₂CO) are provided using the Becke fuzzy atoms partition for PBE/6-31+G(d) obtained from single point (SP) calculations on the B3LYP/6-311++G(d,p) optimized geometries. We thus optimized the geometries at the B3LYP/TZP level, followed by SP at the PBE/DZP one and QTAIM condensation (it should be noted that, due to the use of Slater basis sets in ADF rather than Gaussian ones, it is not possible to exactly reproduce the same computational protocol) without

<figure><img src="image_15.png" /><figcaption>FIG. 11. Correlation between the eigenvalues obtained at our level of theory with those reported by Sablon et al.⁶¹</figcaption></figure>

relativistic corrections. Figure 11 shows a high correlation ($R^2$ = 0.96) for the obtained non-null eigenvalues between both approaches.

However, about the role of the computational parameters, among the expected dominant contributions in Eq. (12), there is the one involving HOMO and LUMO as the occupied and unoccupied pair. It is then controlled by the HOMO-LUMO gap energy at the denominator. It is well-known that this gap may depend on the use exchange-correlation functional (in general, it tends to increase with the mixing of exact exchange). However, this is not the purpose of this paper to investigate the influence of the functional choice on the numerical results. We nevertheless believe that the conclusions obtained within the same family of compounds are less dependent of this choice.

## REFERENCES

¹P. Geerlings, F. de Proft, and W. Langenaeker, "Conceptual density functional theory," *Chem. Rev.* **103**, 1793-1874 (2003).

²D. Chakraborty and P. K. Chattaraj, "Conceptual density functional theory based electronic structure principles," *Chem. Sci.* **12**, 6264-6279 (2021).

³P. Geerlings, E. Chamorro, P. K. Chattaraj, F. de Proft, J. L. Gázquez, S. Liu, C. Morell, A. Toro-Labbé, A. Vela, and P. Ayers, "Conceptual density functional theory: Status, prospects, issues," *Theor. Chem. Acc.* **139**, 36 (2020).

⁴R. G. Parr, R. A. Donnelly, M. Levy, and W. E. Palke, "Electronegativity: The density functional viewpoint," *J. Chem. Phys.* **68**, 3801-3807 (1978).

⁵R. G. Pearson, "Chemical hardness and density functional theory," *J. Chem. Sci.* **117**, 369-377 (2005).

⁶W. Yang and R. G. Parr, "Hardness, softness, and the fukui function in the electronic theory of metals and catalysis," *Proc. Natl. Acad. Sci. U. S. A.* **82**, 6723-6726 (1985).

⁷R. G. Parr and W. Yang, "Density functional approach to the frontier-electron theory of chemical reactivity," *J. Am. Chem. Soc.* **106**, 4049-4050 (1984).

⁸R. T. Sanderson, "An interpretation of bond lengths and a classification of bonds," *Science* **114**, 670-672 (1951).

⁹R. G. Pearson, "The principle of maximum hardness," *Acc. Chem. Res.* **26**, 250-255 (1993).

¹⁰F. A. Carey and R. J. Sundberg, *Advanced Organic Chemistry Part A* (Springer, 2007).

¹¹G. Klopman, "Chemical reactivity and the concept of charge- and frontier-controlled reactions," *J. Am. Chem. Soc.* **90**, 223-234 (1968).

¹²L. Salem, "Intermolecular orbital theory of the interaction between conjugated systems. I. General theory," *J. Am. Chem. Soc.* **90**, 543-552 (1968).

¹³K. Fukui, "Role of frontier orbitals in chemical reactions," *Science* **218**, 747-754 (1982).

980
981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

999

1000

1001

1002

1003

1004

1005

1006

1007

1008

1009

1010

1011

1012

1013

1014

1015

1016

1017

1018

1019

1020

1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147

¹⁴K. Fukui and H. Fujimoto, "An MO-theoretical interpretation of the nature of the chemical reaction. I. Partitioning analysis of the interaction energy," *Bull. Chem. Soc. Jpn.* **41**, 1989 (1968).

¹⁵P. W. Ayers and R. G. Parr, "Variational principles for describing chemical reactions: The fukui function and chemical hardness revisited," *J. Am. Chem. Soc.* **122**, 2010–2018 (2000).

¹⁶R. G. Parr and W. Yang, *Density-Functional Theory of Atoms and Molecules, International Series of Monographs on Chemistry* (Oxford University Press, 1994).

¹⁷P. K. Chattaraj, "Chemical reactivity and selectivity: Local HSAB principle versus frontier orbital theory," *J. Phys. Chem. A* **105**, 511–513 (2001).

¹⁸B. Mussard and J. G. Ángyán, "Relationships between charge density response functions, exchange holes and localized orbitals," *Comput. Theor. Chem.* **1053**, 44–52 (2015), part of Special Issue: Understanding Structure and Reactivity from Topology and Beyond.

¹⁹P. Geerlings, S. Fias, Z. Boisdenghien, and F. de Proft, *Chem. Soc. Rev.* **43**, 4989–5008 (2014).

²⁰N. Sablon, F. de Proft, and P. Geerlings, "The linear response kernel of conceptual DFT as a measure of electron delocalisation," *Chem. Phys. Lett.* **498**, 192–197 (2010).

²¹Z. Boisdenghien, C. Van Alsenoy, F. de Proft, and P. Geerlings, "Evaluating and interpreting the chemical relevance of the linear response kernel for atoms," *J. Chem. Theory Comput.* **9**, 1007–1015 (2013).

²²Z. Boisdenghien, S. Fias, F. Da Pieve, F. de Proft, and P. Geerlings, "The polarisability of atoms and molecules: A comparison between a conceptual density functional theory approach and time-dependent density functional theory," *Mol. Phys.* **113**, 1890–1898 (2015).

²³A. Scherrer, V. Verschinin, and D. Sebastiani, "Eigensystem representation of the electronic susceptibility tensor for intermolecular interactions within density functional theory," *J. Chem. Theory Comput.* **8**, 106–111 (2012).

²⁴C. Dresler and D. Sebastiani, "Reduced eigensystem representation of the linear density–density response function," *Int. J. Quantum Chem.* **120**, e26085 (2020).

²⁵A. Scherrer and D. Sebastiani, "Moment expansion of the linear density–density response function," *J. Comput. Chem.* **37**, 665–674 (2016).

²⁶P. Senet, "Variational principle for eigenmodes of reactivity in conceptual density functional theory," *ACS Omega* **5**, 25349–25357 (2020).

²⁷H. C. Longuet-Higgins, "The electronic states of composite systems," *Proc. R. Soc. London, Ser. A* **235**, 537–543 (1956).

²⁸M. Franco-Pérez, F. Heidar-Zadeh, P. W. Ayers, F. De Proft, A. Vela, J. L. Gázquez, and P. Geerlings, "Temperature and external fields in conceptual density functional theory," *Chem. Sci.* **15**, 20090–20121 (2024).

²⁹M. Berkowitz and R. G. Parr, "Molecular hardness and softness, local hardness and softness, hardness and softness kernels, and relations among these quantities," *J. Chem. Phys.* **88**, 2554–2557 (1988).

³⁰S. Liu, T. Li, and P. W. Ayers, "Potentialphilicity and potentialphobicity: Reactivity indicators for external potential changes from density functional reactivity theory," *J. Chem. Phys.* **131**, 114106 (2009).

³¹F. Guégan, T. Pigeon, F. de Proft, V. Tognetti, L. Joubert, H. Chermette, P. W. Ayers, D. Luneau, and C. Morell, "Understanding chemical selectivity through well selected excited states," *J. Phys. Chem. A* **124**, 633–641 (2020).

³²F. Guégan, V. Tognetti, J. I. Martínez-Araya, H. Chermette, L. Merzoud, A. Toro-Labbé, and C. Morell, "A statistical thermodynamics view of electron density polarisation: Application to chemical selectivity," *Phys. Chem. Chem. Phys.* **22**, 23553–23562 (2020).

³³V. Tognetti, F. Guégan, D. Luneau, H. Chermette, C. Morell, and L. Joubert, "Structural effects in octahedral carbonyl complexes: An atoms-in-molecules study," *Theor. Chem. Acc.* **136**, 85 (2017).

³⁴S. Menant, F. Guégan, V. Tognetti, L. Merzoud, L. Joubert, H. Chermette, and C. Morell, "Polarisation of electron density and electronic effects: Revisiting the carbon–halogen bonds," *Molecules* **26**, 6218 (2021).

³⁵P. Geerlings, S. Fias, Z. Boisdenghien, and F. De Proft, "Conceptual DFT: Chemistry from the linear response function," *Chem. Soc. Rev.* **43**, 4989–5008 (2014).

³⁶S. Fias, P. W. Ayers, F. De Proft, and P. Geerlings, "Properties of the density functional response kernels and its implications on chemistry," *J. Chem. Phys.* **157**, 114102 (2022).

³⁷A. Vela and J. L. Gazquez, "A relationship between the static dipole polarizability, the global softness, and the fukui function," *J. Am. Chem. Soc.* **112**, 1490–1492 (1990).

³⁸T. K. Ghanty and S. K. Ghosh, "Correlation between hardness, polarizability, and size of atoms, molecules, and clusters," *J. Phys. Chem.* **97**, 4951–4953 (1993).

³⁹P. Fuentealba, "A local model for the hardness kernel and related reactivity parameters in density functional theory," *J. Chem. Phys.* **103**, 6571–6575 (1995).

⁴⁰C. Morell, A. Grand, and A. Toro Labbé, "New dual descriptor for chemical reactivity," *J. Phys. Chem. A* **109**, 205–212 (2005).

⁴¹V. Tognetti, C. Morell, P. W. Ayers, L. Joubert, and H. Chermette, "A proposal for an extended dual descriptor: A possible solution when Frontier molecular orbital theory fails," *Phys. Chem. Chem. Phys.* **15**, 14465–14475 (2013).

⁴²P. Hohenberg and W. Kohn, "Inhomogeneous electron gas," *Phys. Rev.* **136**, B864–B871 (1964).

⁴³R. F. Nalewajski, J. Korchowiec, and Z. Zhou, "Molecular hardness and softness parameters and their use in chemistry," *Int. J. Quantum Chem.* **34**, 349–366 (1988).

⁴⁴M. Franco-Pérez, P. W. Ayers, J. L. Gázquez, and A. Vela, "Thermodynamic responses of electronic systems," *J. Chem. Phys.* **147**, 094105 (2017).

⁴⁵R. G. Parr and P. K. Chattaraj, "Principle of maximum hardness," *J. Am. Chem. Soc.* **113**, 1854–1855 (1991).

⁴⁶G. Hoffmann, F. Guégan, V. Labet, L. Joubert, H. Chermette, C. Morell, and V. Tognetti, "Expanding horizons in conceptual density functional theory: Novel ensembles and descriptors to decipher reactivity patterns," *J. Comput. Chem.* **45**, 1716–1726 (2024).

⁴⁷J. R. Sabin, S. B. Trickey, S. P. Apell, and J. Oddershede, "Molecular shape, capacitance, and chemical hardness," *Int. J. Quantum Chem.* **77**, 358–366 (2000).

⁴⁸A. Holt, G. Karlström, and R. Lindh, "The charge capacity of the chemical bond," *Chem. Phys. Lett.* **436**, 297–301 (2007).

⁴⁹P. Chaquin, "Absolute electronegativity and hardness: An analogy with classical electrostatics suggests an interpretation of the Parr 'electrophilicity index' as a 'global energy index' leading to the 'minimum electrophilicity principle,'" *Chem. Phys. Lett.* **458**, 231–234 (2008).

⁵⁰P. Politzer, "A relationship between the charge capacity and the hardness of neutral atoms and groups," *J. Chem. Phys.* **86**, 1072–1073 (1987).

⁵¹J. C. Watts and J. E. Huheey, "Mulliken electronegativities and ¹¹⁹Sn Mössbauer shifts," *Chem. Phys. Lett.* **14**, 89–90 (1972).

⁵²G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler, "Chemistry with ADF," *J. Comput. Chem.* **22**, 931–967 (2001).

⁵³C. Adamo and V. Barone, "Toward reliable density functional methods without adjustable parameters: The PBE0 model," *J. Chem. Phys.* **110**, 6158–6170 (1999).

⁵⁴E. v. Lenthe, E. J. Baerends, and J. G. Snijders, "Relativistic regular two-component Hamiltonians," *J. Chem. Phys.* **99**, 4597–4610 (1993).

⁵⁵F. Zielinski, V. Tognetti, and L. Joubert, "Condensed descriptors for reactivity: A methodological study," *Chem. Phys. Lett.* **527**, 67–72 (2012).

⁵⁶R. Bader, *Atoms in Molecules: A Quantum Theory, International Series of Monographs on Chemistry* (Clarendon Press, 1990).

⁵⁷B. Wang, P. Geerlings, C. Van Alsenoy, F. Heidar-Zadeh, P. W. Ayers, and F. De Proft, "Investigating the linear response function under approximations following the coupled-perturbed approach for atoms and molecules," *J. Chem. Theory Comput.* **19**, 3223–3236 (2023).

⁵⁸B. A. McGuire, R. A. Loomis, C. M. Charness, J. F. Corby, G. A. Blake, J. M. Hollis, F. J. Lovas, P. R. Jewell, and A. J. Remijan, "Interstellar carbodiimide (HNCNH): A new astronomical detection from the GBT primos survey via maser emission features," *Astrophys. J. Lett.* **758**, L33 (2012).

⁵⁹B. Wang, P. Geerlings, S. Liu, and F. De Proft, "Extending the scope of conceptual density functional theory with second order analytical methodologies," *J. Chem. Theory Comput.* **20**, 1169–1184 (2024).

⁶⁰S. Fias, F. Heidar-Zadeh, P. Geerlings, and P. W. Ayers, "Chemical transferability of functional groups follows from the nearsightedness of electronic matter," *Proc. Natl. Acad. Sci. U. S. A.* **114**, 11633–11638 (2017).

⁶¹N. Sablon, F. D. Proft, P. W. Ayers, and P. Geerlings, "Computing second-order functional derivatives with respect to the external potential," *J. Chem. Theory Comput.* **6**, 3671–3680 (2010).