Ionic forces and stress tensor in all-electron DFT calculations using enriched finite element basis

Consider a periodic unit cell $\Omega$ containing $N_{e}$ electrons and $N_{a}$ nuclei with ionic position vectors $\boldsymbol{\textbf{R}}=\{\boldsymbol{\textbf{R}}_{1},\boldsymbol{\textbf{R}}_{2},\cdots,\boldsymbol{\textbf{R}}_{N_{a}}\}$. Neglecting spin, the free energy of the system in Kohn-Sham density functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965) at finite temperature Mermin (1965) is given by

$$\begin{split}\mathcal{F}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\mathbf{R})=T_{\mathrm{s}}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\mathbf{R})+&E_{\mathrm{xc}}(\rho)+E_{\mathrm{el}}(\rho,\mathbf{R})-\\ &E_{\mathrm{ent}}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}})\>,\end{split}$$

(1)

where $\mathbf{u}_{\boldsymbol{\textbf{k}}}=\{u_{1,\boldsymbol{\textbf{k}}}(\boldsymbol{\textbf{x}}),u_{2,\boldsymbol{\textbf{k}}}(\boldsymbol{\textbf{x}}),\cdots,u_{N,\boldsymbol{\textbf{k}}}(\boldsymbol{\textbf{x}})\}\>\>(N>N_{e}/2)$ represent electronic wavefunctions corresponding to k, a point in the reciprocal space. These wavefunctions are considered to be non-orthogonal, in general. $\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}$ represents the matrix corresponding to the single-particle density operator ($\hat{\Gamma}$) expressed in the non-orthogonal basis $\mathbf{u}_{\boldsymbol{\textbf{k}}}$. Hence, the elements of the matrix are given by $\Gamma_{pq}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}=\sum_{r=1}^{N}{{S^{\boldsymbol{\textbf{k}}}}^{-1}_{pr}}\langle u_{r,\boldsymbol{\textbf{k}}}|\hat{\Gamma}|u_{q,\boldsymbol{\textbf{k}}}\rangle$, where $S^{\boldsymbol{\textbf{k}}}$ is the overlap matrix whose elements are given by

$$S^{\boldsymbol{\textbf{k}}}_{m,n}=\int_{\Omega}u_{m,\boldsymbol{\textbf{k}}}^{*}(\boldsymbol{\textbf{x}})\,u_{n,\boldsymbol{\textbf{k}}}(\boldsymbol{\textbf{x}})d\boldsymbol{\textbf{x}}\>.$$

(2)

The superscript ‘^∗’ in the above expression denotes complex conjugate. The electron density $\rho(\boldsymbol{\textbf{x}})$ in Eq. (1) can be written in terms of the density matrix and non-orthogonal wavefunctions as follows:

$$\rho(\mathbf{x})=2\fint_{\Omega_{\text{BZ}}}\left\{\sum_{p,q,r=1}^{N}\Gamma_{pq}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\,\,{{S^{\boldsymbol{\textbf{k}}}}^{-1}_{qr}}u_{r,\boldsymbol{\textbf{k}}}^{*}(\boldsymbol{\textbf{x}})\,u_{p,\boldsymbol{\textbf{k}}}(\boldsymbol{\textbf{x}})\right\}d\boldsymbol{\textbf{k}}\>,$$

(3)

where the average integral over the Brillouin zone, in practice, is computed as a discrete sum over k-points lying in the Brillouin zone ($BZ$) with an associated weight $w_{\boldsymbol{\textbf{k}}}$. A common choice for the $k$-point grid is the Monkhorst-Pack (MP) scheme Monkhorst and Pack (1976), which is adopted in this work. In the present work, we restrict our analysis to spin-independent systems. However, all the ideas discussed subsequently can be generalized to spin-dependent systems.

The first constituent of the free energy functional in Eq. (1) is the kinetic energy of the non-interacting electrons, $T_{s}$, which is given by

$$T_{\text{s}}\left(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\mathbf{u}_{\boldsymbol{\textbf{k}}}\right)=2\fint_{\Omega_{\text{BZ}}}\left\{\sum_{p,q,r=1}^{N}\int_{\Omega}\Gamma_{pq}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}{S^{\boldsymbol{\textbf{k}}}}^{-1}_{qr}u_{r,\boldsymbol{\textbf{k}}}^{*}(\mathbf{x})\left(-\frac{1}{2}(\nabla+i\boldsymbol{\textbf{k}})^{2}\right)u_{p,\boldsymbol{\textbf{k}}}(\mathbf{x})d\mathbf{x}\right\}d\boldsymbol{\textbf{k}}\>.$$

(4)

The exchange-correlation functional in the local density approximation (LDA)Ceperley and Alder (1980); Perdew and Zunger (1981) adopted in this work is given by

$$E_{\text{xc}}(\rho)=\int_{\Omega}F(\rho)d\boldsymbol{\textbf{x}}=\int_{\Omega}\varepsilon_{xc}(\rho)\rho(\boldsymbol{\textbf{x}})d\boldsymbol{\textbf{x}}\,.$$

(5)

$E_{\text{el}}$ represents the classical electrostatic interaction energy between the electrons and nuclei and can be written as follows:

$$E_{\text{el}}(\rho,\mathbf{R})=E_{\text{H}}(\rho)+E_{\text{ext}}(\rho,\mathbf{R})+E_{\text{ZZ}}(\boldsymbol{\textbf{R}})\>,$$

