# An Artificial Neural Network-based Density Functional Approach for Adiabatic Energy Differences in Transition Metal Complexes

João Paulo Almeida de Mendonça, Lorenzo Antonio Mariano, Emilie

Devijver, Noel Jakse, Roberta Poloni

► To cite this version:

João Paulo Almeida de Mendonça, Lorenzo Antonio Mariano, Emilie Devijver, Noel Jakse, Roberta Poloni. An Artificial Neural Network-based Density Functional Approach for Adiabatic Energy Differences in Transition Metal Complexes. Journal of Chemical Theory and Computation, 2023, 19 (21), pp.7555-7566. 10.1021/acs.jctc.3c00600. hal-04247521

HAL Id: hal-04247521

https://hal.science/hal-04247521

Submitted on 18 Oct 2023

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.
An Artificial Neural Network-based Density Functional Approach for
Adiabatic Energy Differences in Transition Metal Complexes

João Paulo Almeida de Mendonça,*,† Lorenzo Antonio Mariano,† Emilie Devijver,‡ Noel Jakse,†
and Roberta Poloni*,†

†Université Grenoble Alpes, CNRS, Grenoble INP, SIMaP, 38000 Grenoble, France

‡Université Grenoble Alpes, CNRS, Grenoble INP, LIG, 38000 Grenoble, France

Received September 26, 2023; E-mail: joao-paulo.almeida-de-mendonca@grenoble-inp.fr; roberta.poloni@grenoble-inp.fr

**Abstract:** During the past decades, approximate Kohn-Sham density-functional theory schemes garnered many successes in computational chemistry and physics; yet the performance in the prediction of spin state energetics is often unsatisfactory. By means of a machine-learning approach, an enhanced exchange and correlation functional is developed to describe adiabatic energy differences in transition metal complexes. The functional is based on the computationally efficient revision of the regularized strongly constrained and appropriately normed (r2SCAN) functional and improved by an artificial neural-network-correction trained over a small dataset of electronic densities, atomization energies and/or spin state energetics. The training process, performed using a bio-inspired non gradient-based approach adapted for this work from Particle Swarm Optimization, is analyzed and discussed extensively. The resulting machine learned meta-GGA functional is shown to outperform most known density functionals in the prediction of adiabatic energy differences for a diverse set of transition metal complexes with varying local coordinations and metal choices.

# 1 Introduction

Density-functional theory (DFT) has revolutionized the field of materials science by providing a powerful tool for predicting the electronic and structural properties of materials,¹⁻³ enabling their design and discovery for a wide range of applications.⁴,⁵ Despite the many successes, the implementations of DFT still face several challenges that limit its predictive power and applicability.⁶ As an example, one can mention the calculation of spin-state energetics in transition metal complexes which represents a challenge for any modern electronic structure ab initio method.⁷⁻¹² This has strong implications in the assistance and guidance alongside experimental efforts towards, for example, the design of novel spin crossover (SCO) complexes. The vast majority of SCO materials, which are of great interest for applications such as molecular spintronics, molecular electronics, and sensors¹³⁻¹⁶ are octahedrally coordinated molecular complexes. For these complexes, the thermodynamics of the spin transition (i.e., the transition temperature T₁/₂) between the low-spin (LS) and high-spin (HS) states is related to the adiabatic energy difference, Δ*E*<sub>H-L</sub>=*E*<sub>HS</sub>-*E*<sub>LS</sub>, i.e., the energy difference between the two spin states computed at their corresponding geometry.

The computational chemistry community has devoted
a substantial effort in understanding the performance of
approximate Khon-Sham density-functional theory (KS-

DFT) schemes.<sup>10,17-21</sup> Large deviations in the values of the spin splitting energies are found among different families of exchange and correlation functionals. Semilocal functionals, such as the generalized gradient approximation (GGA), for example, tend to overstabilize LS,<sup>17,22-26</sup> while Hartree-Fock (HF) overstabilizes HS.<sup>27</sup> In this context, hybrid functionals may be viewed as a possibility to alleviate the problem.<sup>20,27,28</sup> Hubbard *U*-corrected DFT methods, named DFT+U, developed to address the limitations of standard DFT in accurately describing systems with localized electrons,<sup>2,29-31</sup> was also studied in this respect.<sup>32-36</sup> Recently, it was shown that the Hubbard *U* correction in DFT+U may yield large errors in the spin-state energy difference of octahedral Fe-II complexes owing to the bias toward high-spin states imposed by the Hubbard term in the total energy.<sup>37</sup> Later, a non-self consistent density-corrected scheme adopting the Hubbard *U* density was shown to yield a good description of adiabatic energy differences for Fe(II)-based complexes.<sup>38</sup>

Recently, several studies have shown how density func-
tional theory can benefit from machine learning (ML) tech-
niques.<sup>39-50</sup> Pioneer works like the one from Snyder et
al. proved that is possible to learn the kinetic energy
of 1D fermionic systems.<sup>50</sup> Brockherde et al. developed
a ML scheme to learn the density via the external poten-
tial/density Hohenberg-Kohn map, so that self consistency
can be bypassed.<sup>48</sup> Learning the exchange and correlation
functional itself, as demonstrated by the pioneer work by Na-
gai et al.,<sup>44</sup> has also been studied.<sup>39-41,43-45</sup> In this class of
methods, the exchange and correlation is replaced, adjusted,
or corrected by an artificial neural network that receives as
input functions of the electronic density.

These machine learning applications to KS-DFT use the
supervised approach philosophy where the functional "learns
itself" from a given high quality set of data.<sup>49,51</sup> An impor-
tant aspect worth to be mentioned here is that the knowl-
edge of the physical constraints that the exact exchange and
correlation functional must fulfill can be exploited to cre-
ate a strongly transferable functional. In this spirit, Nagai,
Akashi, and Sugino<sup>39</sup> trained their ANN on the atomization
energies and densities from single point calculations of three
molecules, and showed that their resulting metal-GGA func-
tional performed well in describing atomization energies of
molecules not seen in the training, and also on computing
equilibrium lattice constants for solids.

The aim of the present work is to extend upon the latter
approach and train an exchange and correlation functional
that exhibit high accuracy in the prediction of both elec-
tronic densities and adiabatic energy differences of transi-
tion metal complexes, as well as atomization energies and
densities of simpler molecules. An important novelty of this
work, is that the training process is explicitly shown and discussed in detail, thus allowing the reader to understand the optimization process and its limitations. The ML functional is trained using a bio-inspired non gradient-based approach adapted from Particle Swarm Optimization (PSO).⁵² This class of optimization algorithms has shown success in handling intricate nonlinear loss functions⁵³ owing to the collective intelligence ingredient that improves the efficiency of the optimization process. These results show that by training a correction over r2SCAN using three light molecules and three diatomic (metal-nonmetal) transition molecules, the prediction of adiabatic energy differences are improved compared to the state-of-the-art in approximate KS-DFT at no expenses for the performance on atomization energies. Furthermore, the analysis of the non-gradient based training shows that the $\Delta E_{H-L}$ term dominates the optimization and that the electronic density plays an auxiliary role by accelerating the convergence of the loss function. Compared to the known density functional approaches tested here for comparison, the exchange and correlation functional optimized in this work exhibits the best performance for both the Fe complexes and the generality test.

# 2 Methodology

## 2.1 Density Functional Theory and Neural Network Functionals

In the KS-DFT framework, the total energy, $E_{tot}$, and the electronic density, $\rho(\mathbf{r})$, can be computed by solving the KS equations for non-interacting fermions⁵⁴

