1 Introduction↩︎

Many fields have benefited from Kohn-Sham density functional theory (KSDFT) for electronic structure calculations [1]. Thousands of papers using KSDFT are published every year. Despite the progress made in computing science since the method first appeared, ab initio calculations for large systems or at high temperatures remain elusive, costly, and sometimes out of reach. This limitation primarily stems from the orthogonalization of the wave functions (also called orbitals¹) associated with the energy levels of electrons [3]. As a rule of thumb, the calculation scales as \(N^3\), where \(N\) is the number of orbitals. The more atoms there are to simulate, the more electrons are present, and the more orbitals must be computed. At finite temperatures, the situation is exacerbated because electrons occupy higher energy levels, resulting in significantly more orbitals to compute [4].

This scaling issue prevents extensive applications of KSDFT. Indeed, many physical situations of interest require simulating a very large number of atoms. This need is driven by the high structural complexity of molecular chemistry and materials physics (see [5] and references therein). Phenomena outside local thermodynamic equilibrium conditions, such as mass diffusion at an interface or non-Newtonian fluids, also demand large simulation cells. Conversely, while only a few fields require an accurate account of temperature effects on electronic structure, they are crucial for advancing our knowledge. Astrophysics is one of them, encompassing a vast array of conditions ranging from giant planets to stellar interiors and supernova explosions [6]. Temperature-dependent electronic structure is also critical in another, more technological field that has recently seen impressive breakthroughs: inertial confinement fusion (ICF) [7]–[9].

These fields share a common need for hydrorad simulations. Consequently, there is a constant demand for new data to describe material behavior, including equations of state (EOS), opacities, and transport coefficients. In the context of ICF, these quantities must be provided on short timescales to keep pace with rapid technological advances. They are required over an extremely wide range of conditions, spanning densities from \(10^{-5}\) to \(10^5\) g/cm\(^3\) and temperatures from \(300\) K to \(10^9\) K, including the particularly challenging regime known as warm dense matter (WDM) [10].

Accurate data at these very high temperatures can be obtained thanks to recent and ongoing efforts to extend KSDFT applications to higher temperatures [11]–[17]. Unfortunately, these advances require significant computing resources,
which are incompatible with the rapid turnaround required in the ICF technology cycle between experiments and simulations. In the present work, we use one of these KSDFT extensions as a reference. This extension bypasses the \(N^3\) scaling bottleneck by replacing the highest orbitals with plane waves excluded from the orthogonalization computation [12]. This approach is implemented as an orbital-based option in the Abinit code called Extended DFT (originally known as Extended FPMD) [13], [18]. Other alternatives have also been proposed and are currently being tested [11], [14]–[17].

To meet the material property demands of ICF, a fair trade-off must be struck between accuracy and efficiency. Another way to circumvent the KSDFT scaling issue is to use orbital-free density functional theory (OFDFT) coupled with Born-Oppenheimer molecular dynamics (OFMD) [19], [20]. The OFDFT approach relies primarily on a local density approximation for the kinetic energy [21]. With this approximation, orbitals are no longer needed. This speeds up calculations while maintaining accuracy, at least at high temperatures. For low temperatures, a more refined approximation of the non-interacting kinetic energy functional is required to restore accuracy. Semi-local corrections involving the density gradient exist [22], as well as non-local functionals utilizing a kernel that connects the density at two different points [23], [24]. These investigations into kinetic energy functionals parallel those for exchange and correlation functionals [25]. With OFDFT, since only the electron density is computed instead of many orbitals, the computation time scales as \(N\) rather than \(N^3\) (see Fig. 1). As noted in [26], the advantages of a fast code expand the domain of DFT applications (e.g., reaching higher temperatures and improving statistics). It opens up new possibilities, such as preprocessings of KSDFT simulations to reach equilibrium configurations, rapid generation of low-fidelity data for machine learning [27], and the calculation of dynamic properties (like viscosity) that require many atoms and timesteps [28].

Codes implementing OFDFT can be categorized by their treatment of temperature effects. For zero-\(T\) conditions and large systems, profess [29] and DFTpy [30] implement periodic boundary conditions, but real-space implementations also exist to handle more general boundary conditions [31]–[37]. Unfortunately, these codes cannot be utilized for ICF applications, where temperature effects are paramount. At finite temperatures, there are fewer OFMD codes available: profess@qe [38], ofmd [39], and more recently dragon [26].

In this paper, we present Archê, a finite-temperature OFMD code. The core design of the code is based on the approach described in [40], with several notable differences: the use of more recent programming languages and paradigms;
GPU support; the inclusion of unit and integration tests; a scalable three-dimensional Fourier transform; the adoption of a self-consistent field (SCF) method instead of free energy minimization via a conjugate-gradient algorithm (which, to our knowledge, is employed by all existing OFDFT codes); initialization of the electronic density using the converged density from the previous MD timestep and a density-mixing strategy following [41], which both significantly accelerate SCF convergence. Some of these developments arise from general physical considerations and can be readily transferred to other modeling frameworks.

Archê is a parallel code written in C++,
CUDA, and Python. It is linked to HeFFTe [42], [43], a scalable fast Fourier transform (FFT) library. This code serves a dual purpose: first, to quickly produce EOS for WDM; second, to explore numerical methods and implementations that enhance performance, as well as to test other non-interacting free energy functionals to improve accuracy [22], [44]. A particular effort has been made to build a flexible and modular framework so that ideas can be tested without compromising performance. Following a free-function–based programming approach [45], all energy contributions and functionals are implemented as independent functions that can be freely combined or replaced. To substantiate the presentation of the code, we detail its validation process and assess its performance on both CPU and GPU architectures.

The article is organized as follows. Section 2 summarizes the OFMD theoretical framework. Section 3 describes its implementation in Archê. Section 4 validates the code by comparing its results with more accurate models. Finally, Section 5 highlights the performance of Archê through various benchmarks.

2 Theoretical background↩︎

In this section, we present the model as implemented in the Archê code. Unless otherwise stated, equations are written in atomic units (\(e = 4 \pi \varepsilon_0 = \hbar = m_e = 1\)).

2.1 Molecular dynamics, the main loop↩︎

The purpose of a molecular dynamics simulation is to evolve a system of \(N_i\) atoms over time. For each atom, this is achieved by solving Newton’s equation for the force \(\vec{F}_j\) acting on the atomic nucleus \(j\): \[\vec{F}_j = m_j \vec{a}_j,\]

where \(m_j\) is the mass of the nucleus and \(\vec{a}_j\) is its acceleration.

The specificity of each MD model lies in the accuracy or scale at which the forces are computed. The common aspect in MD models is time discretization. There are several schemes, each with its own properties and an associated statistical ensemble in which the simulation is performed [46], [47]. We chose the Verlet integration algorithm, which enables simulations in the microcanonical ensemble (the isokinetic ensemble algorithm, also implemented in Archê, follows similar lines [48]).

The Verlet scheme consists of two steps. The first step determines the new position of the nucleus \(\vec{r}_j\) at time \(t + \Delta t\), where \(\Delta t\) is the MD timestep of the simulation:

\[\vec{r}_j(t + \Delta t) = 2 \vec{r}_j(t) - \vec{r}_j(t - \Delta t) + \frac{\vec{F}_j(t)}{m_j} \Delta t^2\, .\]

The second step computes the velocities \(\vec{v}_j\) at time \(t\) from the positions at times \(t + \Delta t\) and \(t - \Delta t\): \[\vec{v}_j(t) = \frac{\vec{r}_j(t + \Delta t) + \vec{r}_j(t - \Delta t)}{2 \Delta t}\, .\]