(6)

where $E_{\text{H}}$ and $E_{\text{ZZ}}$ are the Hartree energy and the nuclear repulsive energy, respectively, and are given by

$$E_{\mathrm{H}}(\rho)=\frac{1}{2}\int_{\Omega}\int_{\mathbb{R}^{3}}\frac{\rho(\boldsymbol{\textbf{x}})\rho(\boldsymbol{\textbf{x}}^{\prime})}{|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{x}}^{\prime}|}d\boldsymbol{\textbf{x}}^{\prime}d\boldsymbol{\textbf{x}}\>,$$

(7)

and

$$E_{\mathrm{zz}}(\boldsymbol{\textbf{R}})=\frac{1}{2}\sum_{I}\sum_{J,\,J\neq I}\frac{Z_{I}Z_{J}}{\left|\mathbf{R}_{I}-\mathbf{R}_{J}\right|}\>.$$

(8)

Note that, in the above equation, the sum over $J$ represents all nuclei in $\mathbb{R}^{3}$ while the sum over $I$ includes only those in the domain $\Omega$. $E_{\text{ext}}$ in Eq. (6) represents the interaction between electrons and nuclei and is given by

$$E_{\mathrm{ext}}(\rho,\boldsymbol{\textbf{R}})=-\sum_{J}\int_{\Omega}\rho(\mathbf{x})\frac{Z_{J}}{\left|\mathbf{x}-\mathbf{R}_{J}\right|}d\mathbf{x}\>.$$

(9)

Lastly, the entropic energy of the electrons in Eq. (1) is given by

$$E_{\text{ent}}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}})=-2\sigma\fint_{\Omega_{\text{BZ}}}\text{tr}\left[\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\ln\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}+\left(\boldsymbol{I}-\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\right)\ln\left(\boldsymbol{I}-\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\right)\right]d\boldsymbol{\textbf{k}}\>,$$

(10)

where $\sigma=k_{B}T$ with $k_{B}$ denoting the Boltzmann constant and $T$ denoting the electronic temperature.

The electronic ground-state, for given nuclear positions, is given by the following variational problem:

$$\min_{\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\in\mathbb{R}^{N\times N}}\min_{\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\in(H_{\text{per}}^{1}(\Omega))^{N}}\mathcal{F}_{c}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\mathbf{R})\>,\quad\forall\,\boldsymbol{\textbf{k}}\in BZ$$

(11)

where $\mathcal{F}_{c}$ is a constrained free-energy functional given by

$$\mathcal{F}_{c}=\mathcal{F}-\mu\left[2\operatorname{tr}\left(\fint_{\Omega_{\text{BZ}}}\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}d\boldsymbol{\textbf{k}}\right)-N_{e}\right]\>.$$

(12)

In the above equation, $\mu$, which is the Lagrange multiplier corresponding to the constraint on the number of electrons, also denotes the Fermi level. We note that $H_{\text{per}}^{1}(\Omega)$ denotes the Hilbert space of periodic functions such that the functions and their first-order derivatives are square-integrable in $\Omega$.

Here we present the local real-space formulation that forms the basis for the derivation of the configurational forces expression. To begin with, the electrostatic interaction energy Eq. (6) presented in the last section are extended in real-space. However, these extended interactions can be recast into a local variational formulation using auxiliary electrostatic potentials (cf. Das et al. (2015)) as:

$$\begin{split}E_{\mathrm{H}}(\rho)+E_{\mathrm{ext}}(\rho,\boldsymbol{\textbf{R}})+E_{\mathrm{zz}}(\boldsymbol{\textbf{R}})=\max_{\phi\in H_{\text{per}}^{1}(\Omega)}\left\{\int_{\Omega}\Big{[}(\rho(\boldsymbol{\textbf{x}})+b(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}}))\phi(\boldsymbol{\textbf{x}})-\frac{1}{8\pi}\left|\nabla\phi(\boldsymbol{\textbf{x}})\right|^{2}\Big{]}d\boldsymbol{\textbf{x}}\right\}-\\ \sum_{I=1}^{N_{a}}\max_{V_{I}\in H^{1}\left(\mathbb{R}^{3}\right)}\left\{\int\Big{[}b_{I}(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}})V_{I}(\boldsymbol{\textbf{x}})-\frac{1}{8\pi}\left|\nabla V_{I}(\boldsymbol{\textbf{x}})\right|^{2}\Big{]}d\boldsymbol{\textbf{x}}\right\}\>.\end{split}$$

(13)

In the above equation, $b(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}})=\sum_{I}b_{I}(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}})=\sum_{I}-Z_{I}\tilde{\delta}\left(\left|\boldsymbol{\textbf{x}}-\mathbf{R}_{I}\right|\right)$ is the nuclear charge density represented using regularized Dirac-delta functions centered at the nucleus. Also, the second term in Eq. (13) serves the purpose of removing the nuclear self-interaction contribution contained in the first term.

In the present work, however, we adopt the smeared charge representation for the nuclear charge Pask et al. (2012) which is computationally more efficient in a discrete setting. To this end, the nuclear charge density is written as

$$b_{s}(\boldsymbol{\textbf{x}})=\sum_{I}b_{s,I}(\boldsymbol{\textbf{x}})=\sum_{I}-Z_{I}g(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{c,I})\,,$$

(14)