$$ \left(-\frac{\hbar}{2m} \nabla^2 + v_{\text{eff}}\right) \phi_i = \varepsilon_i \phi_i. \qquad (1) $$

In Eq. 1, the effective potential reads

$$ v_{\text{eff}}(\mathbf{r}) = v_{\text{ext}}(\mathbf{r}) + \int \frac{\rho(\mathbf{r})}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}' + \frac{\delta E_{xc}[\rho(\mathbf{r})]}{\delta \rho(\mathbf{r})}, \quad (2) $$

where $\rho(\mathbf{r}) = \sum_i \phi_i^* \phi_i$, $v_{\text{ext}}(\mathbf{r})$ is an external potential and $E_{xc}$ is the functional of exchange and correlation. The exact form of $E_{xc}$ is unknown and therefore must be approximated. By solving the eigen-equation presented in Eq. 1, one obtains the eigenvalues, $\varepsilon_i$, and eigenfunctions, $\phi_i$, that correspond to the Kohn-Sham energies and orbitals. Because $\rho(\mathbf{r})$ is needed to build $v_{\text{eff}}$ and therefore compute $E_{tot}$, the set of equations above is solved self-consistently, i.e., $\rho(\mathbf{r})$ is recomputed at each iteration until convergence on the total energy and charge density is achieved. The effective potential $v_{\text{eff}}[\rho]$ is at the core of the mean-field KS-DFT since it allows the non-interacting single particle scheme to mimic the electronic correlation of the many electrons system. After solving the KS equations, the total energy can be written as

$$ E = \sum_i \varepsilon_i - \frac{e^2}{2} \iint \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}' d\mathbf{r} + E_{xc}[\rho(\mathbf{r})] \\ - \int \frac{\delta E_{xc}[\rho(\mathbf{r})]}{\delta \rho(\mathbf{r})} \rho(\mathbf{r}) d\mathbf{r}. \qquad (3) $$

The functional of exchange and correlation can be written as

$$ E_{xc} = \int \varepsilon_{xc}[\rho(\mathbf{r})]\rho(\mathbf{r})d\mathbf{r}, \qquad (4) $$

where $\varepsilon_{xc}$ is the exchange and correlation energy density. The idea developed in this work, summarized in Figure 1, is to obtain a $\varepsilon_{xc}$ as a correction over r2SCAN functional,⁵⁵ such that

$$ \varepsilon_{xc} = F_x \varepsilon_x^{\text{r2SCAN}} + F_c \varepsilon_c^{\text{r2SCAN}}. \qquad (5) $$

Recently, Nagai and co-workers⁴⁴ developed a new functional, defined as a correction over the strongly constrained and appropriately normed (SCAN),⁵⁶ trained on energies and densities of a small subset (three molecules only made up by first and second row elements) of the G2 dataset.⁵⁷ The result showed an improvement over SCAN in the prediction of ionization potentials and atomization energies on the complete G2 dataset. Yet, previous reports on the literature pointed out some major numerical instabilities in some meta-GGA functionals⁵⁸ including SCAN itself, for electron densities rapidly changing in space. These manifest as sharp features of the energy derivatives and make the functional unstable over grid changes thus resulting in bad gradients estimation during the self-consistent field process.⁵⁹,⁶⁰ To solve this problem, Bartok and co-workers developed a regularized version of SCAN (rSCAN) which exhibits improved smoothness properties and improves its numerical performance.⁶¹ However, this improvement came at the expense of some of the physical limits originally obeyed by SCAN. r2SCAN⁵⁵ follows this implementation path, recovering most of the physical limits that SCAN originally had. The 17 exact conditions applicable to meta-GGA (listed in the SI from Sun et al.⁵⁶) were imposed in the original definition of SCAN and most of these conditions are present on r2SCAN. Only the fourth-order gradient expansion condition, lost in the rSCAN implementation, is not recovered.⁵⁵

In this work, in order to address the case of transition metal complexes, where rapid fluctuations of the electronic density are expected, a new ANN-based functional defined as a correction over r2SCAN is proposed. This choice allows us to circumvent the instabilities of SCAN that manifest themselves in large changes of the energy upon tiny changes of the density, as shown in Figure S1 for the FeO molecule. Such instability (observed for FeO and not for H₂O) is found when training a correction over SCAN. Finally, we impose that $\varepsilon_{xc}$, defined in eq. 5, satisfies the physical limits imposed in r2SCAN using Lagrange multipliers, as done by Nagai and co-workers in their more recent work.³⁹

Figure 2 shows a diagrammatic representation of the feed-forward neural network adopted here. The architecture reflects the choice to adopt a small number of fitted parameters. Recent efforts on ML functionals are moving in the direction of complex and deeper neural networks which results in a large amount of parameters to be fitted. For example, Nagai et al. use an architecture with hidden layers composed of 100 neurons each, resulting in a total number of parameters surpassing $10^4$.³⁹ The recent DM21 functional from DeepMind also report a total of roughly $4 \times 10^5$ parameters.⁴¹ Here, the choice was made to move in the opposite direction, adopting a smaller neural network with precisely 1506 parameters, which is more efficient to train.

The local and semilocal functions that are provided to the input neurons of the ANN are consistent with a meta-GGA functional. Evaluated for each point in the integration grid,
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Schematic representation of a single point calculation performed in the present approach. The step where the ANN correction is explicitly applied over r2SCAN is highlighted in blue.</figcaption></figure>

the inputs are defined as

$$ 
\begin{aligned}
\rho &:= \sum_{i}^{\text{occ}} |\phi_i|^2, \\
\zeta &:= \frac{\rho_{\uparrow} - \rho_{\downarrow}}{\rho}, \\
s &:= \frac{|\nabla \rho|}{2\sqrt[3]{3\pi\rho^4}}, \\
\tau &:= \frac{1}{2} \sum_{i}^{\text{occ}} |\nabla \phi_i|^2.
\end{aligned} 
\quad (6) 
$$

These quantities correspond to the electronic density, the spin density, the scaled density gradient, and the kinetic-energy density associated to the electronic state of the system, respectively. Pre-processing of these functions is performed with a tanh based filter before they enter the ANN.<sup>39,44</sup> These filters are built in such a way to fit the possible values of each quantity inside the (-1, 1) interval, also mimicking some of the LDA energy features for the $\rho$ and $\zeta$ normalization and translating the kinetic energy of the uniform electron gas to 1. Also, this filtering allows us to handle infinities with finite values in the Lagrange multipliers, e.g., for $\rho \to \infty$ we have $\bar{n}_s \to 1$. Those transformations are defined as: $\bar{n}_s = \tanh\left(\left(\frac{\rho}{\rho_c}\right)^{\frac{1}{3}}\right)$, $\bar{\zeta} = \tanh\left(\frac{1}{2}\left\{(1+\zeta)^{\frac{4}{3}} - (1-\zeta)^{\frac{4}{3}}\right\}\right)$, $\bar{s} = \tanh(s)$, and $\bar{\tau} = \tanh\left(\frac{\tau-\tau_{\text{unif}}}{\tau_{\text{unif}}}\right)$; where $\tau_{\text{unif}} = \frac{3}{10}(3\pi^2)^{\frac{2}{3}}\left(\frac{\rho}{\rho_c}\right)^{\frac{5}{3}}$.