These steps are sufficient to evolve the system, provided that the forces are known. Since the propagating particles are charged nuclei, the forces \(\vec{F}_j\) that set the \(N_i\) nuclei in motion are of Coulombic origin. These forces can be broken down into two contributions, \(\vec{F}_j^{ii}\) and \(\vec{F}_j^{ie}\), which are derived from the corresponding potentials: \(U^{ii}\) and \(U^{ie}\) for the nucleus-nucleus and nucleus-electron interactions, respectively.

At each time \(t\), each \(\vec{F}_j(t)\) depends on the positions of the entire system \(R(t) = (\vec{r}_1(t), \ldots, \vec{r}_{N_i}(t))\) and is given by \[\begin{align} \label{eq:force} &\vec{F}_j(t) &= \vec{F}_j^{ii}(t) + \vec{F}_j^{ie}(t) \\ & &= -\nabla_{\vec{r}_j}(U^{ii}(R(t)) + U^{ie}(R(t))).\nonumber \end{align}\tag{1}\] In the following, we assume the Born-Oppenheimer approximation, in which the electronic response is instantaneous, ensuring that the electronic density is in equilibrium with the nuclear positions at any given time [49]. The electron-nucleus interaction \(U^{ie}\) uses the electronic density, computed within the density functional theory framework summarized below. The nucleus-nucleus interaction \(U^{ii}\) is evaluated using Ewald summation (see Sec.2.6).

2.2 Orbital-free model↩︎

The electronic density \(n\) is a continuous field over the domain \(\Omega\) that contains a system of \(N_e\) electrons, defined as \[\label{defN} N_e = \int_\Omega n(\vec{r}) d\vec{r}.\tag{2}\]

Hohenberg and Kohn proved in 1964 [50] that for an external potential \(V_{\text{ext}}(\vec{r})\), there exists a functional \(\mathcal{E}[n]\) that yields the total ground-state energy \(E_0\) via: \[\label{minE} E_0 = \underset{n}{\text{min}}\left(\mathcal{E}[n] + E_{\text{ext}}[n]\right).\tag{3}\]

where \(E_{\text{ext}}[n] = \int_\Omega n(\vec{r}) V_{\text{ext}}(\vec{r}) d\vec{r}\) represents the energy of the electrons in the potential \(V_{\text{ext}}\).

When solving Eq.(3 ), KSDFT uses orbitals (which are more precise but computationally expensive), whereas OFDFT models work solely with \(n\) (see the comprehensive review in [19] for a detailed comparison). In both cases, \(n\) must integrate to \(N_e\), and the problem is equivalent to extremizing the Lagrangian subject to the constraint: \[\begin{align} &\mathcal{L}[n(\vec{r})] = &\mathcal{E}[n(\vec{r})] + E_{\text{ext}}[n(\vec{r})] \\ & &- \mu \left(\int n(\vec{r}) d\vec{r} - N_e\right),\nonumber \end{align}\] where \(\mu\) represents the chemical potential, i.e., the Fermi energy at \(T=0\).

At finite temperatures, the energy functional \(\mathcal{E}\) is replaced by a temperature-dependent functional \(\mathcal{F}\), and the free energy \(F\) is obtained by extremizing the following Lagrangian [51]: \[\begin{align} \label{Lagrangian} & \mathcal{L}[n(\vec{r})] = &\mathcal{F}[n(\vec{r}), T] + E_{\text{ext}}[n(\vec{r})] \\ & &- \mu \left(\int n(\vec{r}) d\vec{r} - N_e\right).\nonumber \end{align}\tag{4}\]