where $g(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{c,I})$ denotes a smeared charge which is localized within $|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|<r_{c,I}$ and integrates to unity. We employ the following form for the unit smeared charge Pask et al. (2012)

$$g(r,r_{c})=\begin{cases}\frac{-21(r-r_{c})^{3}(6r^{2}+3rr_{c}+r_{c}^{2})}{5\pi r_{c}^{8}},&0\leq r\leq r_{c},\\ 0,&r>r_{c}\end{cases}$$

(15)

The $r_{c,I}$’s are chosen to be the largest possible values that avoid overlap between two neighboring smeared charges. We note that, while not explicitly indicated, the smeared charge sphere is considered to wrap around the periodic boundary should it breach the boundaries of the computational domain. The local reformulation for the electrostatic energy using smeared nuclear charge is given by

$$\begin{split}E_{\text{el}}(\rho,\mathbf{R})&=E_{\mathrm{H}}(\rho)+E_{\mathrm{ext}}(\rho,\boldsymbol{\textbf{R}})+E_{\mathrm{zz}}(\boldsymbol{\textbf{R}})\\ &=\max_{\phi\in H_{\text{per}}^{1}(\Omega)}\mathcal{L}_{\text{el}}(\phi,\rho,\mathbf{R})\>,\end{split}$$

(16)

where

$$\begin{split}\mathcal{L}_{\text{el}}(\phi,\rho,\mathbf{R})=&\int_{\Omega}\Big{[}(\rho(\boldsymbol{\textbf{x}})+b_{s}(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}}))\phi(\boldsymbol{\textbf{x}})-\frac{1}{8\pi}\left|\nabla\phi(\boldsymbol{\textbf{x}})\right|^{2}\Big{]}d\boldsymbol{\textbf{x}}+\sum_{I}\int_{\Omega_{I}}\rho(\boldsymbol{\textbf{x}})\Big{[}V_{I}\left(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|\right)-V_{s,I}\left(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{c,I}\right)\Big{]}d\boldsymbol{\textbf{x}}\\ &-\sum_{I}\frac{1}{2}\,\int_{\Omega_{I}}b_{s,I}(\boldsymbol{\textbf{x}})V_{s,I}\left(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{c,I}\right)d\boldsymbol{\textbf{x}}\>.\end{split}$$

(17)

In the above, $(\rho(\boldsymbol{\textbf{x}})+b_{s}(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}}))$ constitutes a charge-neutral system. The second term is a correction term containing the difference between the exact nuclear potential ($V_{I}$) corresponding to the point nuclear charge and the smeared nuclear potential ($V_{s,I}$) corresponding to $b_{s,I}$. These are given by

$$V_{I}(r)=-\frac{Z_{I}}{r}\,,$$

(18)

$$V_{s,I}(r,r_{c,I})=-Z_{I}v_{g}(r,r_{c,I})\,,$$

(19)

where $v_{g}(r,r_{c})$ is the potential corresponding to the $g(r,r_{c})$ and is given by

$$v_{g}(r,r_{c})=\begin{cases}\frac{9r^{7}-30r^{6}r_{c}+28r^{5}r_{c}^{2}-14r^{2}r_{c}^{5}+12r_{c}^{7}}{5r_{c}^{8}},&0\leq r\leq r_{c}\\ \frac{1}{r},&r>r_{c}\,.\end{cases}$$

(20)

We note that since $V_{I}$ and $V_{s,I}$ are identical for $r>r_{c,I}$, the correction for each nucleus is evaluated only inside the sphere of radius $r_{c,I}$ around $\boldsymbol{\textbf{R}}_{I}$, whose domain is denoted by $\Omega_{I}$. The last term in Eq. (17) removes the self-interaction energy corresponding to the smeared nuclear charges.

Finally, the electronic ground-state energy for a given positions of the nuclei R can be computed using Eq. (11) and Eq. (16) as the following variational problem:

$$E_{0}(\mathbf{R})=\mathcal{L}(\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},{\bar{\mathbf{u}}_{\boldsymbol{\textbf{k}}}},\bar{\phi},\mathbf{R})=\min_{\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\in\mathbb{R}^{N\times N}}\min_{\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\in(H_{\text{per}}^{1}(\Omega))^{N}}\max_{\phi\in H_{\text{per}}^{1}(\Omega)}\mathcal{L}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\phi,\mathbf{R})\>,$$

(21)

where

$$\mathcal{L}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\phi,\mathbf{R})=T_{\mathrm{s}}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\boldsymbol{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\mathbf{R})+E_{\mathrm{xc}}(\rho)+\mathcal{L}_{\text{el}}(\phi,\rho,\mathbf{R})-E_{\mathrm{ent}}(\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}})-\mu\left[2\operatorname{tr}\left(\fint_{\Omega_{\text{BZ}}}\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}d\boldsymbol{\textbf{k}}\right)-N_{e}\right]\>.$$

(22)

In the above, $\mathcal{L}$ denotes the local reformulation of the constrained free-energy functional $\mathcal{F}_{c}$, and $\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\bar{\mathbf{u}}_{\boldsymbol{\textbf{k}}}$ and $\bar{\phi}$ denote the extremizers of the variational problem.