The corrections $F_x$ and $F_c$ are computed separately so that different asymptotic limits are applied to each of them.<sup>39,56</sup> For the correct uniform coordinate density-scaling condition<sup>62</sup> and the exact spin scaling relation<sup>63</sup> to be satisfied simultaneously, the exchange energy should not depend explicitly in $\rho$ and $\zeta$, so this values are not included as inputs in the calculation of $F_x$. The weights $w_{3c}$ are set as a 32 × 16 matrix, using the second hidden layer from both the exchange and the correlation networks to compute the third hidden layer of the correlation network. This improves the training process<sup>39</sup> and is only possible in this direction, since the inputs for $F_x$ are a part of the inputs for $F_c$. The activation function used here is a version of the softplus function modulated so it has value and derivative as 1 at the origin. In this way, if all parameters on the NN where set to zero, the functional obtained would be exactly r2SCAN.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Architecture of the artificial neural network used in the present work. The w's and b's represent the weights and bias used on the feed-forward process and the products in parenthesis express the number of parameters in each of those elements.</figcaption></figure>

## 2.2 Training of the Neural Network Functional

### 2.2.1 Loss Function

The loss function is the sum of errors of adiabatic energy differences, atomization energies, and electron densities

$$ 
\begin{aligned}
\Delta_{\text{err}} = & \sum_{i}^{G2} \left[ \frac{1}{E_0} |AE_i^{\text{REF.}} - AE_i^{\text{DFT}}| + \right. \\
& \left. + \frac{1}{N_i} \int |\rho_i^{\text{REF.}} - \rho_i^{\text{DFT}}| dV \right] + \\
& \sum_{j}^{\text{TMC}} \left[ \frac{1}{E_0} |\Delta E_{\text{H-L},j}^{\text{REF.}} - \Delta E_{\text{H-L},j}^{\text{DFT}}| + \right. \\
& \left. \frac{1}{N_j} \int \left( |\rho_{\text{H},j}^{\text{REF.}} - \rho_{\text{H},j}^{\text{DFT}}| + |\rho_{\text{L},j}^{\text{REF.}} - \rho_{\text{L},j}^{\text{DFT}}| \right) dV \right],
\end{aligned} 
\quad (7) 
$$
where the first summation is performed over three small molecules (NO, H₂O, NH₃) taken from the original G2 set with experimental reference values of atomization energies, AE<sup>REF.</sup>, extracted from Curtiss et al.⁵⁷ To avoid confusion, from now on italic G2 refers to this set of three molecules, while the (non-italic) G2 refers to the complete set by Curtiss et al.⁵⁷

The second summation is performed over three diatomic transition metal complexes (TMC) with reference values for adiabatic energy difference, ΔE<sup>REF.</sup><sub>H-L</sub>, taken from experiments<sup>64</sup> (and close to the CCSD(T) values within 0.02 eV). These molecules are FeO, CuF, CrH (see upper panel of Figure 3). The error in the density for both the G2 set and the TMC is computed using the L<sup>1</sup> metric, by taking the CCSD(T) reference for the electronic density. Finally, E<sub>0</sub> = 1 Hartree, and N<sub>j</sub> is the number of electrons in the system j, both used to make the loss function dimensionless.

In order to evaluate the effects of the database, six different training sets were considered: G2 ([NO + H₂O + NH₃]), G2+FeO, G2+CuF, G2+CrH, TMC ([FeO + CuF + CrH]), and ALL (G2 + TMC). The terms of the loss function only apply when at least one of the related molecules is present in the training set. The Spearman correlation matrices and the principal component analysis of the loss function terms during the optimization are shown in the SI (Figures S2-S19).

The inclusion of the first term in eq. 7 was shown to yield good results for energetics of light element-molecules. 43,44 Besides, it was shown that training sets with different types of data per sample (i.e., sparse training data) can significantly improve the performance of a machine learned-density functional. 65

The experimental values of ΔE<sub>H-L</sub> are 1.81 eV for CuF (spin change: S = 1 → 3),<sup>66,67</sup> 0.14 eV for FeO (S = 5 → 7),<sup>68</sup> and −1.39 eV for CrH (S = 4 → 6).<sup>69,70</sup> The choice of these specific molecules is done after considering four points: (i) the values of ΔE<sub>H-L</sub> span a large range from negative to positive values, (ii) the diversity in terms of composition, (iii) the adiabatic energy difference involves different spin states in each case, and (iv) the configuration interaction weight of the dominant electronic configuration computed using CASSCF is different in the three cases pointing to different degrees of static correlations (C<sup>2</sup><sub>0</sub> = 0.810, 0.783, and 0.938; for the low spin configuration of CuF, FeO, and CrH respectively).<sup>71</sup> We do not claim here that our approach is capable to treat systems with multiconfigurational character. Rather, that the ANN-DFT yields a solution similar to the coupled cluster one, which gives an excellent agreement with experiment for these molecules.<sup>64</sup>

### 2.2.2 Optimization Algorithm

Since the equilibrium energies and densities are computed after converging the SCF process with a given set of parameters, the loss function depends on quantities that are related in a non-trivial way to the parameters of the ANN-functional. Gradient based optimizations of loss functions that depend on quantities computed upon solving the SCF cycle have been performed in the past by adopting symmetry arguments in order to minimize the grid size and/or by fixing a small number of SCF steps.<sup>43,46,65,72</sup> When these conditions do not apply, such as in this case, non-gradient based approaches have been employed.<sup>39,44</sup> Besides, recent studies<sup>53</sup> have shown that smart non-gradient based optimization techniques can outperform gradient-based ones for

complicated loss functions, where this is not linearly or even continuously dependent on the parameters. This is exactly the nature of loss function here. To handle this complex parameters-loss function landscape, an adaptation of Particle Swarm Optimization (PSO)<sup>52</sup> is used as the training algorithm. PSO is built over a bio-inspired meta-heuristic: solutions (points on the parameters space of our neural network) will be treated as particles in a swarm (a collection of solutions) that forage for resources (that search for the minima of the loss function). The application of the PSO algorithm to the problem of functional optimization is shown in the flowchart of Figure 4(a). In the initial step, a search region in the parameters space is defined, and a swarm of particles with positions x<sub>i</sub>(t = 0) is randomly generated. These are initialized with velocities v<sub>i</sub>(t = 0) = 0 in this space. At each step t in the PSO evolution, all the particles are submitted to a collective loss function evaluation. From this, the best solution ever visited by the swarm x<sub>i,best</sub>(t) and the best solution ever visited by the each particle x<sub>i,best</sub>(t) are saved. The velocities and positions on the swarm at t+1

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Representation of the different sets of molecular complexes studied in this work.</figcaption></figure>
are then updated as follows:

$$ \mathbf{v}_i(t+1) = \omega \mathbf{v}_i(t) + c_1 \mathrm{Rand}(0, 1)[(\mathbf{x}_{i,best}(t) - \mathbf{x}_i(t))] \quad (8) $$

$$ + c_2 \mathrm{Rand}(0, 1)[(\mathbf{x}_{best}(t) - \mathbf{x}_i(t))], $$

$$ \mathbf{x}_i(t+1) = \mathbf{x}_i(t) + \mathbf{v}_i(t+1), \quad (9) $$

where Rand(0, 1] are random numbers independently generated in the (0, 1] interval and $\omega$, $c_1$, and $c_2$ are constants that ponder the different behaviors of the particles on the swarm: $\omega$ is the weight for the term that mimics the inertia of a flying particle and is associated with the pink component of the velocity in Figure 4(b). $c_1$ gives to the particle a velocity component that points toward the best solution known by the particle itself and it is the red component in Figure 4(b). $c_2$ gives a velocity component that points toward the best solution known by the swarm, which implements the concept of collective intelligence and it is shown by the blue component in Figure 4(b).