\(\mathcal{F}\) can be broken down into three different contributions: a non-interacting free-energy functional \(\mathcal{F}_0\), the Hartree energy functional \(E_H\), and the exchange and correlation free-energy functional \(\mathcal{F}_{\text{xc}}\) [20], [52]. Minimizing Eq. 4 with respect to \(n\) leads to the following Euler equation: \[\label{Eulerian} \frac{\delta \mathcal{F}_0}{\delta n} + V_{\text{eff}} = \mu\tag{5}\] with \(V_{\text{eff}}:= V_H + V_{\text{xc}} + V_{\text{ext}}\), where \(V_{\text{xc}}:= \frac{\delta \mathcal{F}_{\text{xc}}}{\delta n}\) and \(V_H := \frac{\delta E_H}{\delta n} = \int_\Omega \frac{n(\vec{r}\,')}{|\vec{r} - \vec{r}\,'|} d\vec{r}\,'.\)

A century ago, Thomas and Fermi independently found an analytic solution to Eq.(5 ) for an ideal gas of non-interacting electrons in a mean-field effective potential \(V_{\text{eff}}\) (see Ref. [53] for the derivation): \[\label{SCF1} n(\vec{r}, \mu) = \frac{(2 k_B T_e)^{3/2}}{2 \pi^2} I_{1/2} \left( \frac{\mu - V_{\text{eff}} (\vec{r}) }{k_B T_e}\right) .\tag{6}\] Here, \(I_{1/2}(.)\) is the Fermi-Dirac integral of order one-half defined as \(\operatorname {I} _{j}(x) = \int _{0}^{\infty }\!{\frac{t^{j}}{\mathrm {e} ^{t-x}+1}}\;\mathrm {d} t\), \(k_B\) is the Boltzmann constant, and \(T_e\) is the electronic temperature. For EOS purposes, \(V_{\text{ext}} := V^{\text{ie}}\), where \(V^{\text{ie}}\) is the Coulomb potential of the nuclei, and \(V_{\text{H}}\), denoted as \(V^{\text{ee}}\), is the Hartree potential solution to the Poisson equation: \[\label{SCF2} \Delta V^{\text{ee}}(\vec{r}) = -4 \pi n(\vec{r}, \mu).\tag{7}\]

Within the Thomas-Fermi (TF) model, two main sources of error can be identified. The first arises from the assumption of a uniform electron gas for the functional \(\mathcal{F}_0\), while the second is associated with the mean-field approximation. Several refinements of this model exist, including the von Weizsäcker gradient correction to the uniform electron gas [54] and exchange-correlation effects within the mean-field approximation [55].

2.3 Discretization of the periodic cell↩︎

Exploiting the three-dimensional periodicity of the system, the Poisson equation is solved in Fourier space using FFTs, which requires a regular grid representation in both real and reciprocal spaces.

We study a periodic system composed of an infinite repetition of the unit cell \(\Omega\). In the real spatial domain, it is a box of dimensions \([L_x \times L_y \times L_z]\), and \(\widehat{\Omega}\) is the reciprocal cell in Fourier space. \(\Omega_n\) (resp. \(\widehat{\Omega}_n\)) is the discretized \([N_x \times N_y \times N_z]\) grid of \(\Omega\) (resp. \(\widehat{\Omega}\)) on points \(\vec{r}_{\vec{n}}\) (resp. \(\vec{k}_{\vec{n}}\)) defined by \[\vec{r}_{\vec{n}} = n_x \frac{L_x}{N_x} \, \vec{e}_x + n_y \frac{L_y}{N_y} \, \vec{e}_y + n_z \frac{L_z}{N_z} \, \vec{e}_z \,\] with \(n_\mu \in \left[0, N_\mu\right]\) and \[\vec{k}_{\vec{n}} = n_x \frac{1}{L_x} \, \vec{e}_x + n_y \frac{1}{L_y} \, \vec{e}_y + n_z \frac{1}{L_z} \, \vec{e}_z \,\] with \(n_\mu \in \left[-N_\mu / 2 + 1, N_\mu/ 2 \right]\).

The unitary Fourier transform and its inverse, in ordinary frequency, connecting \(\Omega\) and \(\widehat{\Omega}\) are given by: \[\begin{align} \widehat{f}(\vec{k}) &= \int_\Omega f(\vec{r}) \, e^{-2 \pi i \vec{k} \cdot \vec{r}} d \vec{r}\,, \\ f(\vec{r}) &= \int_{\widehat{\Omega}} \widehat{f}(\vec{k}) \, e^{2 \pi i \vec{r} \cdot \vec{k}} d \vec{k}\,. \end{align}\]

In Fourier space, the solution to the Poisson equation Eq.@eq:SCF2 is: \[\begin{align} \widehat{V}^{\text{ee}}(\vec{k}) &= \frac{\widehat{n}(\vec{k}, \mu)}{\pi |\vec{k}|^2}\,,\quad \vec{k} \neq \vec{0} \\ \widehat{V}^{\text{ee}}(\vec{0}) &= 0 \, . \end{align}\]

2.4 Pseudopotential↩︎

The Coulomb potential of the nuclei, \(V^{\text{ie}}\), which enters Eq.@eq:SCF1 , is the sum of the contributions of all \(N_i\) nuclei present in the simulation cell \(\Omega\) and all their periodic images in the replicas of \(\Omega\): \[V^{\text{ie}}(\vec{r}) = \sum_{\vec{R}} \sum_s \sum_{j \in s} V_{s}(\vec{r} - \vec{r}_j + \vec{R})\, ,\] where \(V_{s}(\vec{r})\) is the Coulomb potential of the nuclei of species \(s\) to which nucleus \(j\) belongs, \(V_s(\vec{r}) = \frac{Z_s}{|\vec{r}\,|}\,\) and the outer sum over \(\vec{R}\) covers all lattice translation vectors \(\vec{R} = n_x L_x \vec{e}_x + n_y L_y \vec{e}_y + n_z L_z \vec{e}_z\) for all relative integers \((n_x, n_y, n_z)\).

However, \(\Omega_n\), the discretized grid of \(\Omega\), does not allow for an accurate resolution of the SCF cycle with the Coulomb potential because it is too stiff close to the nucleus. We therefore employ a smoother norm-conserving pseudopotential \(\widetilde{V}_{s}\) [56].

The potential \(V^{\text{ie}}\) is then obtained in the reciprocal cell \(\widehat \Omega_n\) in Fourier space by: \[\label{eq:iepseudo} \widehat V^{\text{ie}}(\vec{k}) = \sum_s \widehat{\widetilde{V}_{s}}(\vec{k}) \, \widehat S_{s}(\vec{k})\, ,\tag{8}\] where \(\widehat{S_s}\) is the structure factor of species \(s\): \[\widehat{S_s}(\vec{k}) = \sum_{j \in s} e^{-2 \pi i \vec{k} \cdot \vec{r_j}} .\]

Eq.@eq:eq:iepseudo is valid only for \(\vec{k} \neq \vec{0}\). Indeed, the pseudopotential behaves like a Coulombic potential at long range, and so \(\widehat{\widetilde{V}_{s}}(\vec{0})\) diverges. The solution is to simply set \(\widehat{\widetilde{V}_{s}}(\vec{0}) = 0\), as noted by Payne [57]: "[\(\ldots\)] the total ionic potential at \(k=0\) is infinite, so the electron-ion energy is infinite. However, there are similar divergences in the Coulomb energies due to the electron-electron interactions and the ion-ion interactions. The Coulomb \(k=0\) contributions to the total energy from the three interactions cancel exactly."

For each density and temperature, we compute the pseudopotential following the approach proposed by Lambert in Ref. [40]. This method is used in OFMD codes such as dragon [26] and ofmd [39]. First, an electronic density \(n(r)\) is obtained within an average-atom model [21]. This density is a non-linear response to the stiff Coulomb potential. The norm-conserving method consists of imposing smoother behavior close to the nucleus while conserving the charge within a radius \(r_{\text{cut}}\) and the regularity at \(r_{\text{cut}}\). The pseudo-density \(\tilde{n}(r)\) is then given by:

\[\tilde{n}(r) = \begin{cases} \exp(a + b r^2 + c r^4), & r < r_{\text{cut}} \\ n(r), & r \geq r_{\text{cut}} \end{cases}\]

The parameters \(a,b,c\) are determined by the following constraints: \[\begin{align} & \tilde{n}(r_{\text{cut}}) = n(r_{\text{cut}}) \\ &\left. \frac{\partial \tilde{n}}{\partial r} \right|_{r_{\text{cut}}} = \left. \frac{\partial n}{\partial r} \right|_{r_{\text{cut}}} \\ &\int_{0}^{r_{\text{cut}}} dr \, 4\pi r^2 \tilde{n}(r) = \int_{0}^{r_{\text{cut}}} dr \, 4\pi r^2 n(r) \end{align}\]

Once the pseudo-density is known, the pseudopotential is obtained by inverting Eq.@eq:SCF1 such that \(\tilde{n}\) minimizes the Lagrangian in Eq.@eq:Lagrangian when the Coulomb potential of the nuclei is replaced by the pseudopotential:

\[\begin{align} & &\widetilde{V}_{s}(r) = \mu - k_{\mathrm{B}} T_{\mathrm{e}} I_{1/2}^{-1} \left( \frac{2 \pi^2 \tilde{n}(r)}{(2 k_{\mathrm{B}} T_{\mathrm{e}})^{3/2}} \right) ~~~~~~\\ & &- \frac{1}{r} \int_{0}^{r} dx \, 4\pi x^2 \tilde{n}(x) - \int_{r}^{r_{\text{cut}}} dx \, 4\pi x \tilde{n}(x) \nonumber \\ & &- V_{\text{xc}}[\tilde{n}],\nonumber \end{align}\]

where \(I_{1/2}^{-1}\) is the inverse Fermi-Dirac function.

As shown in Fig. 2, the pseudo-density is less polarized in the vicinity of the nucleus than the exact density. Both the pseudo-density and the pseudopotential deviate from their exact counterparts inside the cutoff radius \(r_{\text{cut}}\). This cutoff radius is chosen as a fraction of the Wigner–Seitz radius \(r_{\text{ws}}\). This radius defines the sphere enclosing the average atom and is related to the ionic number density \(\rho\) through \[\frac{4}{3}\pi \rho \, r_{\text{ws}}^{3} = 1\,.\]

Consequently, the pseudopotential is consistent with both the density and the temperature of the simulation. Within the average-atom framework, the pressure depends only on the value of the density \(n\) at \(r_{\text{ws}}\); it is therefore unaffected by the use of a norm-conserving pseudopotential. The same behavior is expected in OFMD simulations.

In contrast, the total energy depends on the density throughout the entire Wigner–Seitz sphere and is therefore influenced by the pseudo-density. Because the exact density is more localized and more strongly bound near the nucleus, a systematic energy difference arises between the pseudo and exact descriptions. This correction is evaluated within the average-atom model and subsequently applied to the energies obtained from OFMD simulations employing pseudopotentials.

For instance, in the average-atom framework, the energy values obtained for aluminum at \(\rho_0\) and 100 eV (corresponding to Fig. 2) are \(-6550\) eV for the Coulombic potential and \(-2971\) eV for the pseudopotential.

2.5 Self-consistent method↩︎

Solving Eq.@eq:SCF1 directly is not possible a priori because the chemical potential \(\mu\) is unknown, and the potential \(V^{\text{ee}}\) depends on \(n\) via Eq.@eq:SCF2 . It is feasible, however, using a self-consistent field (SCF) method combined with Newton’s algorithm. This process is described below.

It starts with an input guess: \(n_i\), which is \(n_0 \equiv N_e / | \Omega |\) for the first MD timestep (see also Sec.2.5.1 for subsequent MD timesteps). This guess serves as a source term to solve the Poisson equation, Eq.@eq:SCF2 , for \(V^{\text{ee}}\). Once \(V^{\text{ee}}\) is obtained, \(V_i\) is the first input potential used as the effective potential \(V_{\text{eff}}\) in Eq.@eq:SCF1 . The latter equation computes an active field \(n_a\) with the current value of the chemical potential \(\mu\).

The following equation is then solved to find the value of \(\mu\):

\[\label{newtonmu} \int_\Omega n_a(\vec{r}, \mu) d\vec{r} - N_e = 0.\tag{9}\]

To this end, Newton’s method is applied starting from an initial guess \(\mu_0\). Given the previously computed \(V_i\), it converges within a few iterations (typically fewer than five) to the chemical potential \(\mu_a\) and the corresponding electronic density \(n_a\).

The electronic density \(n_a\) is combined with the input guess \(n_i\) (see Sec.2.5.2) to generate an SCF iteration density output \(n_o\), which is then used as an updated input value \(n_i\) for the next SCF iteration. The process is repeated until the difference \(||n_i - n_a||_{L^2} < \varepsilon\) is sufficiently small. In this method, the density is never negative or null and always complies with the electroneutrality constraint required by Eq. 9 .

2.5.1 Initializing electron density↩︎

It is possible to speed up the convergence of the SCF loop by setting an initial density close to the solution. During the development of Archê, we observed that from one timestep to the next, the electronic densities were only slightly different. Using this information, our first attempt was to use the last computed electronic density as an initial guess, but this did not improve convergence. What does improve convergence is initializing the density as a function of the nuclei positions, as follows.

We assume that the electron density \(n_i(\vec{r})\) is the superposition of as many functions \(Z_s g_s(\vec{r})\) as there are species in the system, centered on each nucleus (which is a realistic approximation in a very dilute gas):

\[n_i(\vec{r}) = \sum_{\vec{R}} \sum_s \sum_{j \in s} Z_s\,g_s(\vec{r} - \vec{r}_j + \vec{R})\, ,\]

where the lattice translation vectors \(\vec{R} = n_x L_x \vec{e}_x + n_y L_y \vec{e}_y + n_z L_z \vec{e}_z\) for all relative integers \((n_x, n_y, n_z)\). The Fourier representation in the reciprocal cell \(\widehat \Omega_n\) is given by: \[\widehat{n_i}(\vec{k}) = \sum_s Z_s\,\widehat{g_s}(\vec{k})\, \widehat{S}_s(\vec{k}) ,\] where \(\widehat{S_s}\) is the structure factor of species \(s\). Now, we assume that all the normalized profiles \(g_s(\vec{r})\) are equal to a single profile \(g(\vec{r})\). This assumption yields: \[\widehat{n_i}(\vec{k}) = \widehat{g}(\vec{k})\, \widehat{S}_Z(\vec{k}) ,\] with the charge-averaged structure factor \(\widehat{S}_Z = \sum_s Z_s\,\widehat{S}_s\). Consequently, \(\widehat{g} = \widehat{n_i}/\widehat{S}_Z\).

Therefore, the algorithm consists of computing and storing \(\widehat{g^{n}}\) at time \(n\) to initialize the density at time \(n+1\) using:

\[\widehat{n_i^{n+1}} = \widehat{S_Z^{n+1}} \widehat{g^{n}}\]

These additional computations are offset by the resulting speedup in convergence (Sec. 2.5.3).

2.5.2 Mixing electron density↩︎

The mixing algorithm implemented in Archê is an adaptation of the algorithm described in the appendix of More et al. [41]. Originally designed for an average-atom model, we applied it to OFMD in Archê, where it significantly improved convergence (see Sec. 2.5.3). To our knowledge, this has never been done for OFMD simulations. Here, we elaborate on the formulas from More et al.’s paper while maintaining the same notation.

The algorithm computes an output density \(n_o\) using the active density \(n_a\) computed during the SCF cycle (see the beginning of Sec. 2.5) and the input density \(n_i\). Setting \(\Delta n = n_i - n_a\), we have \(n_o = n_a + w \Delta n\), and this mixture of the two electronic densities serves as the new input density for the next iteration of the SCF loop. In practice, \(w\) must often be chosen conservatively — for example, \(n_o = 0.95 n_i + 0.05 n_a\) — which usually requires many iterations.

More et al. utilized all available information. The goal of the algorithm is to select, at each step, an input density that is as close to the converged density as possible. This density extremizes the Lagrangian in Eq.(4 ) with respect to variations in the parameter \(w\), using densities of the form \(n = n_a + w \Delta n\). This density conserves the electroneutrality by construction. Therefore, minimization is achieved when \(\frac{d F[n_o(\vec{r})]}{dw} = 0\). It is reasonable to approximate the free energy \(F\) using a functional expansion in the small variation \(\Delta n\). At second order, this allows us to determine a value of \(w\) that minimizes this approximation of \(F\): \[\begin{align} & &F[n(\vec{r})] \simeq F[n_a(\vec{r})] + w \int_\Omega \frac{\delta F}{\delta n(\vec{r})} \Delta n(\vec{r}) d\vec{r}\nonumber \\ & & + \frac{1}{2} w^2 \int_\Omega \frac{\delta^2 F}{\delta n(\vec{r}) \delta n(\vec{r'})} \Delta n(\vec{r}) \Delta n(\vec{r'}) d\vec{r} d\vec{r'}\nonumber \\ & & + \mathcal{O}(\Delta n^3)\, , \end{align}\] where the functional derivatives are evaluated at \(n_a\). This leads to \(w = U / (K + U)\), with: \[U = -\int_\Omega \frac{\delta F}{\delta n(\vec{r})} \Delta n(\vec{r}) d\vec{r},\] \[U + K = \int_\Omega \frac{\delta^2 F}{\delta n(\vec{r}) \delta n(\vec{r'})} \Delta n(\vec{r}) \Delta n(\vec{r'}) d\vec{r} d\vec{r'}.\] The first functional derivative reads: \[\frac{\delta F}{\delta n(\vec{r})} = \frac{\delta \mathcal{F}_0}{\delta n(\vec{r})} + V_{\text{eff}}.\] For simplicity, we omit the \(V_{\text{xc}}\) contribution. Its inclusion is straightforward and does not alter the conclusions. The second functional derivative then reads: \[\frac{\delta^2 F}{\delta n(\vec{r}) \delta n(\vec{r'})} = \frac{\delta^2 \mathcal{F}_0}{\delta n(\vec{r}) \delta n(\vec{r'})} + \frac{1}{{|\vec{r} - \vec{r}\,'|}}.\] The evaluation of \(V_{\text{eff}}\) at \(n_a\) is denoted \(V_a\). However, \(n_a\) is computed with a different potential, \(V_i\), which is treated as an external potential. The Euler equation used for the computation of \(n_a\) is: \[\frac{\delta \mathcal{F}_0}{\delta n(\vec{r})} + V_i = \mu,\] leading to: \[\begin{align} & \dfrac{\delta^2 \mathcal{F}_0}{\delta n(\vec{r}) \delta n(\vec{r'})} &= \dfrac{d\mu}{dn}\, \delta(\vec{r} - \vec{r}\,') \\ & &= \left(\dfrac{dn}{d\mu}\right)^{-1}\, \delta(\vec{r} - \vec{r}\,') .\nonumber \end{align}\] We use the second right-hand side of this equation since \(\dfrac{dn}{d\mu}\) is computed during Newton’s method to find the value of \(\mu\) that ensures electroneutrality. The complete evaluation of the functional derivatives at \(n_a\) gives: \[\begin{align} &\dfrac{\delta F}{\delta n(\vec{r})} &= \mu_a - V_i + V_a \\ & &= \mu_a - \int_\Omega \frac{\Delta n(\vec{r'})}{|\vec{r} - \vec{r'}|}d\vec{r'},\nonumber \end{align}\] and \[\frac{\delta^2 F}{\delta n(\vec{r}) \delta n(\vec{r'})} = \left(\dfrac{dn_a}{d\mu_a}\right)^{-1} \delta(\vec{r} - \vec{r}\,') + \frac{1}{{|\vec{r} - \vec{r}\,'|}}.\]