The enriched finite element (EFE) basis used in this work is constructed by augmenting the finite element (FE) basis with atom-centered enrichment functions. The enrichment functions, which attempt to capture the highly oscillatory behavior of the wavefunctions around the nuclei, obviate the need for highly refined finite element mesh near the nuclei, thus greatly reducing the computational cost. In this work, we derive configurational forces in the context of enriched finite element basis presented by Kanungo and Gavini Kanungo and Gavini (2017). We note that using the EFE basis tends to yield ill-conditioned matrices in both the generalized eigenvalue problem and the linear system of equations (electrostatics). This affects the robustness of the self-consistent field iterations. There are many approaches Schweitzer (2011); Albrecht et al. (2018) proposed in recent years to address this problem. In particular, Rufus et. al. Rufus et al. (2021) have proposed the orthogonalized enriched FE (OEFE) basis, in which, the atomic solutions are orthogonalized with respect to the underlying FE basis before augmenting as enrichments. The OEFE basis yields well-conditioned matrices ensuring numerical robustness. However, we note that the EFE basis presented in Kanungo and Gavini (2017) and the OEFE basis presented in  Rufus et al. (2021) span the same function space. Thus, in this work we employ the following strategy to compute ionic forces and cell stresses. First, we compute the electronic ground-state using the OEFE basis. We subsequently transform the fields from the OEFE to the EFE basis. Finally, we evaluate the forces and stresses using the expressions presented in this work. We derive the configurational force expressions in the EFE basis to avoid accounting for the additional contributions arising from the orthogonalizing terms in the OEFE basis.

Next, we provide a brief overview of the EFE basis. In the EFE basis, the wavefunctions $u_{\alpha,\boldsymbol{\textbf{k}}}(\boldsymbol{\textbf{x}})$ are approximated as

$$\begin{split}u_{\alpha,\mathbf{k}}(\boldsymbol{\textbf{x}})\approx\,&u_{\alpha,\mathbf{k}}^{h}(\boldsymbol{\textbf{x}})\\ =&\underbrace{\sum_{i=1}^{n_{h}}N^{C}_{i}(\boldsymbol{\textbf{x}})u_{\alpha,\boldsymbol{\textbf{k}},i}^{C}}_{\text{Classical FE}}+\underbrace{\sum_{I=1}^{N_{a}}\sum_{j=1}^{n_{I}}N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})u_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}}_{\text{Enrichment}}\,.\end{split}$$

(23)

In the above equation, the superscript $h$ indicates a discrete field, and the superscript $C$ and $E$ are used to distinguish the classical and the enrichment components, respectively. In particular, $N^{C}_{i}(\boldsymbol{\textbf{x}})$ denotes the $i\textsuperscript{th}$ classical FE basis function, and $u_{\alpha,\boldsymbol{\textbf{k}},i}^{C}$ denotes the expansion coefficient of $N^{C}_{i}(\boldsymbol{\textbf{x}})$ for $u_{\alpha,\mathbf{k}}$. Similarly, $N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}})$ denotes the $k$-point dependent enrichment function for $u_{\alpha,\mathbf{k}}~{}(\forall\alpha)$. The index $I$ runs over all the atoms ($N_{a}$) in the system, and the index $j$ runs over all the atomic Kohn-Sham orbitals ($n_{I}$) we include for the atom $I$. In other words, the $I$^th atom, situated at $\boldsymbol{\textbf{R}}_{I}$, contributes $n_{I}$ enrichment functions, each centered around $\boldsymbol{\textbf{R}}_{I}$. $u_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}$ represents the expansion coefficient of $N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})$ corresponding to $u_{\alpha,\mathbf{k}}$. The reader is referred to Kanungo and Gavini (2017); Rufus et al. (2021) for the form of $N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})$.

In a similar vein, the electrostatic potential $\phi$ is given as

$$\phi(\boldsymbol{\textbf{x}})\approx\phi^{h}(\boldsymbol{\textbf{x}})=\underbrace{\sum_{j=1}^{n_{h}}N^{C}_{j}(\boldsymbol{\textbf{x}})\phi_{j}^{C}}_{\text{Classical FE}}+\underbrace{\sum_{I=1}^{N_{a}}N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\phi_{I}^{E}}_{\text{Enrichment}}\,.$$

(24)

As with the discretization of $u_{\alpha,\mathbf{k}}$ (Eq.( 23)), $N^{C}_{j}(\boldsymbol{\textbf{x}})$ denotes the $j\textsuperscript{th}$ classical FE basis function and $\phi_{j}^{C}$ denotes its corresponding coefficient. Similarly, $N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})$ is the $I\textsuperscript{th}$ enrichment function with a corresponding coefficient $\phi_{I}^{C}$.

Thus, in the finite-dimensional EFE basis, the electronic ground-state energy is given by