The standard PSO approach⁵² is modified here by changing the search region between PSO runs as described below. This solves the observed instability when trying to evaluate the loss function in points that are too far from solutions that we know already. At the beginning of the optimization process, the particles are initialized inside a box of size $s$ around the origin ($\epsilon_{x/c} = F_{x/c}(0) \times \epsilon_{x/c}^{2SCAN} = \epsilon_{x/c}^{2SCAN}$). Once a given number of PSO steps is done, the searching box is moved so that its center corresponds to the best solution found until then as shown in Figure 4(c). The particles are redistributed and the PSO search restarts. The PSO steps are performed using $c_1 = 0.1$, $c_2 = 0.2$, and $\omega = 0.1$, that satisfies the known criteria for stability on PSO.⁷³ The search box in each migration is defined using $s = 0.05$ (with periodic boundary conditions) and 10 PSO steps are taken in each migration. A total of 32 particles are used in each step and a total of 50 migrations is performed for every run.

## 2.3 Technical Details

Our functional is implemented using an adaptation of the xc_pcNN import provided by Nagai and co-workers:⁷⁴ the neural network architecture is changed to the one proposed here and the physical constraints and corrections were adapted to the case of r2SCAN (imported from libxc⁷⁵).

The training code is an in-house application developed in PYTHON3, combining the functionalities of pyscf⁷⁶ and ase⁷⁷ to compute the DFT SCFs and pyswarm⁷⁸ to perform the optimization inside each migration during the training process. For each particle (each parametrization of the ANN-DFT) the full SCF cycle is performed in order to compute the loss function. The 6-311++G(3df,3pd) and the cc-pCVTZ basis sets are used for the G2 and the transition metal complexes, respectively. The number of iterations within the SCF is set to 50 and the convergence criteria on the energy change is set to $10^{-6} E_h$. The calculations that failed to converge within this criterion were not considered during the optimization process.

The geometries of the molecules in the training set are taken as follows: for the G2 set these are taken from ase.collections.g2 (originally obtained via G2(MP2,SVP) by Curtiss et al.⁵⁷) while for the transition metal complexes the MRCI geometries are considered.³²

The reference densities are obtained via orbital optimized-

CCSD(T), as implemented in ORCA 5.0⁷⁹ using 6-311++G(3df,3pd) and cc-pVTZ basis sets for the G2 set and the transition metal complexes, respectively. All DFT calculations reported in Table 1 (PBE, SCAN, B3LYP, r2SCAN) and Table 2 (TPSS, r2SCAN, PBE) are performed using pyscf with cc-pVTZ basis sets, except for the PBE[U] results. The calculations performed using NAS⁷⁴ and DM21⁸⁰ functionals use the codes provided by their authors, and both are adaptations of pyscf. The PBE[U] values are obtained using the Hubbard U-density corrected scheme in Quantum ESPRESSO.⁸¹ We use ultrasoft GBRV pseudopotentials⁸² with wavefunction and charge density cutoffs of 80 Ry and 800 Ry, respectively. The Markov-Payne correction⁸³ is used to treat charged systems under periodic boundary conditions. Calculations in Table 1 are taken from ref.³⁸ Calculations in Table 2 are performed using the same computational scheme, meaning that the PBE total energy is evaluated using the PBE+U electronic density, and the U is the self-consistent linear response $U_{sc}$ computed separately for each molecule in each spin state.³⁸ The $U_{sc}$ parameters are obtained through an iterative process using the hp.x routine. For FeL2 the $U_{sc}$ are 8.92 eV (2S+1=1), 9.18 eV (2S+1=3), and 8.81 eV (2S+1=5); for MnL2 6.86 eV (2S+1=2), 6.74 eV (2S+1=4), and 4.80 eV (2S+1=6); for FeL2SH 7.38 eV (2S+1=4), and 6.25 eV (2S+1=6); for [Co(NCH)₆]²⁺ 7.83 eV (2S+1=2), and 6.76 eV (2S+1=4); for NiCp(acac) 9.66 eV (2S+1=1), and 7.32 eV (2S+1=3); for MnCp₂ 6.93 eV (2S+1=2), and 4.60 eV (2S+1=6). The geometries reported by Pierloot et al.⁸⁴ are used for all the calculations in Table 2 in order to maintain consistency with the CASPT2/CC reference values. The calculations performed to evaluate the performance of different functionals on the atomization energies of the G2 set (shown in section S3 the SI) use geometries form ase.collections.g2 and the 6-311++G(3df,3pd) basis set.

# 3 Results

## 3.1 Trainings

In Figure 5(a), we draw the evolution of the loss function as defined in eq. 7, where only the best solution (particle) at each migration is reported. After 30 migrations, the solution is reaches a plateau for every run except for G2. For this training set, a monotonic decrease of the loss function is observed until 36 migrations. A different behavior is observed for the training sets including TMC, where the slow monotonic decrease is followed by a sharp decrease, and then by a stationary region. This change in behavior is associated with solutions that cannot be improved via migrations. At this point we declare the solution has converged. The migrations selected as converged for each training set are: #23 for ALL, #20 for TMC, #19 for G2+FeO, #18 for G2+CuF, #14 for G2+CrH, and #36 for G2.

The training for the G2 follows path of error reduction of the atomization energy at the expense of a small increase in the error of the density (see Figures S2-S3). The loss function stops decreasing after migration 36 and this is associated to a limit in the error of the atomization energy for NO and H₂O. When the transition metal complexes are included in the training, a fast optimization of the adiabatic energy difference is observed simultaneously with the corresponding
<figure><img src="image_4.png" /><figcaption>**Figure 4.** (a) Flowchart for the functional training approach as described in the text, including the PSO loop (green) and the migrations scheme (red). (b) Schematic representation of a PSO step. Magenta, red and blue vectors are associated with the right-hand side terms of equation 8, respectively. (c) Migration of the search region: all particles (grey and blue circles) are removed from the previous box (grey square) and reinitialized randomly inside a new one (blue square).</figcaption></figure>

densities, highlighting the role of accurate electronic densities in the evaluation of adiabatic energy differences.²⁰,³⁸ This is seen in the Spearman correlation matrices computed during the optimization for the first 2000 functionals visited for the G2+FeO training (see Figure S6) and for the ALL training for FeO and CuF (see Figure S15). In these training runs (when one or more transition metal complexes are included), the stationary region is associated with one of the $\Delta E_{H-L}$ reaching limit accuracy (i.e., the error associated to this quantity reaches zero). After this point, the process searches ways to optimize the other degrees of freedom (densities or atomization energies) and this leads to an increased dispersion on the error values associated with the $\Delta E_{H-L}$.

Figures 5(c) and (d) report the heatmap for the Spearman correlation matrix between the partial errors that compose the loss function for particles sampled at migration #15, i.e., during the optimization, and at migration #22, where the optimization stops, for the ALL training ({NO, H₂O, NH₃, FeO, CuF, CrH}). At migration #15, in Figure 5(c) we see that the errors on energies and densities of FeO and CuF, as well as the densities on the G2 molecules, are optimized simultaneously. The loss function part associated with $\Delta E_{H-L}^{CrH}$ and $\rho_{CrH}$ is negatively correlated with those values meaning that their error increase while the total loss function decreases.

In the case of the ALL and TMC training sets, $\Delta E_{H-L}^{CuF}$ was the first to reach a converged value close to zero. At this point of the optimization, the particles dispersion was enough to allow to find a new optimizing path for $\Delta E_{H-L}^{FeO}$, as seen in the top left green block in Figure 5(d) (positive correlation between these quantities and the total loss function).

The second block in Figure 5(d), showing strong positive, correlation is negatively correlated with $\Delta_{err}$.

## 3.2 Performance Tests

### 3.2.1 The performance on Fe complexes