This derivation led More et al. to select this specific value \(w = U / (K + U)\) for the mixing parameter, where: \[U = \int_\Omega \frac{\Delta n(\vec{r})\Delta n(\vec{r'})}{|\vec{r} - \vec{r'}|}d\vec{r}d\vec{r'},\] \[K = \int_\Omega \frac{[\Delta n(\vec{r})]^2}{\partial n_a(\vec{r}) / \partial \mu_a} d\vec{r}.\]

In practice, \(w\) is computed and updated only after a few SCF iterations, once the chemical potential fluctuations are small. Despite the additional computations, this method pays off by reducing the overall wall-clock time of the simulations.

2.5.3 Convergence speed up↩︎

Fig.3 provides a numerical example of the improvements achieved through the new density mixing and initialization algorithms. This example displays the total number of SCF iterations required for a simulation of 101 molecular dynamics timesteps. Algorithm (A) represents the naive setup, which initializes the electron density at the start of the SCF loop with a constant value and mixes the input and active densities in fixed proportions at each iteration to generate an output density. Setup (B) incorporates only the improved mixing algorithm based on More’s approach
(Sec.2.5.2). Setup (C) includes only the initialization improvement using a superposition approximation with an average one-center density profile (Sec.2.5.1). Finally, version (D) is the complete implementation utilizing both improvements. This example, like many others, demonstrates that when both improvements are combined, the speedup can approach a factor of 6. Furthermore, under many conditions, the mixing algorithm not only accelerates convergence but actively ensures it. This enhancement is highly valuable for both the speed and the robustness of the code — a combination rare enough to be worth emphasizing.

2.6 Forces and thermodynamic quantities↩︎

Once the total energy of the system is known as a function of nuclear positions and velocities, along with the electron density contribution, the forces acting on the nuclei can be computed as the derivative of the energy with respect to the nuclear positions, according to Eq.@eq:eq:force . The pressure is computed as the trace of the stress tensor, whose contributions are obtained either from the free energy via a functional derivative with respect to the deformation tensor [58] or from the forces using a generalized virial theorem [59].

The different contributions to the total energy are the kinetic energies of the nuclei (\(K^\text{i}\)) and electrons (\(K^\text{e}\)), the interaction between nuclei (\(U^{\text{ii}}\)), the interaction between electrons (\(U^{\text{ee}}\)), and the interaction between nuclei and electrons (\(U^{\text{ie}}\)). Therefore, we split the total energy \(E\) as follows: \[E = K^\text{i} + K^\text{e} + U^\text{ii} + U^{\text{ie}} + U^{\text{ee}} + \Delta E_\text{AA}.\] A correction \(\Delta E_\text{AA}\) is added to compensate for the regularization of the pseudopotential (see Sec. 2.4).

The forces involve only the contributions \(U^\text{ii}\) and \(U^\text{ie}\), which depend explicitly on the nuclear positions. As previously stated: \[\begin{align} &\vec{F}_j(t) &= \vec{F}_j^{ii}(t) + \vec{F}_j^{ie}(t) \\ & &= -\nabla_{\vec{r}_j}(U^{ii}(R(t)) + U^{ie}(R(t))).\nonumber \end{align}\]

The pressure can be derived from the trace of the so-called average virial stress tensor (also known as the pressure tensor; see [60], [61]): \[P = - \frac{1}{3} \text{tr}(\boldsymbol{\Sigma}).\] There are as many contributions to the stress tensor as there are to the energy (excluding the correction \(\Delta E_\text{AA}\), which arises from comparing the Coulombic and pseudopotentials within the average atom model): \[\boldsymbol{\Sigma} = \boldsymbol{\sigma^i} + \boldsymbol{\sigma^e} + \boldsymbol{\sigma^{ii}} + \boldsymbol{\sigma^{ee}} + \boldsymbol{\sigma^{ie}}.\]

2.6.1 Kinetic energy of nuclei↩︎

The kinetic contributions follow from the principles of statistical mechanics [58]. The kinetic energy is given by: \[\label{def::K95i} K^i = \frac{1}{2}\sum^{N_i}_i m_i \boldsymbol{v}^2_i\, ,\tag{10}\] and its contribution to the stress tensor reads: \[\label{def::sigma95i} \sigma^i_{\mu \nu} = -\frac{2}{|\Omega|} \sum^{N_i}_i \frac{m_i v_{\mu} v_{\nu}}{2}\, .\tag{11}\]

2.6.2 Kinetic energy of electrons↩︎

The TF kinetic energy is given by [40], [53], [62]: \[K_e = k_B T_e \int \mathrm{d}\vec{r}\, c_{\text{TF}} \, \mathrm{I}_{3/2}[\eta(\vec{r})]\] with \(c_{\text{TF}} = \frac{\sqrt{2}}{\pi^2} (k_B T_e)^{3/2}\) and \[\eta(\vec{r}) = \left[\mu - V_{\text{eff}}(\vec{r}) \right] / (k_B T_e).\]

Since the non-interacting free energy density \(f_0\), defined by \(\mathcal{F}_0 = \int \mathrm{d}\vec{r}\;f_0(n)\), is a local functional of the density \(n\) in the TF model, the stress is diagonal and reads: \[\begin{align} &\sigma^e_{\mu \nu} &= \frac{1}{|\Omega|}\int \mathrm{d}\vec{r}\,\left[f_0[n] - n \frac{\partial f_0}{\partial n} \right] \delta_{\mu \nu} \nonumber\\ & &= -\frac{2 K_e}{3 |\Omega|}\, \delta_{\mu \nu} . \end{align}\] Corrections to the TF model incorporate the density gradient through semi-local functionals of the von Weizsäcker type [19], [20], [22], as well as non-local functionals involving density values evaluated at two different spatial positions [19], [20], [23], [24]. Expressions for the corresponding contributions to the kinetic energy and stress tensor can be found in the supplemental material of Ref.[63].

2.6.3 Interaction energy between nuclei↩︎

The interaction between nuclei requires summing the long-range Coulomb interactions between the particles in the simulation box and all their infinite periodic images. This is evaluated using Ewald summation, and is therefore split into three components: short-range (particle-particle), long-range (particle-mesh), and a neutralizing jellium background. The separation between the long- and short-range components is controlled by the Ewald parameter \(\alpha\) [2], [64]. \[U^\text{ii} = U^\text{pp} + U^\text{pm} + U^\text{je}\]

We selected a value for \(\alpha\) that allows the short-range interaction to be evaluated using the minimum image convention [47].

The short-range energy \(U^\text{pp}\) reads: \[U^\text{pp} = \sum^{N_i}_{j = 1} \sum^{N_i}_{k > j} Z_j Z_k \frac{\operatorname{erfc} (\alpha |\vec{r}_{jk}|)}{|\vec{r}_{jk}|}\,.\] The long-range energy \(U^{pm}\) reads: \[\begin{align} \label{eq:ewald:longrange} &U^\text{pm} &= \frac{1}{2 |\Omega|} \sum_{\vec{k}\,\in\,\widehat{\Omega}_n,\,\vec{k} \neq 0} \widehat{V}_{\alpha}(\vec{k}) | \widehat{S}_Z(\vec{k}) |^2 \nonumber\\ & &- \frac{\alpha}{\sqrt{\pi}} \sum^{N_i}_{j = 1} Z_j^2\,, \end{align}\tag{12}\] with: \[\widehat{V}_{\alpha}(\vec{k}) = \frac{1}{\pi} \frac{e^{- ({\pi k}/{\alpha})^2}}{k^2} .\]

The value \(\vec{k} = 0\) is excluded from the summation because the system is neutralized by a jellium of opposite charge \(Q = \sum_j^{N_i} Z_j\), as explained in Appendix F of Ref.[46]. The second term on the right-hand side of Eq. 12 corrects for the self-interaction energy that was artificially introduced in the first one.

The energy contribution due to the jellium is given by: \[E^\text{je} = - \frac{\pi}{2 \alpha^2 |\Omega|} Q^2 .\]

The jellium does not contribute to the forces, which are split into short-range (pp) and long-range (pm) components as follows: \[\vec{F}_j^{ii}(t) = \vec{F}_j^\text{pp} + \vec{F}_j^\text{pm} ,\]

\[\begin{align} & &\vec{F}^\text{pp}_j = - \sum^{N_i}_{k \neq j} Z_j Z_k \\ & & \times \left( \frac{\operatorname{erfc} (\alpha |\vec{r}_{jk}|)}{|\vec{r}_{jk}|} + \frac{2 \alpha}{\sqrt{\pi}} \operatorname{exp}(-\alpha^2 |\vec{r}_{jk}|^2) \right)\frac{\vec{r}_{jk}}{ |\vec{r}_{jk}|^2} ,\nonumber \end{align}\] and \[\begin{align} &\vec{F}^\text{pm}_j = &Z_j\,\dfrac{2 \pi}{|\Omega|}~ \sum_{\vec{k}\,\in\,\widehat{\Omega}_n,\,\vec{k} \neq 0} \widehat{V}_{\alpha}(\vec{k}) \\ & &\times \,\text{Im}\left( e^{i 2 \pi \vec{k} \cdot \vec{r}_j} \widehat{S}_Z(\vec{k}) \right) \vec{k} .\nonumber \end{align}\]

The jellium contributes to the stress tensor, which is also split into short-range (pp) and long-range (pm) components: \[\boldsymbol{\sigma^\text{ii}} = \boldsymbol{\sigma^\text{pp}} + \boldsymbol{\sigma^\text{pm}} + \boldsymbol{\sigma^\text{je}} ,\] with \[\begin{align} & &\sigma^\text{pp}_{\mu \nu} = -\frac{1}{|\Omega|} \sum^{N_i}_{i = 1} \sum^{N_i}_{j > i} Z_i Z_j \\ & &\times\left( \frac{\operatorname{erfc} (\alpha |\vec{r}_{ij}|)}{|\vec{r}_{ij}|} + \frac{2 \alpha}{\sqrt{\pi}} \operatorname{exp}(-\alpha^2 |\vec{r}_{ij}|^2) \right) \frac{r_{ij}^\mu r_{ij}^\nu}{ |\vec{r}_{ij}|^2} ,\nonumber \end{align}\] \[\begin{align} &\sigma^\text{pm}_{\mu \nu} &= \frac{1}{2 |\Omega|^2} \sum_{\vec{k}\,\in\,\widehat{\Omega}_n,\,\vec{k} \neq 0} \widehat{V}_{\alpha}(\vec{k})\\ & & \times \left| \widehat{S}_Z(\vec{k}) \right|^2\left[ 2 k_\mu k_\nu \left( \frac{\pi^2}{ \alpha^2} + \frac{1}{k^2} \right) - \delta_{\mu \nu}\right] ,\nonumber \end{align}\] and \[\sigma^\text{je}_{\mu \nu} = \frac{\pi}{2 \alpha^2 |\Omega|^2} Q^2\, \delta_{\mu \nu}\;.\]

2.6.4 Interaction energy between electrons↩︎

The electrostatic energy can be computed in several equivalent ways:

\[\begin{align} &U^\text{ee} &= \frac{1}{2 } \int_{\Omega} n (\vec{r}) \, V^\text{ee} (\vec{r}) \, d\vec{r}\\ & &= \frac{1}{2 } \int_{\widehat \Omega} \widehat n (\vec{k}) \, \overline{\widehat V^\text{ee}} (\vec{k}) \, d\vec{k} \,,\nonumber \end{align}\] with their discretized versions: \[\begin{align} &U^\text{ee} &= \frac{|\Omega|}{2N} \sum_{\vec{r}\,\in\,{\Omega}_n}\, n (\vec{r}) \,V^{\text{ee}} (\vec{r}) \\ & &= \frac{1}{2 |\Omega|} \sum_{\vec{k}\,\in\,\widehat{\Omega}_n,\,\vec{k} \neq 0} \widehat{n} (\vec{k}) \,\overline{\widehat V^\text{ee}} (\vec{k}) \,. \nonumber \end{align}\]

The tensor formula is derived from Appendix G of Ref.[46]:

\[\sigma^\text{ee}_{\mu \nu} = \frac{1}{2 |\Omega|^2} \sum_{\vec{k}\,\in\,\widehat{\Omega}_n,\,\vec{k} \neq 0} \widehat{n}(\vec{k}) \overline{ \widehat{V}^\text{ee} }(\vec{k}) \left[ 2 \frac{k_\mu k_\nu }{k^2} - \delta_{\mu \nu}\right]\, .\]

Quantum corrections due to exchange and correlation are not considered here. They can take the form of local density functionals or semi-local functionals that depend on the density gradient [20].

2.6.5 Interaction energy between nuclei and electrons↩︎

The general formula incorporating \(V^{\text{ie}}\), the complete pseudopotential defined in Eq.@eq:eq:iepseudo , is: \[\begin{align} &E^\text{ie} &= \int_\Omega n(\vec{r}) V^{\text{ie}} (\vec{r}) \,d\vec{r} \\ & &-\langle n \rangle \sum_{s} \int_0^{\infty} N_s \left( \widetilde{V}_s(r) - \frac{Z_s}{r} \right) 4 \pi r^2 dr.\nonumber \end{align}\] where \(\langle n \rangle = Q / |\Omega|\) and \(N_s\) is the total number of ions of species \(s\). The second term on the right-hand side accounts for the non-Coulomb part of the pseudopotential at \(k=0\) (see Sec. II.D.3 of Ref.[57]).

For each nucleus \(j\) of species \(s\), the total electron force acting upon it involves the pseudopotential \(\widetilde{V}_s\), according to (see Ref. [46] for the derivation and [29] for an example implementation in the profess code): \[\begin{align} \label{Fie} &\vec{F}^{ie}_{j,s} &= -\nabla_{\vec{r}_j}(U^{ie}(R(t))) \\ & &= 2\pi \sum_{\vec{k}\,\in\,\widehat{\Omega}_n} \widehat{\widetilde{V}}_s(\vec{k})\operatorname{Im}\left(\,\overline{\widehat{n}(\vec{k})} e^{-2\pi i \vec{k}\cdot \vec{r_j}}\right) \vec{k}\,.\nonumber \end{align}\tag{13}\]

Accordingly, the stress is given by [46]: \[\begin{align} & &\sigma^\text{ie}_{\mu \nu} = -\frac{E^\text{ie} }{|\Omega|} \delta_{\mu \nu}\\ & &-\frac{1}{|\Omega|^2} \sum_{\vec{k} \neq 0} \widehat{n}(\vec{k}) \frac{k_\mu k_\nu}{|\vec{k} |} \sum_s \sum_{j \in s} \frac{\partial \widehat{\widetilde{V_s}} (|\vec{k}|)}{\partial|\vec{k}|} \, e^{ 2 \pi i \vec{k} \cdot \vec{r}_j} \nonumber \end{align}\]

3 Archê implementation↩︎

Archê is the successor to the Fortran code ofmd developed by Lambert [40]. The most computationally demanding parts are written in C++ and CUDA, while Python is used for pre- and post-processing. It is parallelized using MPI and relies on the HeFFTe fast Fourier transform library [42]. It can run on both CPU and GPU architectures, achieving an average ten-fold speedup on a single GPU compared to 256 CPU cores. It is interfaced with the libxc library for exchange-correlation functionals [65] at zero temperature and includes its own implementation of the finite-temperature exchange-correlation functional KSDT [55], translated from its Fortran implementation in the Abinit code [66]. Table [CodeNumber] provides key statistics regarding the codebase.

The flowchart in Fig. 4 summarizes the overall architecture of the Archê code. Calculation loops are represented by three rectangles enclosing their respective routines. It illustrates that at each MD timestep, the self-consistent solution is computed within a nested SCF loop, which itself contains another nested loop: the Newton’s algorithm used to determine the chemical potential. SCF iterations account for over 95% of the total computation time on CPUs. Consequently, the number of SCF iterations is the primary factor determining execution time. For this reason, our initial development efforts focused on reducing this number (see Sections 2.5.1 and 2.5.2).

The SCF loop is implemented in two versions. The first is written in standard C++ for CPU architectures, while the second is implemented in C++/CUDA for GPUs. FFT calls are handled through a wrapper that dispatches them to the appropriate FFT library depending on the architecture (CPU or GPU). Both versions are parallelized using MPI. Only the electronic density calculation is parallelized via domain decomposition in both real and reciprocal spaces, primarily to leverage parallel FFTs. Communication costs are the main bottleneck for the scalability of the Fourier transforms, and consequently for the entire computation, on both CPU and GPU architectures. This effect is even more pronounced on GPUs.

Archê was validated through term-by-term comparisons with Lambert’s ofmd code. Development was guided by these comparisons, and updates were accepted only when relative differences were below \(10^{-5}\) for all force components. After successive iterations, excellent agreement between ofmd and Archê was achieved, and unit tests were implemented to prevent any future regressions in the results.

4 Validation↩︎

Archê was validated by comparing its internal energy and pressure results for aluminum with those obtained using the Extended KSDFT
model implemented in Abinit [13]. This reference model, hereafter referred to as Abinit Extended, explicitly accounts for electronic orbitals.

Simulations were run on GPU using a discrete 3D spatial grid of \(256^3\) points and the Thomas-Fermi functional. The timestep \(\Delta t\) was chosen to be less than one-tenth of the smaller of the inverse plasma frequency \(\tau_{\mathrm{p}}\) and the ratio of the system length \(L\) to the thermal velocity \(v_{\mathrm{th}}\): \[\Delta t < 0.1\min\left(\tau_{\mathrm{p}}, L/v_{\mathrm{th}}\right).\]

To extract an equation of state from a simulation, the contributions to the internal energy and pressure are averaged over time. As shown in Fig.5, this must be done after a thermalization phase, which, in our simulations, represents around 20% of the total duration of 3000 timesteps. Several system sizes were tested \(N_i \in \{32, 108,\) \(128, 256\}\), and the simulation yielding the smallest variance was selected. Variances were computed using non-correlated values, following [67]. In the present study, relative errors are less than \(1\%\).

In all comparison figures, Archê results are shown as solid lines with ‘+’ markers, while those obtained with Abinit Extended are represented by ‘o’ circles. Agreement is considered good when the ‘+’ is centered within the circle, corresponding to a relative error of less than 5%. In some configurations (e.g., at low temperatures), Archê cannot match the reference, and the discrepancies can become quite significant (such as an order of magnitude for energies). Because the present study relies solely on a Thomas-Fermi model, agreement is naturally good at high temperatures but deteriorates at lower temperatures, as expected.

Fig. 6 compares the pressures obtained by Archê and Abinit Extended. Fig. 7 compares the energies along the isotherms with the Abinit Extended results. Since energy is defined only up to an additive constant, the Archê energies had to be shifted. We applied a uniform correction equal to the difference between the lowest energy value obtained with Abinit Extended and the corresponding Archê value. This same correction was added to all Archê energies. Visible differences appear at low temperatures, but for temperatures above 10 eV (\(\sim 10^5\) K), the energies obtained are close (within 1%).

In addition to aluminum, boron and iron were also compared with Abinit Extended [68], [69], yielding similar results: good agreement at high temperatures, but reduced accuracy at lower temperatures. Incorporating a gradient correction and, most notably, the KSDT exchange-correlation functional significantly improves agreement at low temperatures. A more in-depth analysis of these results will be presented in a future paper.

5 Benchmarks↩︎

Computations were performed on a TGCC supercomputer [70] located in Bruyères-le-Châtel, France. Each node is equipped with two AMD Milan EPYC 7763 processors (64 cores each),
256 GB of memory, and four NVIDIA A100
(80 GB) GPU accelerators.

5.1 CPU benchmarks↩︎

The following figures illustrate how the total simulation time varies as a function of specific physical parameters (temperature, density, and number of nuclei \(N_i\)) as well as numerical parameters (number of grid points \(N_g\) and number of CPU cores). These benchmarks provide insights into the algorithmic complexity and scalability of the code. The algorithmic complexity of the various code components is presented in Fig.8. The most computationally intensive component is the SCF loop, which includes the FFT resolution of the Poisson equation that has a complexity of \(\mathcal{O}(N_g \log N_g)\) and the evaluation of the Fermi-Dirac integrals across the \(N_g\) grid points. Both procedures are executed during each SCF iteration. Consequently, the total number of SCF iterations significantly contributes to the overall computation time. Outside the SCF loop, the other time-consuming components are the calculations of the ion-ion (\(ii\)) and ion-electron (\(ie\)) forces, which have complexities of \(\mathcal{O}(N_i^2)\) and \(\mathcal{O}(N_i N_g)\), respectively.

5.1.1 Scalability↩︎

Fig. 9 displays the total simulation time as a function of the number of CPU cores under typical simulation conditions, spanning various temperatures and numbers of atoms. Scalability remains nearly ideal up to 512 cores, but decreases slightly beyond that point. This highlights the communication overhead affecting computation time, as parallelization efficiency drops when the number of subdomains increases.

Fig. 10 plots the total simulation time against the number of grid discretization points \(N_g\) for different temperatures. These simulations were performed exclusively to measure performance, as the SCF convergence for the electronic density \(n_e\) under these thermodynamic conditions had already been achieved with smaller grids (\(64\times64\times64 = 2^{18}\)). Since the complexity of the dominant algorithm in this simulation is governed by the Fourier transform \(\mathcal{O}(N_g \log N_g)\), the computation time should theoretically increase slightly faster than linearly with the number of grid points, ignoring communication overhead. However, because the FFT implementation is highly optimized, the execution time appears to scale linearly with the number of grid points for the sizes tested, above \(\approx 2^{24}\).

Fig. 11 presents the total simulation time as a function of the number of nuclei \(N_i\) at different temperatures. The algorithmic complexity of the SCF loop is independent of \(N_i\) (see Fig. 8). However, \(N_i\) plays a role in computing positions, velocities, \(\vec{F}^{ie}\) forces, and in initializing the SCF loop. Furthermore, while the calculation of \(\vec{F}^{ii}\) forces is proportional to \(N_i^2\), this cost is negligible for the number of atoms used in these simulations. Consequently, the observed scalability with respect to the number of atoms \(N_i\) is effectively linear.

5.1.2 Trends with physical parameters↩︎

Fig. 12 shows the total simulation time as a function of temperature for various compressions and numbers of atoms. An interesting trend emerges: the hotter the system, the shorter the time required to simulate it. This behavior is precisely the opposite of what is observed in KSDFT simulations, where higher temperatures lead to a greater number of occupied orbitals and, consequently, longer computation times.

In Archê, this reduction in execution time occurs because the self-consistent loop converges more rapidly. In fact, temperature has no direct impact on the number of operations required per iteration. However, at high temperatures, the electron density becomes more uniform than at low temperatures, allowing the algorithm to accurately converge in significantly fewer iterations.

Fig. 13 illustrates the total simulation time as a function of density across different temperatures. As noted previously, the number of SCF iterations dominates the overall runtime. Density has no direct effect on the number of operations per iteration; nonetheless, as density increases, the SCF convergence requires slightly more iterations to be achieved.

5.2 GPU speedup↩︎

To fully exploit the A100 processors, the most computationally demanding parts of Archê were ported to GPUs. In the CPU version, the majority of the runtime is spent on the Fourier transform due to its individual cost and the high number of calls. Recognizing that the CUDA toolkit provides a highly efficient FFT library called
CuFFT, we replaced our CPU library with
NVIDIA’s single-GPU implementation, which generated most of the results presented here. Since then, we have integrated the HeFFTe library [42] using the CuFFT backend. This enables multi-GPU execution and is mandatory when handling grid sizes of one billion points or more.

To maximize GPU acceleration, it is critical to minimize memory transfers between the host (CPU) and the device (GPU). For this reason, when computing the electronic density, all fields associated with the grid are allocated and evaluated directly on the GPU. By keeping all grid arrays on the GPU, memory transfers per timestep are reduced to just two operations: copying the nuclear positions (input) and retrieving the
electron-nucleus forces \(\vec{F}^{ie}\) (output), which amounts to \(2 \times 3 \times N_i\) floating-point numbers. Consequently, the entire SCF loop — including its initialization and the final force calculations, which account for up to 95% of the CPU computation time — is executed on the GPU with minimal data movement.

In the results presented below, while GPU efficiency varies depending on the specific conditions, running Archê on a GPU consistently provides a massive advantage over a CPU. Currently, the primary limitation is the memory capacity of a single GPU. To simulate larger systems requiring larger grids, we employ the multi-GPU version utilizing domain decomposition, where each subdomain is processed by a single GPU.

Figure 14 presents the ratio of CPU to GPU simulation times as a function of temperature for various compressions and system sizes. GPU efficiency slightly decreases at higher temperatures because fewer iterations are required for the electronic density to converge. This relatively increases the weight of the code sections still executed on the CPU, such as the Ewald summation and the updates of particle positions and velocities. Nevertheless, even in the least favorable scenario of this benchmark, a single GPU outperforms 256 CPU cores by a factor of roughly eight.

Fig. 15 shows the ratio of CPU to GPU simulation times as a function of the number of atoms in the simulation box for different temperatures. The general trend suggests a slight decrease in efficiency as the number of atoms \(N_i\) grows. This can be attributed to increased CPU-GPU data transfers and a larger fraction of time spent on CPU-side computations. However, the curve corresponding to simulations at 300 eV does not exhibit a clear, monotonic trend.

Fig. 16 illustrates the strong scaling of the multi-GPU version of the code across different grid sizes. This version introduces inter-GPU communications, which inevitably degrade performance. Because the simulations were run on nodes containing four GPUs each, communications between different nodes penalize performance even further. For this reason, it is always preferable to use a single GPU when its memory is sufficient to hold the simulation grid. The results indicate that scaling for the \(256^3\) grid is disappointing; ideally, a satisfactory scaling curve would resemble Fig. 9. The explanation is straightforward: GPUs excel at processing massive blocks of data, whereas inter-GPU communications — especially inter-node ones — carry a very high latency cost on our architecture.

6 Conclusion↩︎

We have presented Archê, an OFMD code designed to compute equations of state (EOS) for plasmas, offering an optimal balance between accuracy and efficiency.

The core design is detailed alongside the governing equations, which are presented as closely as possible to their implementation in the source code. Unlike other OFMD codes, Archê relies on a self-consistent field (SCF) approach to calculate the electron density, rather than minimizing the free energy via a conjugate gradient method. This choice enabled the implementation of two algorithms that accelerate SCF convergence by up to a factor of six. Complete derivations for both improvements are provided.

The first improvement to the SCF scheme involves using the converged density field from the previous MD timestep to extract an average one-center density profile. This profile is then superposed onto the new nuclear positions to initialize the density for the current SCF cycle.

The second improvement adapts an algorithm originally proposed for average-atom models by More et al. [41] for use in molecular dynamics. At each SCF iteration, the initial and active densities are mixed to accelerate convergence. The mixing proportion is dynamically determined to minimize a second-order approximation of the free energy.

We also detailed the introduction of a norm-conserving pseudopotential, derived from an average atom model, to describe the electron-nucleus interaction. Particular emphasis was placed on the energy correction required to yield an accurate EOS.

Using the simple Thomas-Fermi kinetic functional, the code was validated for aluminum by comparing its results with Kohn-Sham molecular dynamics simulations performed using Abinit Extended [13]. At high temperatures (above 10 eV), the pressure and internal energy values show agreement typically within 1% to 5%. Larger discrepancies appear at lower temperatures; these can be mitigated by incorporating exchange and correlation corrections, as well as by utilizing semi-local or non-local kinetic energy functionals. These refinements will be explored in future work.

Archê’s performance was evaluated across various architectures, including multi-CPU, single-GPU, and multi-GPU setups. For EOS production, the single-GPU configuration proves to be the most efficient, running an order of magnitude faster than a 256-core CPU architecture. While the multi-GPU architecture allows for larger memory allocation, it currently does not offer the same scalability as its multi-CPU counterpart.

Performance trends as a function of both numerical and physical parameters were analyzed and explained in light of the code’s algorithmic complexity. In summary, Archê exhibits overall linear complexity with respect to the number of atoms and the number of real and reciprocal grid points. Execution time is weakly dependent on density; however, interestingly, it decreases as temperature increases — a stark contrast to the behavior of Kohn-Sham orbital-based simulations.

We believe that Archê establishes a flexible software foundation for testing state-of-the-art kinetic energy functionals at finite temperatures. This addresses the need for accurate EOS and transport coefficients.

Declaration of competing interest↩︎

The authors declare that they have no competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. AI and large language models (LLMs) were not systematically used in the writing of this paper, and were strictly excluded from the research process.

Data availability↩︎

Data other than cited are confidential.

Acknowledgments↩︎

We would like to acknowledge Gilles Kluth, Gwenaël Salin, Etienne Jaupard, Mikael Tacu, Vanina Recoules, and Gérald Faussurier for their valuable advice and fruitful discussions.

References↩︎

---

1. Orbital is the term used "to denote a solution of the wave equation for a system of only one electron". A system of \(N\) electrons can be approximated "by assigning the \(N\) electrons to \(N\) different orbitals, where each orbital is a solution of a wave equation for one electron" [2].↩︎