$$\begin{split}E_{0}^{h}(\boldsymbol{\textbf{R}})&=\mathcal{L}\bigg{(}\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\bigg{\{}u^{h}_{\alpha,\boldsymbol{\textbf{k}}}\Big{(}\{N^{C}_{i}(\boldsymbol{\textbf{x}})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}\}\Big{)}\bigg{\}},\phi^{h}\Big{(}\{N^{C}_{j}(\boldsymbol{\textbf{x}})\},\{\bar{\phi}_{j}^{C}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{\bar{\phi}_{I}^{E}\}\Big{)},\mathbf{R}\bigg{)}\\ &=\min_{\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}\in\mathbb{R}^{N\times N}}\min_{\{u_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\}\in\mathbb{C}^{N\times{n_{h}}}}\min_{\{u_{\alpha,\boldsymbol{\textbf{k}},\nu}^{E}\}\in\mathbb{C}^{N\times{n_{E}^{u_{\boldsymbol{\textbf{k}}}}}}}\max_{\{\phi_{j}^{C}\}\in\mathbb{R}^{n_{h}}}\max_{\{\phi_{j}^{E}\}\in\mathbb{R}^{N_{a}}}\\ &\,\,\mathcal{L}\bigg{(}\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\bigg{\{}u^{h}_{\alpha,\boldsymbol{\textbf{k}}}\Big{(}\{N^{C}_{i}(\boldsymbol{\textbf{x}})\},\{u_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{u_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}\}\Big{)}\bigg{\}},\phi^{h}\Big{(}\{N^{C}_{j}(\boldsymbol{\textbf{x}})\},\{\phi_{j}^{C}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{\phi_{I}^{E}\}\Big{)},\mathbf{R}\>\bigg{)}\>.\end{split}$$

(25)

In the above equation, $\{\,.\,\}$ are used to denote lists of functions or scalar coefficients running over all indices. $\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}\},\{\bar{\phi}_{j}^{C}\}\text{ and }\{\bar{\phi}_{I}^{E}\}$ denote the discrete ground-state solution (stationarity point of $\mathcal{L}$) for given nuclear positions R and basis functions. These quantities can be computed by solving the Kohn-Sham equations in the EFE basis. In particular, in the context of canonical eigenfunctions, $\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\}$ and $\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},\nu}^{E}\}$ are the solution of the generalized Kohn-Sham eigenvalue problem; $\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}}$ contains the fractional occupancies of the eigenstates given by the Fermi-Dirac distribution; $\{\bar{\phi}_{j}^{C}\}\text{ and }\{\bar{\phi}_{I}^{E}\}$ are solutions to the Poisson problem governing the electrostatic interactions.

We now derive expressions for the configurational force corresponding the electronic ground-state free energy in Eq. (21). We employ the approach of inner variations, where the generalized force corresponding to a perturbation of underlying space is evaluated. To this end, we define a bijective mapping $\boldsymbol{\tau}^{\varepsilon}\,:\mathbb{R}^{3}\mapsto\mathbb{R}^{3}$ which represents an infinitesimal perturbation of the underlying space, mapping a material point x to a new point $\boldsymbol{\textbf{x}}^{\mathbf{\varepsilon}}=\boldsymbol{\tau}^{\varepsilon}(\boldsymbol{\textbf{x}})$. Using a Taylor series expansion in $\varepsilon$, this mapping can be written as

$$\boldsymbol{\tau}^{\varepsilon}(\boldsymbol{\textbf{x}})=\boldsymbol{\textbf{x}}+\varepsilon\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})+\mathcal{O}(\varepsilon^{2})\>,$$

(26)

where $\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})=\left.\frac{d}{d\varepsilon}\boldsymbol{\tau}^{\varepsilon}(\boldsymbol{\textbf{x}})\right|_{\varepsilon=0}$ is the generator of the spatial perturbation. We note that $\boldsymbol{\textbf{x}}^{0}=\boldsymbol{\textbf{x}}$ denotes the reference configuration.

The configurational force corresponding to a prescribed $\boldsymbol{\tau}^{\varepsilon}(\boldsymbol{\textbf{x}})$ is given by

$$\left.\frac{d}{d\varepsilon}\mathcal{L}^{\varepsilon}(\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}},\varepsilon},\bar{\mathbf{u}}^{\varepsilon}_{\boldsymbol{\textbf{k}}},\bar{\phi}^{\varepsilon},\mathbf{R}^{\varepsilon})\right|_{\varepsilon=0}\,,$$

(27)

where $\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}},\varepsilon}$, $\bar{\mathbf{u}}^{\varepsilon}_{\boldsymbol{\textbf{k}}}$ and $\bar{\phi}^{\varepsilon}$ are the extremizers of constrained free energy functional in the perturbed configuration ($\mathcal{L}^{\varepsilon}$). In other words, the configurational force is the Gâteaux derivative of $\mathcal{L}^{\varepsilon}$ in the direction of $\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})$ with all the electronic fields set to the electronic ground-state.

In the discrete setting with the enriched finite element basis, the configurational force is given by

$$\begin{split}\hat{F}^{h}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))=\frac{d}{d\varepsilon}\mathcal{L}^{\varepsilon}\bigg{(}\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}},\varepsilon},\bigg{\{}u_{\alpha,\boldsymbol{\textbf{k}}}^{h}\Big{(}\{N^{C,\varepsilon}_{i}(\boldsymbol{\textbf{x}}^{\varepsilon})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C,\varepsilon}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}_{I}^{\varepsilon})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E,\varepsilon}\}\Big{)}\bigg{\}},\\ \phi^{h}\Big{(}\{N^{C,\varepsilon}_{j}(\boldsymbol{\textbf{x}}^{\varepsilon})\},\{\bar{\phi}_{j}^{C,\varepsilon}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}_{I}^{\varepsilon})\},\{\bar{\phi}_{I}^{E,\varepsilon}\}\Big{)},\mathbf{R^{\varepsilon}}\bigg{)}\bigg{|}_{\varepsilon=0}\,.\end{split}$$

(28)