A test set of six Fe complexes, previously studied by our group,³⁷,³⁸ is chosen here to assess the performance of the six ANN functionals. This set of molecules is shown in the middle panel of Figure 3. These molecules exhibit a octahedrally coordinated Fe with the associated ligand field strength varying such that the reference spin splitting energies change from negative to positive values (see Table 1).

An additional functional, named here $ANN_{bias}^{ALL}$, is tested. Starting from the collection of functionals generated during the training performed using the ALL training set (see previous section), we select the one that gives the lowest mean absolute errors (MAE) of the adiabatic energy differences computed for the six molecules with respect to the CASPT2/CC³⁸ energies. For each migration computed using the ALL training set, we take the best solution in the loss function to compute this new metric. The evolution of this MAE metric along the training is shown in the magenta curve of Figure 5(b). In the same figure the evolution of the loss function for the ALL training set is reported in black for comparison (this is the same curve as in Figure 5(a)). The functional that best performs on the loss function is called $ANN_{lf}^{ALL}$. We stress here that the same set of parameters
<figure><img src="image_5.png" /><figcaption>**Figure 5.** (a) Evolution of the loss function values during the training process, for the six different training sets considered here. (b) Comparison between the evolution of the loss function and that of the performance in the Fe complexes test along the migrations in the training with the complete ALL training set. (c) and (d) show Spearman's rank correlation coefficients between the total loss function and each one of the L¹ distances that compose it. The two correspond, respectively, to the samples visited along migrations #15 and #22. The correlation matrix has converted in a heatmap.</figcaption></figure>

is used for the magenta and the black curves of Figure 5(b), at each migration.

The error on the performance test of the six octahedral molecules reaches a minimum at step #18. The functional that performs the best on this test is called $ANN_{bias}^{ALL}$, since, instead of considering only the minimization of the loss function, we are biasing it to best perform for the case of Fe complexes. This would be similar to an early stop procedure, validated by the Fe complexes set. This approach is not a traditional one, and it is introduced here to study the impact of such a choice. We emphasize that the more traditional approach of separating part of the training data aside to be used as test is not used here, given the small size of the training data. Table 1 gathers the $ΔE$ computed using the seven ANN-functionals (from the six training sets along with the biased one) for this set of Fe complexes. The performance of other DFT functionals, together with two recent ANN-based functionals, are reported for comparison in Table 1.

The two ANN-based functions are the one from ref. 39 named hereafter NAS from Nagai-Akashi-Sugino, and the DM21.⁴¹

The performance of PBE, PBE[U], SCAN, and B3LYP was already discussed in previous studies.²⁰,³⁸ Here, consistent with previous results, MAE of 0.92, 0.79 and 0.70 eV are found. The values computed using the Hubbard U-density corrected scheme using with linear-response values of U (named PBE[U]) are also reported for comparison. These are taken from ref. ³⁸ where the method was shown to yield the smallest error among 19 different density functional approaches. The PBE[U] associated MAE is 0.44 eV, and larger errors are observed for TPSSh (MAE=0.49 eV), PBE0 (MAE=0.64 eV) and M06-2X (MAE=2.41 eV).

It is worth noting here that r2SCAN (MAE=0.49 eV) yields a significant improvement over SCAN (MAE=0.79 eV) in terms of $ΔE$, with a performance close to PBE[U]. This can be related to the improved numerical stability of the
<table><caption>**Table 2.** Performance of different exchange and correlation functionals for the prediction of adiabatic energy differences for the test set. Some spin configurations could not be converged and thus a few values from the table are missing, and marked as "-".</caption><thead><tr><th rowspan="2"></th><th rowspan="2">2S+1</th><th colspan="7">ΔE<sub>H-L</sub> (eV)</th></tr><tr><th>ANN<sub>lf</sub><sup>ALL</sup></th><th>ANN<sub>bias</sub><sup>ALL</sup></th><th>r2SCAN</th><th>TPSS</th><th>PBE</th><th>PBE[U]</th><th>DM21</th><th>REF.<sup>b</sup></th></tr></thead><tbody><tr><td rowspan="2">FeL<sub>2</sub></td><td>1 → 5</td><td>-1.439</td><td>-0.944</td><td>-0.611</td><td>0.227</td><td>-</td><td>-1.096</td><td>-1.968</td><td>-1.487</td></tr><tr><td>3 → 5</td><td>0.193</td><td>0.513</td><td>0.746</td><td>1.259</td><td>1.245</td><td>0.756</td><td>0.814</td><td>0.213</td></tr><tr><td rowspan="2">MnL<sub>2</sub></td><td>2 → 6</td><td>-2.078</td><td>-1.427</td><td>-0.957</td><td>-0.055</td><td>-0.111</td><td>-0.961</td><td>-</td><td>-1.782</td></tr><tr><td>4 → 6</td><td>-0.061</td><td>0.207</td><td>0.414</td><td>0.835</td><td>0.919</td><td>0.491</td><td>-0.739</td><td>-0.455</td></tr><tr><td rowspan="2">FeL<sub>2</sub>SH</td><td>2 → 6</td><td>-0.464</td><td>0.081</td><td>0.470</td><td>-0.416</td><td>1.087</td><td>-</td><td>-</td><td>0.399</td></tr><tr><td>4 → 6</td><td>-0.488</td><td>-0.130</td><td>0.123</td><td>0.588</td><td>0.521</td><td>0.245</td><td>-</td><td>-0.017</td></tr><tr><td rowspan="2">[Co(NCH)<sub>6</sub>]<sup>2+</sup></td><td>2 → 4</td><td>-0.372</td><td>-0.170</td><td>-0.027</td><td>0.31</td><td>0.476</td><td>0.165</td><td>0.057</td><td>-0.581</td></tr><tr><td>1 → 3</td><td>-0.204</td><td>-0.177</td><td>-0.157</td><td>0.046</td><td>0.163</td><td>0.095</td><td>0.414</td><td>0.117</td></tr><tr><td>MnCp<sub>2</sub></td><td>2 → 6</td><td>-0.220</td><td>-0.354</td><td>0.184</td><td>0.649</td><td>0.981</td><td>0.432</td><td>-</td><td>0.304</td></tr><tr><td>MAE<sup>b</sup></td><td></td><td>0.349</td><td>0.406</td><td>0.474</td><td>0.945</td><td>0.885<sup>a</sup></td><td>0.482<sup>a</sup></td><td>0.460<sup>a</sup></td><td></td></tr></tbody></table>

case of the FeL₂SH molecule. The deviation from the reference is 0.9 eV for the 2 → 6 transition and 0.5 eV for the 4 → 6 transition, while it is only approximately 0.1 eV for r2SCAN (for either transitions). We attribute this issue to the challenging contribution from the sulfur atom also not seen during the training. When our ANN<sub>lf</sub><sup>ALL</sup> functional is tested on a similar molecule without S, i.e. FeL₂OH, the deviation from the CASPT2/CC decreases to 0.3 eV and 0.1 eV, respectively.

Finally, to assess how the choice of the geometry affects the ANN<sub>lf</sub><sup>ALL</sup> values, we compute the adiabatic energy differences of Table 2 by employing geometries optimized using r2SCAN (see results in Table S2). The values change by up to 0.5 eV for FeL2. Yet, the overall performance (MAE) only marginally changes with respect to the geometries from Pierloot and coworkers.⁸⁴ The MAE is 0.401 eV for the r2SCAN geometries versus 0.349 eV in Table 2.

## 4 Conclusion

