Machine Learning Framework for CIDER Functionals

CIDER functionals are trained with Gaussian-process regression (GPR). This section describes the form used for local energy-density models and integrated electronic-energy observations.

Gaussian Process Regression

Gaussian process regression is a nonparametric Bayesian regression method. A general treatment is given by Rasmussen and Williams. [1] The construction below follows its application to density functionals in CIDER23X and CIDER24X. [2][3]

Let \(f(\mathbf{x})\) denote the learned function, with \(\mathbf{x}\) as its feature vector. Let \(\mathbf{X}\) be a feature matrix, where each row \(\mathbf{x}_i\) is the feature vector for training point \(i\), and let \(\mathbf{y}\) be a target vector, where each element \(y_i\) is the observed value of the target function for training point \(i\).

The covariance kernel \(k(\mathbf{x}, \mathbf{x}')\) defines the covariance between the predictive function at two points: \(k(\mathbf{x}, \mathbf{x}')=\operatorname{Cov}(f(\mathbf{x}), f(\mathbf{x}'))\). The matrix of covariances between training points is denoted \(\mathbf{K}\) with matrix elements \(K_{ij}=k(\mathbf{x}_i, \mathbf{x}_j)\). Observation standard deviations \(\sigma_i\) specify the assumed noise or prior uncertainty. Their variances form the diagonal matrix \(\boldsymbol{\Sigma}_\mathrm{noise}\). The predictive function for a test point \(\mathbf{x}_*\) is [1]

\[f(\mathbf{x}_*) = \sum_i k(\mathbf{x}_*, \mathbf{x}_i) \alpha_i\]

with weight vector

\[\boldsymbol{\alpha} = \left(\mathbf{K} + \boldsymbol{\Sigma}_\text{noise}\right)^{-1} \mathbf{y}\]

For molecular and solid-state training, the observations are extensive quantities obtained by integrating a local predictive function over real space. This changes the covariance construction because each label couples all quadrature points in one chemical system.

Fitting Total Energy Data

Consider an extensive quantity \(F\) that contributes to the electronic energy of a chemical system. In CiderPress, \(F\) can be an exchange, correlation, or exchange-correlation contribution. The regression learns a dimensionless function that is integrated with a baseline energy density:

\[F = \int \mathrm{d}^3\mathbf{r}\, e^{\mathrm{base}}(\mathbf{r})\,f\left(\mathbf{x}(\mathbf{r})\right)\]

Here \(f\left(\mathbf{x}(\mathbf{r})\right)\) is the dimensionless enhancement factor learned by the regression, and \(e^{\mathrm{base}}\) is the multiplicative baseline energy density: the LDA exchange energy density for an exchange kernel, and the PBE correlation energy density for a correlation kernel. For a chemical system indexed by \(m\), numerical quadrature gives

\[F^m = \sum_{g\in m} w_g^m e_g^{\mathrm{base}} f\left(\mathbf{x}_g^m\right)\]

where \(g\) indexes quadrature points and \(w_g^m\) are the quadrature weights. The covariances between the numerical integrals \(F^m\) and \(F^n\) can be written as

\[\text{Cov}(F^m, F^n) = \sum_{g \in m} \sum_{h \in n} w_g^m w_h^n e_g^{\mathrm{base}} e_h^{\mathrm{base}} k(\mathbf{x}_g^m, \mathbf{x}_h^n)\]

where \(k(\mathbf{x}, \mathbf{x}')\) is the covariance kernel for \(f(\mathbf{x})\). CIDER avoids the direct double sum by introducing a set of control points \(\tilde{\mathbf{x}}_a\) and applying a Nyström/resolution-of-the-identity approximation:

\[\begin{split}\operatorname{Cov}(F^m, F^n) = K_{mn} &\approx \tilde{\mathbf{k}}_m^{\mathsf T} \tilde{\mathbf{K}}^{-1}\tilde{\mathbf{k}}_n, \\ \left(\tilde{\mathbf{K}}\right)_{ab} &= k(\tilde{\mathbf{x}}_a, \tilde{\mathbf{x}}_b), \\ \left(\tilde{\mathbf{k}}_m\right)_a &= \sum_{g\in m} w_g^m e_g^{\mathrm{base}} k(\mathbf{x}_g^m, \tilde{\mathbf{x}}_a).\end{split}\]

Let \(\mathbf{K}_\mathrm{sys}\) be the resulting covariance matrix between the integrated system observations. The predictive mean follows from

\[\begin{split}\mathbf{K}_\mathrm{sys} &= [K_{mn}], \\ \boldsymbol{\beta} &= \left(\mathbf{K}_\mathrm{sys} +\boldsymbol{\Sigma}_\mathrm{noise}\right)^{-1}\mathbf{y}, \\ \boldsymbol{\alpha} &= \tilde{\mathbf{K}}^{-1} \sum_m \tilde{\mathbf{k}}_m\,\beta_m, \\ f(\mathbf{x}_*) &= \sum_a k(\mathbf{x}_*,\tilde{\mathbf{x}}_a)\alpha_a.\end{split}\]

The vector \(\mathbf{y}\) contains the integrated training labels. This construction is implemented by ciderpress.models.train.MOLGP.fit(): Kmm is \(\tilde{\mathbf{K}}\), the columns of Kmn are \(\tilde{\mathbf{k}}_m\), alpha_new is \(\boldsymbol{\beta}\), and kernel.alpha contains \(\boldsymbol{\alpha}\).

If several independent GP components contribute to one observation, their system covariance matrices are added before the shared \(\boldsymbol{\beta}\) solve. Each component then obtains its own control-point coefficients. CIDER26XC uses this construction for the learned exchange and correlation terms. [4]

Fitting Eigenvalues

The Gaussian process scheme discussed above can be extended to fit orbital-occupation derivatives. [3] At a stationary generalized Kohn–Sham solution, the derivative of the total energy with respect to the occupation \(f_i^m\) gives the corresponding eigenvalue \(\epsilon_i^m\): [5]

\[\epsilon_i^m = \frac{\partial E}{\partial f_i^m}\]

The exact particle-number derivatives at an integer electron number \(N\) are \(-I\) from below and \(A\) from above, where \(I\) and \(A\) are the ionization potential and electron affinity. In pure Kohn–Sham DFT, the difference between the neutral system’s lowest unoccupied orbital energy and \(A\) contains the XC derivative discontinuity.

Occupation-derivative labels can provide both valence- and conduction-edge information. CIDER24X fits derivatives of the target energy contribution, \(\partial F/\partial f_i^m\).

For exact exchange \(E_\mathrm{x}^\mathrm{exact}\), the occupation derivative can be evaluated explicitly for a given set of orbitals. A full exchange-correlation reference derivative can be constructed from a consistent ionization-potential, electron-affinity, or band-edge target. The total energy in generalized Kohn–Sham DFT is a functional of the one-particle density matrix \(n_1\) (and therefore of its orbitals and occupations):

\[\begin{split}E[n_1] &= T[n_1] + V[n] + U[n] + E_\mathrm{xc}[n_1] \\ &= E_0[n_1] + E_\mathrm{xc}[n_1].\end{split}\]

Here \(E_0[n_1]\) is the sum of the kinetic (\(T\)), external (\(V\)), and Hartree (\(U\)) energies, all of which can be written explicitly in terms of the generalized Kohn–Sham orbitals. Given a reference total occupation derivative, the XC part follows by subtracting the explicit non-XC derivative:

\[\frac{\partial E_\mathrm{xc}[n_1]}{\partial f_i^m} = \epsilon_i^m - \frac{\partial E_0[n_1]}{\partial f_i^m}.\]

The density matrix, orbitals, reference derivative, and side of the integer particle number must all be defined consistently when constructing such a training label.

The occupation derivative must include the response of both the feature vector and the multiplicative baseline:

\[\frac{\partial F^m}{\partial f_i^m} = \sum_{g\in m}w_g^m\frac{\partial}{\partial f_i^m} \left[e_g^{\mathrm{base}}f(\mathbf{x}_g^m)\right].\]

Including this label requires its covariance with an integrated energy \(F^n\) or another occupation derivative. These covariances are

\[\begin{split}\operatorname{Cov}\left(\frac{\partial F^m}{\partial f_i^m}, F^n\right) &= \tilde{\mathbf{d}}_{mi}^{\mathsf T} \tilde{\mathbf{K}}^{-1} \tilde{\mathbf{k}}_n \\ \operatorname{Cov}\left(\frac{\partial F^m}{\partial f_i^m}, \frac{\partial F^n}{\partial f_j^n}\right) &= \tilde{\mathbf{d}}_{mi}^{\mathsf T} \tilde{\mathbf{K}}^{-1} \tilde{\mathbf{d}}_{nj} \\ \left(\tilde{\mathbf{d}}_{mi}\right)_a &= \sum_{g\in m}w_g^m\frac{\partial}{\partial f_i^m} \left[e_g^{\mathrm{base}} k(\mathbf{x}_g^m,\tilde{\mathbf{x}}_a)\right]\end{split}\]

These covariance blocks place integrated energies and occupation derivatives in the same Gaussian-process fit. The complete occupation-derivative derivation is given in Bystrom et al.[3].