We note that in computing the configurational force in the discrete setting, we choose the generator from the subspace spanned by the FE basis (cf. Sec II.5) as it allows us to take advantage of the isoparametric nature of FE basis functions. In the above, $\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}},\varepsilon},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C,\varepsilon}\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E,\varepsilon}\},\{\bar{\phi}_{j}^{C,\varepsilon}\}$ and $\{\bar{\phi}_{I}^{E,\varepsilon}\}$ represent the ground-state electronic fields in the discrete setting, i.e., the extremizers of $\mathcal{L}^{\varepsilon}$ in subspace spanned by the enriched FE basis. It is important to note here that the enrichment functions $\{N_{j,I}^{E,u_{\boldsymbol{\textbf{k}}}}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}^{\varepsilon})\}$ and $\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}^{\varepsilon})\}$ retain their form even when the underlying space is deformed, as these are a priori computed enrichment functions independent of the finite element discretization.

For the sake of further simplification, we express the constrained free energy functional in terms of the energy density $f$ as

$$\begin{split}&\mathcal{L}\bigg{(}\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\bigg{\{}u_{\alpha,\boldsymbol{\textbf{k}}}^{h}\Big{(}\{N^{C}_{i}(\boldsymbol{\textbf{x}})\},\{u_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{u_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}\}\Big{)}\bigg{\}},\phi^{h}\Big{(}\{N^{C}_{j}(\boldsymbol{\textbf{x}})\},\{\phi_{j}^{C}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{\phi_{I}^{E}\}\Big{)},\mathbf{R}\bigg{)}\\ &=\int_{\Omega}f\bigg{(}\boldsymbol{\Gamma}^{\mathbf{u}_{\boldsymbol{\textbf{k}}}},\bigg{\{}u_{\alpha,\boldsymbol{\textbf{k}}}^{h}\Big{(}\{N^{C}_{i}(\boldsymbol{\textbf{x}})\},\{u_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{u_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}\}\Big{)}\bigg{\}},\phi^{h}\Big{(}\{N^{C}_{j}(\boldsymbol{\textbf{x}})\},\{\phi_{j}^{C}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{\phi_{I}^{E}\}\Big{)},\mathbf{R}\bigg{)}d\boldsymbol{\textbf{x}}\>.\end{split}$$

(29)

The configurational force is decomposed as follows (refer to Appendix for details):

$$\begin{split}\hat{F}^{h}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))&=\frac{d}{d\varepsilon}\Bigg{[}\int_{\Omega^{\varepsilon}}f^{\varepsilon}\bigg{(}\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}},\varepsilon},\bigg{\{}u_{\alpha,\boldsymbol{\textbf{k}}}^{h}\Big{(}\{N^{C,\varepsilon}_{i}(\boldsymbol{\textbf{x}}^{\varepsilon})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C,\varepsilon}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}_{I}^{\varepsilon})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E,\varepsilon}\}\Big{)}\bigg{\}},\\ &\qquad\qquad\qquad\phi^{h}\Big{(}\{N^{C,\varepsilon}_{j}(\boldsymbol{\textbf{x}}^{\varepsilon})\},\{\bar{\phi}_{j}^{C,\varepsilon}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}_{I}^{\varepsilon})\},\{\bar{\phi}_{I}^{E,\varepsilon}\}\Big{)},\mathbf{R^{\varepsilon}}\bigg{)}d\boldsymbol{\textbf{x}}^{\varepsilon}\Bigg{]}\bigg{|}_{\varepsilon=0}\\ &=\hat{F}_{1}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))+\hat{F}_{2}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))\>,\end{split}$$

(30)

where $\hat{F}_{1}$ denotes the contribution to the configurational force that arise from all dependencies on $\varepsilon$ excepting those that arise from the enrichment functions, and $\hat{F}_{2}$ denotes the contribution to the configurational force that arise from the enrichment functions. In particular,

$$\begin{split}\hat{F}_{1}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))&=\frac{d}{d\varepsilon}\Bigg{[}\int_{\Omega^{\varepsilon}}f^{\varepsilon}\bigg{(}\bar{\boldsymbol{\Gamma}}^{\mathbf{u}_{\boldsymbol{\textbf{k}}},\varepsilon},\bigg{\{}u_{\alpha,\boldsymbol{\textbf{k}}}^{h}\Big{(}\{N^{C}_{i}(\boldsymbol{\textbf{x}})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C,\varepsilon}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E,\varepsilon}\}\Big{)}\bigg{\}},\\ &\qquad\qquad\qquad\phi^{h}\Big{(}\{N^{C}_{j}(\boldsymbol{\textbf{x}})\},\{\bar{\phi}_{j}^{C,\varepsilon}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\},\{\bar{\phi}_{I}^{E,\varepsilon}\}\Big{)},\mathbf{R^{\varepsilon}}\bigg{)}d\boldsymbol{\textbf{x}}^{\varepsilon}\Bigg{]}\bigg{|}_{\varepsilon=0}\,.\end{split}$$

(31)

We note that, in writing the above expression, we have made use of the isoparametric nature of the FE basis functions, i.e., $N^{C,\varepsilon}_{i}(\boldsymbol{\textbf{x}}^{\varepsilon})=N^{C}_{i}(\boldsymbol{\textbf{x}})$. We note that the configurational force is of a similar form as the one considered in Motamarri & Gavini Motamarri and Gavini (2018), albeit in a discrete setting. Thus, using the previously derived results Motamarri and Gavini (2018), we note that we can evaluate $\hat{F}_{1}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))$ in terms of the canonical Kohn-Sham eigenfunctions as follows:

$$\begin{split}\hat{F}_{1}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))=\int_{\Omega}\mathbf{E}:\nabla\boldsymbol{\Upsilon}(\mathbf{x})d\mathbf{x}+\mathrm{F}^{K}+\mathrm{F}^{\text{S}}\>,\end{split}$$

(32)

where the Eshelby tensor $\mathbf{E}$ is given by

$$\begin{split}\mathbf{E}=\Bigg{(}\fint_{\Omega_{\text{BZ}}}\Bigg{[}\sum_{\alpha}\bar{f}_{\alpha,\boldsymbol{\textbf{k}}}\>\bigg{\{}\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\cdot\nabla\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})\,-2i\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\boldsymbol{\textbf{k}}\cdot\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})+(|\boldsymbol{\textbf{k}}|^{2}-2\bar{\epsilon}_{\alpha,\mathbf{k}})|\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})|^{2}\bigg{\}}\Bigg{]}d\boldsymbol{\textbf{k}}\>\\ +\varepsilon_{xc}(\bar{\rho}^{h})\bar{\rho}^{h}(\boldsymbol{\textbf{x}})+(\bar{\rho}^{h}(\boldsymbol{\textbf{x}})+b_{s}(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}}))\bar{\phi}^{h}(\boldsymbol{\textbf{x}})-\frac{1}{8\pi}\big{|}\nabla\bar{\phi}^{h}(\boldsymbol{\textbf{x}})\big{|}^{2}\Bigg{)}\boldsymbol{I}\\ -\fint_{\Omega_{\text{BZ}}}\Bigg{[}\sum_{\alpha}\bar{f}_{\alpha,\boldsymbol{\textbf{k}}}\>\bigg{\{}\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\otimes\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})+\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})\otimes\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})-2i\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\big{(}\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})\otimes\boldsymbol{\textbf{k}}\big{)}\bigg{\}}\Bigg{]}d\boldsymbol{\textbf{k}}\\ +\frac{1}{4\pi}\nabla\bar{\phi}^{h}(\boldsymbol{\textbf{x}})\otimes\nabla\bar{\phi}^{h}(\boldsymbol{\textbf{x}})\>,\end{split}$$

(33)

and

$$\mathrm{F}^{K}=\fint_{\Omega_{\text{BZ}}}\Bigg{[}\sum_{\alpha}\bar{f}_{\alpha,\boldsymbol{\textbf{k}}}\int_{\Omega}\>\bigg{\{}-2i\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\left.\frac{d}{d\varepsilon}\boldsymbol{\textbf{k}}^{\varepsilon}\right|_{\varepsilon=0}\cdot\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})+\left.\frac{d}{d\varepsilon}|\boldsymbol{\textbf{k}}^{\varepsilon}|^{2}\right|_{\varepsilon=0}|\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})|^{2}\bigg{\}}d\boldsymbol{\textbf{x}}\Bigg{]}d\boldsymbol{\textbf{k}}\>.$$

(34)

We note that, in the above expression, $\boldsymbol{\textbf{k}}^{\varepsilon}$ depends on the nature of the generator. For a Gaussian-type generator used in the force calculations (cf. Sec. II.5), since there is no change in the lattice vectors, we have

$$\boldsymbol{\textbf{k}}^{\varepsilon}(\boldsymbol{\textbf{k}})=\boldsymbol{\textbf{k}}\implies\left.\frac{d}{d\varepsilon}\boldsymbol{\textbf{k}}^{\varepsilon}\right|_{\varepsilon=0}=\mathbf{0}$$

(35)

However, under affine deformations, this dependence is given by:

$$\boldsymbol{\textbf{k}}^{\varepsilon}(\boldsymbol{\textbf{k}})=(\mathbf{I}-\varepsilon\mathbf{C}^{T})\boldsymbol{\textbf{k}}\implies\left.\frac{d}{d\varepsilon}\boldsymbol{\textbf{k}}^{\varepsilon}\right|_{\varepsilon=0}=-\mathbf{C}^{T}\boldsymbol{\textbf{k}}\>,$$

(36)

where $\mathbf{C}$ is a second order tensor independent of x.
$\mathrm{F}^{\text{S}}$ in Eq.(33) contains contributions from the smeared nuclear charges, and is given by

$$\begin{split}\mathrm{F}^{\text{S}}=\sum_{I}\Bigg{[}\int_{\Omega_{I}}\bar{\rho}^{h}(\boldsymbol{\textbf{x}})\nabla\Big{(}V_{I}\left(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|\right)-V_{s,I}\left(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{c,I}\right)\Big{)}\cdot\Big{(}\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})-\boldsymbol{\Upsilon}(\boldsymbol{\textbf{R}}_{I})\Big{)}d\boldsymbol{\textbf{x}}\Bigg{]}\\ +\sum_{I}\Bigg{[}\int_{\Omega_{I}}\bar{\rho}^{h}(\boldsymbol{\textbf{x}})\Big{\{}V_{I}\left(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|\right)-V_{s,I}\left(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{c,I}\right)\Big{\}}\nabla\cdot\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})d\boldsymbol{\textbf{x}}\Bigg{]}\\ -\sum_{I}\Bigg{[}\int_{\Omega_{I}}\bar{\phi}^{h}(\boldsymbol{\textbf{x}})Z_{I}\nabla g(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{c,I})\cdot\Big{(}\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})-\boldsymbol{\Upsilon}(\boldsymbol{\textbf{R}}_{I})\Big{)}d\boldsymbol{\textbf{x}}\Bigg{]}\>.\end{split}$$