The ANN-based exchange and correlation functional in DFT developed in this work outperforms most known local, semilocal, and meta-GGA functionals in the prediction of adiabatic energy differences for transition metal complexes. It is built as a physically-constrained meta-GGA artificial neural-network correction to r2SCAN and trained over a small dataset of electronic densities, atomization energies and/or adiabatic energy differences. Our results show a good performance of the non gradient-based bioinspired training algorithm PSO where the optimization can reach a hard limit, as in the case of the adiabatic energy differences of CuF, FeO, and CrH. Six functionals were obtained by using different subsets of the training molecules. The functional optimized by using the whole training set, ANN<sub>lf</sub><sup>ALL</sup> (MAE=0.41 eV), and that obtained using the G2 + CuF, ANN<sup>G2+CrF</sup> (MAE=0.42 eV), yield a performance on the Fe complexes test superior to the best performer r2SCAN (MAE=0.44 eV). Such a test performed on octahedral iron complexes shows that our functional substantially improves over recent ANN-based functionals, such as the meta-GGA NAS and the range-separated hybrid DM21. The best performer on the loss function, ANN<sub>lf</sub><sup>ALL</sup>, best performs also

on the generality test which includes molecules exhibiting diverse local coordinations, ligand compositions, and transition metal choices (MAE=0.349 eV). In addition, we show that by changing the criteria used to select the functional (from the loss function itself to the error in a test set), we were able to bias the results to improve the performance on a specific set of data. For example, the best performer on the Fe complexes test, named ANN<sub>bias</sub><sup>ALL</sup>, yields a MAE as small as 0.25 eV, but the associated MAE on the generality test is larger (MAE=0.41 eV) than for the ANN<sub>lf</sub><sup>ALL</sup>. Overall, the results of this study highlight the impact of machine learning approaches to enhance the accuracy of approximate density functional theory through the development of optimized functionals.

**Acknowledgement** Computational resources were provided by the CINES and IDRIS under Project No. INP2227/72914, as well as CIMENT/GRICAD for computational resources. This work was performed within the framework of the Centre of Excellence of Multifunctional Architecture Materials CEMAM-ANR-10-LABX-44-01 funded by the "Investments for the Future" Program. This work has been partially supported by MIAI@Grenoble Alpes (ANR-19-P3IA-0003). The authors thank Lucia Reinig for fruitful discussions. Discussions within the French collaborative network in artificial intelligence in materials science GDR CNRS 2123 (IAMAT) are also acknowledged.

## Supporting Information Available

Preliminary tests performed using SCAN (S1); correlation analysis for the trained functionals (S2); generality test on atomization energies (S3); analysis of the optimized densities of the transition metal molecules (S4); detailed flowchart of the work (S5); performance test on ANN<sub>lf</sub><sup>ALL</sup> using different sets of geometries. This information is available free of charge via the Internet at http://pubs.acs.org

5

References

(1) Poloni, R.; Machon, D.; Fernandez-Serra, M. V.; Le Floch, S.; Pascarelli, S.; Montagnac, G.; Cardon, H.; San-Miguel, A. High-

<footnote>a Average computed over available values.</footnote>

<footnote>b Reference values come from CASPT2/CC.⁸⁴</footnote>
pressure stability of Cs$_6$C$_60$. Phys. Rev. B **2008**, 77, 125413.

(2) Poloni, R.; Mariano, A. L.; Prendergast, D.; Garcia-Barricanal, J. Probing the electric field-induced doping mechanism in YBa$_2$Cu$_3$O$_7$ using computed Cu K-edge x-ray absorption spectra. J. Chem. Phys. **2018**, 149, 234706.

(3) Yang, C.-T.; Kshirsagar, A. R.; Eddin, A. C.; Lin, L.-C.; Poloni, R. Tuning Gas Adsorption by Metal Node Blocking in Photoresponsive Metal-Organic Frameworks. Chem. Eur. J. **2018**, 24, 15167-15172.

(4) Jones, R. O. Density functional theory: Its origins, rise to prominence, and future. Rev. Mod. Phys. **2015**, 87, 897-923.

(5) Jain, A.; Shin, Y.; Persson, K. A. Computational predictions of energy materials using density functional theory. Nat. Rev. Mater. **2016**, 1, 15004.

(6) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Challenges for Density Functional Theory. Chem. Rev. **2012**, 112, 289-320, PMID: 22191548.

(7) Wilbraham, L.; Verma, P.; Truhlar, D. G.; Gagliardi, L.; Ciofini, I. Multiconfiguration Pair-Density Functional Theory Predicts Spin-State Ordering in Iron Complexes with the same Accuracy as Complete Active Space Second-Order Perturbation Theory at a Significantly Reduced Computational Cost. J. Phys. Chem. Lett. **2017**, 8, 2026-2030.

(8) Domingo, A.; Ángels Carvajal, M.; De Graaf, C. Spin Crossover in Fe(II) Complexes: An Ab Initio Study of Ligand σ-Donation. Int. J. Quantum Chem. **2010**, 110, 331-337.

(9) Radoń, M. Benchmarking Quantum Chemistry Methods for Spin-State Energetics of Iron Complexes Against Quantitative Experimental Data. Phys. Chem. Chem. Phys. **2019**, 21, 4854-4870.

(10) Swart, M. Accurate Spin-State Energies for Iron Complexes. J. Chem. Theory Comput. **2008**, 4, 2057-2066.

(11) Phung, Q. M.; Feldt, M.; Harvey, J. N.; Pierloot, K. Toward highly Accurate Spin State Energetics in First-Row Transition Metal Complexes: A Combined CASPT2/CC Approach. J. Chem. Theory Comput. **2018**, 14, 2446-2455.

(12) Reimann, M.; Kaupp, M. Spin-State Splittings in 3d Transition-Metal Complexes Revisited: Toward a Reliable Theory Benchmark. J. Chem. Theory Comput. **2023**, 19, 97-108, PMID: 36576816.

(13) Molnár, G.; Mikolasek, M.; Ridier, K.; Fahs, A.; Nicolazzi, W.; Bousseksou, A. Molecular Spin Crossover Materials: Review of the Lattice Dynamical Properties. Ann. Phys. **2019**, 531, 1900076.

(14) Kumar, K. S.; Ruben, M. Emerging Trends in Spin Crossover (SCO) Based Functional Materials and Devices. Coordin. Chem. Rev. **2017**, 346, 176-205.

(15) Resines-Urien, E.; Burzuri, E.; Fernandez-Bartolome, E.; García García-Tuñón, M. A.; de la Presa, P.; Poloni, R.; Teat, S. J.; Costa, J. S. A switchable iron-based coordination polymer toward reversible acetonitrile electro-optical readout. Chem. Sci. **2019**, 10, 6612-6616.

(16) Resines-Urien, E.; Piñeiro-López, L.; Fernandez-Bartolome, E.; Gamonal, A.; Garcia-Hernandez, M.; Sánchez Costa, J. Cova- lent post-synthetic modification of switchable iron-based coordination polymers by volatile organic compounds: a versatile strategy for selective sensor development. Dalton Trans. **2020**, 49, 7315-7318.

(17) Swart, M.; Groenhof, A. R.; Ehlers, A. W.; Lammertsma, K. Validation of Exchange-Correlation Functional s for Spin States of Iron Complexes. J. Phys. Chem. A **2004**, 108, 5479-5483.

(18) Pierloot, K.; Vancoillie, S. Relative Energy of the High-($^5T_{2g}$) and Low-($^1A_{1g}$) Spin States of $[Fe(H_2O)^6]^{2+}$, $[Fe(NH_3)^6]^{2+}$, and $[Fe(bpy)3]^{2+}$: CASPT2 Versus Density Functional Theory. J. Chem. Phys. **2006**, 125, 124303.

(19) Droghetti, A.; Alfé, D.; Sanvito, S. Assessment of Density Functional Theory for Iron(II) Molecules Across the Spin-Crossover Transition. J. Chem. Phys. **2012**, 137, 124303.

(20) Song, S.; Kim, M.-C.; Sim, E.; Benali, A.; Heinonen, O.; Burke, K. Benchmarks and Reliable DFT Results for Spin Gaps of Small Ligand Fe(II) Complexes. J. Chem. Theory Comput. **2018**, 14, 2304-2311.

(21) Daku, L. M. L.; Aquilante, F.; Robinson, T. W.; Hauser, A. Accurate Spin-State Energetics of Transition Metal Complexes. 1. CCSD(T), CASPT2, and DFT Study of $[M(NCH)6]^{2+}$ (M = Fe, Co). J. Chem. Theory Comput. **2012**, 8, 4216-4231.

(22) Fouqueau, A.; Casida, M. E.; Daku, L. M. L.; Hauser, A.; Neese, F. Comparison of Density Functionals for Energy and Structural Differences Between the High-$^5T_{2g}$: $(t_{2g})^4(e_g)^2$ and Low-$^1A_{1g}$: $(t_{2g})^6(e_g)^0$ Spin States of Iron(II) Coordination Compounds. II. More Functionals and the Hexaminoferrous Cation, $[Fe(NH_3)^6]^{2+}$. J. Chem. Phys. **2005**, 122, 044110.

(23) Radoń, M. Revisiting the Role of Exact Exchange in DFT Spin-State Energetics of Transition Metal Complexes. Phys. Chem. Chem. Phys. **2014**, 16, 14479-14488.

(24) Mortensen, S. R.; Kepp, K. P. Spin Propensities of Octahedral Complexes from Density Functional Theory. J. Phys. Chem. A **2015**, 119, 4041-4050.

(25) Ioannidis, E. I.; Kulik, H. J. Towards Quantifying the Role of Exact Exchange in Predictions of Transition Metal Complex Properties. J. Chem. Phys. **2015**, 143, 034104.

(26) Kulik, H. J. Perspective: Treating Electron Over-Delocalization with the DFT+U Method. J. Chem. Phys. **2015**, 142, 240901.

(27) Reifer, M.; Salomon, O.; Hess, B. A. Reparameterization of hybrid functionals based on energy differences of states of different multiplicity. Theor. Chem. Acc. **2001**, 107, 48-55.

(28) Prokopiou, G.; Kronik, L. Spin-State Energetics of Fe Complexes from an Optimally Tuned Range-Separated Hybrid Functional. Chem.-Eur. J. **2018**, 24, 5173-5182.

(29) Anisimov, V. I.; Zaanen, J.; Andersen, O. K. Band Theory and Mott Insulators: Hubbard *U* Instead of Stoner *I*. Phys. Rev. B **1991**, 44, 943.

(30) Himmetoglu, B.; Floris, A.; de Gironcoli, S.; Cococcioni, M. Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems. Int. J. Quantum Chem. **2014**, 114, 14-49.

(31) Burgess, A. C.; Linscott, E.; O'Regan, D. D. DFT + *U*-type functional derived to explicitly address the flat plane condition. Phys. Rev. B **2023**, 107, L121115.

(32) Kulik, H. J.; Cococcioni, M.; Scherlis, D. A.; Marzari, N. Density Functional Theory in Transition-Metal Chemistry: A Self-Consistent Hubbard *U* Approach. Phys. Rev. Lett. **2006**, 97, 103001.

(33) Scherlis, D. A.; Cococcioni, M.; Sit, P.; Marzari, N. Simulation of Heme Using DFT+U: A Step toward Accurate Spin-State Energetics. J. Phys. Chem. B **2007**, 111, 7384-7391, PMID: 17547444.

(34) Kulik, H. J.; Marzari, N. Systematic Study of First-Row Transition-Metal Diatomic Molecules: A Self-Consistent DFT+U Approach. J. Chem. Phys. **2010**, 133, 114103.

(35) Vela, S.; Fumanal, M.; Cirera, J.; Ribas-Arino, J. Thermal Spin Crossover in Fe(II) and Fe(III). Accurate Spin State Energetics at the Solid State. Phys. Chem. Chem. Phys. **2020**, 22, 4938-4945.

(36) Albavera-Mata, A.; Trickey, S. B.; Hennig, R. G. Mean Value Ensemble Hubbard-U Correction for Spin-Crossover Molecules. J. Phys. Chem. Letters **2022**, 13, 12049-12054, PMID: 36542415.

(37) Mariano, L. A.; Vlaisavljevich, B.; Poloni, R. Biased Spin-State Energetics of Fe(II) Molecular Complexes within Density-Functional Theory and the Linear-Response Hubbard *U* Correction. J. Chem. Theory Comput. **2020**, 16, 6755-6762, PMID: 33108722.

(38) Mariano, L. A.; Vlaisavljevich, B.; Poloni, R. Improved Spin-State Energy Differences of Fe (II) molecular and crystalline complexes via the Hubbard U-corrected Density. J. Chem. Theory Comput. **2021**, 17, 2807-2816.

(39) Nagai, R.; Akashi, R.; Sugino, O. Machine-learning-based exchange correlation functional with physical asymptotic constraints. Phys. Rev. Res. **2022**, 4, 013106.

(40) Li, L.; Hoyer, S.; Pederson, R.; Sun, R.; Cubuk, E. D.; Riley, P.; Burke, K., et al. Kohn-Sham equations as regularizer: Building prior knowledge into machine-learned physics. Phys. Rev. Lett. **2021**, 126, 036401.

(41) Kirkpatrick, J.; McMorrow, B.; Turban, D. H.; Gaunt, A. L.; Spencer, J. S.; Matthews, A. G.; Obika, A.; Thiry, L.; Fortunato, M.; Pfau, D., et al. Pushing the frontiers of density functionals by solving the fractional electron problem. Science **2021**, 374, 1385-1389.

(42) King, D. S.; Truhlar, D. G.; Gagliardi, L. Machine-learned energy functionals for multiconfigurational wave functions. J. Phys. Chem. Lett. **2021**, 12, 7761-7767.

(43) Dick, S.; Fernandez-Serra, M. Highly accurate and constrained density functional obtained with differentiable programming. Phys. Rev. B **2021**, 104, L161109.

(44) Nagai, R.; Akashi, R.; Sugino, O. Completing density functional theory by machine learning hidden messages from molecules. Npj Comput. Mater. **2020**, 6.

(45) Chen, Y.; Zhang, L.; Wang, H.; E, W. Deepks: A comprehensive data-driven approach toward chemically accurate density functional theory. J. Chem. Theory Comput. **2020**, 17, 170-181.

(46) Dick, S.; Fernandez-Serra, M. Machine learning accurate exchange and correlation functionals of the electronic density. Nat. Commun. **2020**, 11, 3509.

(47) Nagai, R.; Akashi, R.; Sasaki, S.; Tsuneyuki, S. Neural-network Kohn-Sham exchange-correlation potential and its out-of-training transferability. J. Chem. Phys. **2018**, 148, 241737.

(48) Brockherde, F.; Vogt, L.; Li, L.; Tuckerman, M. E.; Burke, K.; Müller, K.-R. Bypassing the Kohn-Sham equations with machine learning. Nat. Commun. **2017**, 8, 872.

(49) Li, L.; Snyder, J. C.; Pelaschier, I. M.; Huang, J.; Niranjan, U.-N.; Duncan, P.; Rupp, M.; Müller, K.-R.; Burke, K. Understanding machine-learned density functionals. Int. J. Quantum Chem. **2016**, 116, 819-833.