(37)

We now turn our attention to $\hat{F}_{2}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))$ in Eq.(30), which contains the contributions arising from the $\varepsilon$-dependence of the enrichment functions. We provide here the expression for $\hat{F}_{2}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))$, and refer to the Appendix for the details of the derivation:

$$\begin{split}\hat{F}_{2}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))=&\fint_{\Omega_{\text{BZ}}}\Bigg{[}\sum_{\alpha}\bar{f}_{\alpha,\boldsymbol{\textbf{k}}}\int_{\Omega}\>\bigg{\{}\nabla\,\Big{(}\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon\>*}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\Big{)}\cdot\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})+\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\cdot\nabla\,\Big{(}\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\Big{)}\\ &-2i\,\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon\>*}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\boldsymbol{\textbf{k}}\cdot\nabla\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})-2i\,\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\boldsymbol{\textbf{k}}\cdot\nabla\,\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}+|\boldsymbol{\textbf{k}}|^{2}\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon\>*}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})\\ &\qquad\qquad+|\boldsymbol{\textbf{k}}|^{2}\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\bigg{\}}d\boldsymbol{\textbf{x}}\Bigg{]}d\boldsymbol{\textbf{k}}\\ &+\int_{\Omega}\left.\frac{d}{d\varepsilon}\widetilde{\rho^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\big{[}V_{\text{xc}}(\boldsymbol{\textbf{x}})+\bar{\phi}^{h}(\boldsymbol{\textbf{x}})+\sum_{I}(V_{I}(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|)-V_{s,I}(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|))\big{]}d\boldsymbol{\textbf{x}}\\ &\qquad\qquad+\int_{\Omega}(\bar{\rho}^{h}(\boldsymbol{\textbf{x}})+b_{s}(\boldsymbol{\textbf{x}},\boldsymbol{\textbf{R}}))\frac{d}{d\varepsilon}\widetilde{\phi^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\bigg{|}_{\varepsilon=0}d\boldsymbol{\textbf{x}}-\int_{\Omega}\frac{1}{4\pi}\nabla\big{(}\frac{d}{d\varepsilon}\widetilde{\phi^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\big{)}\bigg{|}_{\varepsilon=0}\cdot\nabla\bar{\phi}^{h}(\boldsymbol{\textbf{x}})d\boldsymbol{\textbf{x}}\\ &-\fint_{\Omega_{\text{BZ}}}\Bigg{[}\sum_{\alpha}\bar{f}_{\alpha,\boldsymbol{\textbf{k}}}\int_{\Omega}\>2\bar{\epsilon}_{\alpha,\mathbf{k}}\bigg{\{}\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon\>*}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h}(\boldsymbol{\textbf{x}})+\bar{u}_{\alpha,\boldsymbol{\textbf{k}}}^{h,*}(\boldsymbol{\textbf{x}})\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}\bigg{\}}d\boldsymbol{\textbf{x}}\Bigg{]}d\boldsymbol{\textbf{k}}\>,\end{split}$$

(38)

where

$$\begin{split}\left.\frac{d}{d\varepsilon}\widetilde{u_{\alpha,\mathbf{k}}^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}&\coloneqq\left.\frac{d}{d\varepsilon}u_{\alpha,\boldsymbol{\textbf{k}}}^{h}\Big{(}\{N^{C}_{i}(\boldsymbol{\textbf{x}})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},i}^{C}\},\{N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}_{I}^{\varepsilon})\},\{\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}\}\Big{)}\right|_{\varepsilon=0}\\ &=\sum_{I=1}^{N_{a}}\sum_{j=1}^{n_{I}}\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}\>\nabla N^{E,u_{\boldsymbol{\textbf{k}}}}_{j,I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\cdot\Big{(}\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})-\boldsymbol{\Upsilon}(\boldsymbol{\textbf{R}}_{I})\Big{)}\>\end{split}$$

(39)

and

$$\begin{split}\left.\frac{d}{d\varepsilon}\widetilde{\phi^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}&\coloneqq\frac{d}{d\varepsilon}\phi^{h}\Big{(}\{N^{C}_{j}(\boldsymbol{\textbf{x}})\},\{\bar{\phi}_{j}^{C}\},\{N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}^{\varepsilon}-\boldsymbol{\textbf{R}}_{I}^{\varepsilon})\},\{\bar{\phi}_{I}^{E}\}\Big{)}\bigg{|}_{\varepsilon=0}\\ &=\sum_{I=1}^{N_{a}}\bar{\phi}_{I}^{E}\>\nabla N^{E,\phi}_{I}(\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I})\cdot\Big{(}\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}})-\boldsymbol{\Upsilon}(\boldsymbol{\textbf{R}}_{I})\Big{)}\end{split}$$

(40)

account for contributions from the enrichment functions. We note that $\left.\frac{d}{d\varepsilon}\widetilde{\rho^{h,\varepsilon}}(\boldsymbol{\textbf{x}}^{\varepsilon})\right|_{\varepsilon=0}$ can be computed using Eq. (39). We remark that $\hat{F}_{2}(\boldsymbol{\Upsilon}(\boldsymbol{\textbf{x}}))$ vanishes when the enrichment functions are absent, i.e., when $\bar{u}_{\alpha,\boldsymbol{\textbf{k}},j,I}^{E}=\bar{\phi}_{I}^{E}=0$.