(50) Snyder, J. C.; Rupp, M.; Hansen, K.; Müller, K.-R.; Burke, K. Finding density functionals with machine learning. Phys. Rev. Lett. **2012**, 108, 253002.

(51) Kalita, B.; Li, L.; McCarty, R. J.; Burke, K. Learning to approximate density functionals. Acc. Chem. Res. **2021**, 54, 818-826.

(52) Kennedy, J.; Eberhart, R. Particle swarm optimization. Proceedings of ICNN'95-International Conference on Neural Networks. 1995; pp 1942-1948.

(53) Tian, Y.; Ha, D. Modern evolution strategies for creativity: Fitting concrete images and abstract concepts. International Conference on Computational Intelligence in Music, Sound, Art and Design (Part of EvoStar). 2022; pp 275-291.
(54) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. *Phys. Rev.* **1965**, *140*, A1133-A1138.

(55) Furness, J. W.; Kaplan, A. D.; Ning, J.; Perdew, J. P.; Sun, J. Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation. *J. Phys. Chem. Lett.* **2020**, *11*, 8208-8215.

(56) Sun, J.; Ruzsinszky, A.; Perdew, J. P. Strongly constrained and appropriately normed semilocal density functional. *Phys. Rev. Lett.* **2015**, *115*, 036402.

(57) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. Assessment of Gaussian-2 and density functional theories for the computation of enthalpies of formation. *J. Chem. Phys.* **1997**, *106*, 1063-1079.

(58) Furness, J. W.; Sun, J. Enhancing the efficiency of density functionals with an improved iso-orbital indicator. *Phys. Rev. B* **2019**, *99*, 041119.

(59) Yang, Z.-h.; Peng, H.; Sun, J.; Perdew, J. P. More realistic band gaps from meta-generalized gradient approximations: Only in a generalized Kohn-Sham scheme. *Phys. Rev. B* **2016**, *93*, 205205.

(60) Sitkiewicz, S. P.; Zalesýn, R.; Ramos-Cordoba, E.; Luis, J. M.; Matito, E. How reliable are modern density functional approximations to simulate vibrational spectroscopies? *J. Phys. Chem. Lett.* **2022**, *13*, 5963-5968.

(61) Bartók, A. P.; Yates, J. R. Regularized SCAN functional. *J. Chem. Phys.* **2019**, *150*, 161101.

(62) Levy, M.; Perdew, J. P. Hellmann-Feynman, virial, and scaling requisites for the exact universal density functionals. Shape of the correlation potential and diamagnetic susceptibility for atoms. *Phys. Rev. A* **1985**, *32*, 2010.

(63) Oliver, G.; Perdew, J. Spin-density gradient expansion for the kinetic energy. *Phys. Rev. A* **1979**, *20*, 397.

(64) Kulik, H. J.; Marzari, N. Systematic study of first-row transition-metal diatomic molecules: A self-consistent DFT+U approach. *J. Chem. Phys.* **2010**, *133*, 114103.

(65) Kasim, M. F.; Iinko, S. M. Learning the Exchange-Correlation Functional from Nature with Fully Differentiable Density Functional Theory. *Phys. Rev. Lett.* **2021**, *127*, 126403.

(66) Dufour, C.; Schamps, J.; Barrow, R. Electronic states of the CuF molecule. II. Nature of the observed states. *J. Phys. B: Atomic Mol. Phys.* (**1968-1987**) **1982**, *15*, 3819.

(67) Brazier, C. R.; Brown, J. M.; Purnell, M. R. Nuclear hyperfine splittings for copper monofluoride in the a3Σ+ state. *J. Mol. Spectrosc.* **1983**, *99*, 279-287.

(68) Drechsler, G.; Boesl, U.; Bäfmann, C.; Schlag, E. Mass selected anion-zero kinetic energy photoelectron spectroscopy (anion-ZEKE): Ground and low excited states of FeO. *J. Chem. Phys.* **1997**, *107*, 2284-2291.

(69) Miller, A. E. S.; Feigerle, C. S.; Lineberger, W. C. Laser photoelectron spectroscopy of CrH⁻, CoH⁻, and NiH⁻: Periodic trends in the electronic structure of the transition-metal hydrides. *J. Chem. Phys.* **1987**, *87*, 1549-1556.

(70) Ram, R.; Jarman, C.; Bernath, P. Fourier transform emission spectroscopy of the A⁶Σ⁺-X⁶Σ⁺ system of CrH: evidence for a A⁴Σ⁺ lowest excited state. *J. Mol. Spectrosc.* **1993**, *161*, 445-454.

(71) Jiang, W.; DeYonker, N. J.; Wilson, A. K. Multireference Character for 3d Transition-Metal-Containing Molecules. *J. Chem. Theory Comput.* **2012**, *8*, 460-467. PMID: 26596596.

(72) Kalita, B.; Pederson, R.; Chen, J.; Li, L.; Burke, K. How Well Does Kohn-Sham Regularizer Work for Weakly Correlated Systems? *J. Phys. Chem. Lett.* **2022**, *13*, 2540-2547. PMID: 35285630.

(73) Cleghorn, C. W.; Engelbrecht, A. P. Particle swarm stability: a theoretical extension using the non-stagnate distribution assumption. *Swarm Intelligence* **2018**, *12*, 1-22.

(74) (ml-electron project), R. N. PySCF version of physically constrained NN-based XC functional. https://github.com/ml-electron-project/pcNN_mol.

(75) Lehtola, S.; Steigemann, C.; Oliveira, M. J.; Marques, M. A. Recent developments in libxc—A comprehensive library of functionals for density functional theory. *SoftwareX* **2018**, *7*, 1-5.

(76) Sun, Q.; Zhang, X.; Banerjee, S.; Bao, P.; Barbry, M.; Blunt, N. S.; Bogdanov, N. A.; Booth, G. H.; Chen, J.; Cui, Z.-H., et al. Recent developments in the PySCF program package. *J. Chem. Phys.* **2020**, *153*, 024109.

(77) Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dulak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C., et al. The atomic simulation environment—a Python library for working with atoms. *J. Phys. Condens. Matter* **2017**, *29*, 273002.

(78) Miranda, L. J. V. PySwarms: a research toolkit for Particle Swarm Optimization in Python. *J. Open Source Softw.* **2018**, *3*, 433.

(79) Neese, F.; Wennmohs, F.; Becker, U.; Riplinger, C. The ORCA quantum chemistry program package. *J. Chem. Phys.* **2020**, *152*, 224108.

(80) DeepMind, DM21 Interface for PySCF. https://github.com/deepmind/deepmind-research/tree/master/density_functionaltec approximationdm21.

(81) Giannozzi, P. et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. *J. Phys. Condens. Matter*

2017, 29, 465901.

(82) Garrity, K. F.; Bennett, J. W.; Rabe, K. M.; Vanderbilt, D. Pseudopotentials for High-Throughput DFT Calculations. *Comp. Mater. Sci.* **2014**, *81*, 446-452.

(83) Makov, G.; Payne, M. C. Periodic Boundary Conditions In Ab Initio Calculations. *Phys. Rev. B* **1995**, *51*, 4014-4022.

(84) Pierloot, K.; Phung, Q. M.; Domingo, A. Spin State Energetics in First-Row Transition Metal Complexes: Contribution of (3s3p) Correlation and its Description by Second-Order Perturbation Theory. *J. Chem. Theory Comput.* **2017**, *13*, 537-553.
Graphical TOC Entry

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