# HFBTHO-AD: Differentiation of a nuclear energy density functional code

Laurent Hascoët, Matt Menickelly, Sri Hari Krishna Narayanan, Jared

O'Neal, Nicolas Schunck, Stefan Wild

► To cite this version:

Laurent Hascoët, Matt Menickelly, Sri Hari Krishna Narayanan, Jared O'Neal, Nicolas Schunck, et al.. HFBTHO-AD: Differentiation of a nuclear energy density functional code. Computer Physics Communications, 2025, 10.1016/j.cpc.2025.109955. hal-05456113

HAL Id: hal-05456113

https://inria.hal.science/hal-05456113v1

Submitted on 13 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# HFBTHO-AD: Differentiation of a nuclear energy density functional code

Laurent Hascoët<sup>a</sup>, Matt Menickelly<sup>b</sup>, Sri Hari Krishna Narayanan<sup>b,*</sup>,
Jared O'Neal<sup>b</sup>, Nicolas Schunck<sup>c</sup>, Stefan M. Wild<sup>d</sup>

<sup>a</sup>INRIA, Sophia-Antipolis, 2004 Route des lucioles, Valbonne, BP 93 06902, France

<sup>b</sup>Argonne National Laboratory, 9700 S. Cass Ave., Lemont, IL, 60439, USA

<sup>c</sup>Lawrence Livermore National Laboratory, 7000 East Ave., Livermore, CA, 94550, USA

<sup>d</sup>Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA, 94720, USA

## Abstract

The HFBTHO code implements a nuclear energy density functional solver to model the structure of atomic nuclei. HFBTHO has previously been used to calibrate energy functionals and perform sensitivity analysis by using derivative-free methods. To enable derivative-based optimization and uncertainty quantification approaches, we must compute the derivatives of HFBTHO outputs with respect to the parameters of the energy functional, which are a subset of all input parameters of the code. We use the algorithmic/automatic differentiation (AD) tool Tapenade to differentiate HFBTHO. We compare the derivatives obtained using AD against finite-difference approximation and examine the performance of the derivative computation.

**Keywords:** Algorithmic differentiation, HFBTHO, Tapenade

**PACS:** Nuclear density functional theory

## 1. Introduction

Nuclear energy density functional (EDF) methods aim to provide complete characterization of properties of atomic nuclei [1]. Similar in spirit to the density functional theory (DFT) of electrons in atoms and molecules [2, 3, 4], it is a computationally efficient reformulation of the quantum many-body problem of A nucleons (protons and neutrons) interacting through intermediate effective potentials such as the Skyrme [5] or Gogny [6] force or their

<footnote>*Corresponding Author</footnote>
relativistic equivalent [7]. Computational implementations of EDF methods are based on either basis expansion methods or lattice techniques. Examples of the former include the HOSPHE [8], HFBTHO [9, 10, 11, 12], HFODD [13, 14, 15, 16, 17, 18, 19, 20], and DIRHB [21] solvers, which expand solutions on the eigenfunctions of the quantum harmonic oscillator; examples of the latter include HFBRAD [22], EV8 [23, 24], SkyAX [25], and solvers that discretize space within a box [26, 27].

Nuclear EDF theory has a few significant differences from electronic DFT [2, 3, 28, 29]. In contrast to the potentials for systems of electrons in atoms or molecules—the Coulomb force, which is of a two-body nature only—the nuclear potential contains in principle two-body, three-body, and up to A-body terms. However, this expansion is not well known. In practice, most nuclear energy density functionals are constructed from empirical models of two-body effective potentials (e.g., the Skyrme and Gogny potentials). In addition, the small number of constituents—protons and neutrons—in atomic nuclei means that beyond mean-field correlations can be significant. The formal consequence is that nuclear EDF methods rely on spontaneous symmetry breaking to effectively incorporate such correlations into the simple EDF scheme—a technique that has also been recently leveraged by ab initio approaches [30, 31, 32, 33, 34, 35]. The practical consequence of symmetry breaking is that the nuclear wavefunction will not obey symmetries of the nuclear Hamiltonian (and will not conserve the associated quantum numbers): nuclei can be deformed, their number of particles may not be constant, and so on. Most of the available EDF solvers differ by which symmetries are built into the code by default. In this work we focus on the HFBTHO solver, which assumes axial and time-reversal symmetry but allows for the breaking of reflection symmetry and particle number.

Because nuclear energy functionals are phenomenological constructs, they depend on a small (on the order of 5–20 depending on the model) number of parameters that must be calibrated on experimental data. These experimental data typically include binding energies of spherical or well-deformed nuclei, proton radii, energy differences between even-even and odd nuclei, the excitation energy of fission isomers or the energy of the giant dipole resonance, etc.; see Chapter 9 of [1] for a discussion. Until about 15 years ago, the objective function encoding the discrepancy between experimental observables and theoretical predictions often included data in simple systems such as nuclear matter or spherical, closed-shell nuclei because of computational limitations. The Universal Nuclear Energy Density Functional (UN-
EDF) collaboration [36] was the first to perform a fully-fledged optimization with an objective function including dozens of deformed nuclei treated the Hartree-Fock-Bogoliubov approximation [37, 38, 39, 40].

Because the outputs of the EDF did not include derivatives of the various observables with respect to the calibration's decision parameters, this prior optimization work used derivative-free optimization [41]. A potential improvement would be to use derivative-based optimization methods [42]; however, using such methods would require derivatives of the outputs with respect to the parameters. Such derivatives could also provide valuable for other higher-level tasks such as sensitivity analysis and uncertainty quantification (UQ) [43, 44].

To formulate the derivatives that must be computed, we consider the $d(= 108)$ experimental values associated with $n_n(= 72)$ different nucleus configurations in the UNEDF0 dataset [37],

$$ \mathcal{D} = \{y_i\}_{i=1}^d, $$

and the weighted least-squares objective function used for the UNEDF0 optimization,

$$ f(\mathbf{x}) = \|\mathbf{r}(\mathbf{x})\|_2^2 = \sum_{i=1}^{d} r_i^2. \qquad (1) $$

Here $\mathbf{r} = (r_1, \dots, r_d) : \mathbb{R}^{n_x} \to \mathbb{R}^d$ is the weighted residual function for $n_x(= 12)$ input parameters, where the individual residuals evaluated at a given parameter point $\mathbf{x} \in \mathbb{R}^{n_x}$ are defined as

$$ r_i(\mathbf{x}) = \frac{y_i - m_i(\mathbf{x})}{w_i}. $$

The $m_i : \mathbb{R}^{n_x} \to \mathbb{R}$ are the model functions, and the $w_i$ are positive weights assigned to each residual; an example setting of the $w_i$ appropriate for the UNEDF0 dataset is provided in [45]. The $d$ many $m_i(\mathbf{x})$ values are computed by executing HFBTHO once for each of the $n_n$ nucleus configurations. Each computed $r_i(\mathbf{x})$ quantifies the normalized deviation between the theoretical value, $m_i(\mathbf{x})$, and the observable value, $y_i$, in $\mathcal{D}$. Thus $f(\mathbf{x})$ provides a scalar value that expresses how well the theory fits the data at the parameter point $\mathbf{x}$.

To enable derivative-based optimization, we construct the Jacobian function $\mathbf{J}$ that returns the Jacobian matrix $\mathbf{J}(\mathbf{x}) \in \mathbb{R}^{d \times n_x}$ of $\mathbf{r}$ at $\mathbf{x}$. Based on
the least-squares structure of the objective function $f$, we can also determine immediately the gradient of the objective function at $\mathbf{x}$ using

$$ \nabla f(\mathbf{x}) = 2\mathbf{J}(\mathbf{x})^T \mathbf{r}(\mathbf{x}). $$

Numerical differentiation methods, such as finite differencing, estimate derivatives using just the original function code. For example, to numerically approximate the partial derivatives of $\mathbf{J}$, we can use forward differences to compute the elements of $\mathbf{J}$ via

$$ J_{ij}(\mathbf{x}) = \frac{\partial r_i(\mathbf{x})}{\partial x_j} \approx \frac{r_i(\mathbf{x} + h_{ij}\mathbf{e}_j) - r_i(\mathbf{x})}{h_{ij}}, \quad (2) $$

where $\mathbf{e}_j$ is the unit vector of all zeros except for a one in the $j$th entry and $h_{ij} > 0$ is the difference parameter. The decision to use only a first-order finite-difference (FD) scheme is motivated by the desire to use the fewest number of residual function evaluations possible.

In addition to the truncation error (associated with nonzero higher-order derivatives), the quality of these approximate derivatives can suffer from round-off noise in the EDF computation. The composition of round-off with finitely terminated iterative procedures (e.g., for eigenvalue approximation) within a residual function calculation results in deterministic computational noise, which can make obtaining reliable finite-difference approximations difficult in practice. If $h_{ij}$ is chosen too small, then the signal of the difference is lost to the noise. On the other hand, if $h_{ij}$ is chosen too large, then the finite-difference estimate fails to approximate the derivative well (due to the nonzero higher-order derivatives). Techniques for choosing $h_{ij}$ to achieve a balance between noise and truncation error have been established (e.g., EC-Noise [46]) but typically require multiple evaluations of $\mathbf{r}$ at different points $\mathbf{x}$.

For large complex codes such as HFBTHO, analytical computation of derivatives can be a time-consuming and error-prone task. Moreover, as the code changes, the code that computes the derivatives must also be updated.

An alternative to approximation methods or laborious hand coding is to utilize algorithmic/automatic differentiation (AD) to compute the derivatives [47, 48, 49]. AD computes the derivatives of the outputs of a mathematical function, expressed as source code, with respect to any independent parameters, provided that a differentiable relationship exists between these parameters and the outputs. AD can be implemented at compile time
or at run time. AD relies on the larger computation being broken down
into intrinsic mathematical operations. By differentiating each operation
and combining their partial derivatives using the chain rule of calculus, the
derivatives for the large computation can be computed. If the partial deriva-
tives are combined from the first operation to the last, it gives rise to the
so-called *forward mode* of AD. When the derivatives are combined from the
last to the first, it gives rise to the so-called *reverse mode* of AD. For a
function $g: \mathbb{R}^N \rightarrow \mathbb{R}^M$, the forward mode computes the derivatives at the
cost of $O(N)$ function evaluations due to the direction of combination of par-
tial derivatives of the operations. The reverse mode similarly computes the
derivatives at the cost of $O(M)$ function evaluations. This implies that the
forward mode is more efficient when $N < M$ and the reverse mode is more
efficient otherwise. The reverse mode, also called backpropagation, therefore,
underpins supervised machine learning techniques [50]. AD can use compiler-
based techniques to generate a new source code in languages such as Fortran
and C [51, 52] or can use techniques such as operator overloading in lan-
guages such as C++ [53, 54]. AD has been implemented at the intermediate
representation level inside the LLVM compiler by the tool Enzyme [55]. AD
has also been implemented in Python-based machine learning frameworks
such as JAX [56], PyTorch [57], and TensorFlow [58] as well as for Julia [59].

The main contribution of this work is to show how to efficiently compute
the derivatives of HFBTHO using AD. We use the open-source AD engine
Tapenade [51] to differentiate the HFBTHO code that computes the model
functions $m_i$, thus generating the $m'_i$ that are then composed using the chain
rule of calculus to produce $\mathbf{J}$. Tapenade is a source transformation AD
tool that takes as input the source code of a mathematical function and
produces a new source code that computes the derivatives of the function.
Tapenade can be used to differentiate Fortran77, Fortran90, Fortran95, and C
codes, with partial extensions to parallel dialects such as MPI, OpenMP, and
CUDA. It has already been used to differentiate several large computational
models [60, 61, 62].

The rest of this paper is organized as follows. Section 2 explains the
process followed for differentiating HFBTHO. Section 3 explains our testing
setup, and Section 4 presents the results. Section 5 concludes the paper
with a brief summary and suggestions for future work. Appendix A de-
tails HFBTHO's physics options that have been disabled. Appendix B
details the external library files containing linear algebra computation that
are differentiated. Appendix C presents instructions for reproducing the
AD results.

## 2. Implementation of HFBTHO-AD

We differentiate the version of HFBTHO that was published in [11] and was used, for example, in [63]. The top-level routine to differentiate is `HF-BTHO_DFT_SOLVER`. In this initial effort we have turned off certain optional physics modules to create a simplified version of the code. These include the collective inertia package and the finite-range Gogny package. The details of these packages and other disabled options are presented in Appendix A.

For the simplified model, the input parameters of interest are `RHO_NM-`, `E_NM`, `K_NM`, `ASS_NM`, `LASS_NM`, `SMASS_NM`, `CrDr`, `CpVO`, `CrDJ`, and `CJ` since they are the parameters of the energy Skyrme functional that were calibrated on experimental measurements [37, 38, 39, 40, 63]. In this initial study we used the HFBTHO variables `eresu`, `ala2` as outputs. These arrays contain various global characteristics of the nucleus—such as the total energy and charge radius—and a subset of these values are used to evaluate the objective function (Eq. 1) in optimization-based parameter calibration. Because we have identified 14 inputs (`CrDr`, `CrDJ`, `CpVO`, `CJ` are arrays of length 2) and 52 outputs (`eresu`, `ala2` are arrays of length 50 and 2, respectively), we use the forward mode of AD, which is typically the more efficient mode when the number of outputs exceeds the number of input parameters. We additionally use vector forward mode AD, which computes the derivatives with respect to all input parameters simultaneously.

We integrated Tapenade into the existing build system of HFBTHO by adding makefile rules to invoke Tapenade on the top-level routine. We found that Tapenade mostly worked as expected, except for two issues discussed below.

### 2.1. OpenMP Shared-Memory Parallelism

HFBTHO uses OpenMP loops as shown in Listing 1. Tapenade supports direct differentiation of OpenMP parallel work-sharing loops [64]. Nevertheless, minor enhancements to Tapenade were required to support our syntax. The enhancements are now part of Tapenade releases.

<caption>Listing 1: Outline of OpenMP usage in HFBTHO</caption>

```
!$OMP Parallel Default (None)
!$OMP& SHARED (var1, var2)
!$OMP& PRIVATE (var3, var4)
```
```
tid = OMP_GET_THREAD_NUM()
!$OMP DO SCHEDULE(DYNAMIC)
    Do ib=1,var1
    Computation
    End Do ! ib
!$OMP End Do
!$OMP End Parallel
```

## 2.2. Linear Algebra Computation

HFBTHO uses numerous linear algebra routines provided by the BLAS and LAPACK libraries [65, 66, 67]. Many AD frameworks such as JAX, PyTorch, and TensorFlow, as well as tools such as Enzyme, provide customized derivatives for a subset of these routines. The formulae and code in Fortran 2003 for reverse mode differentiation of BLAS matrix operations and the Cholesky factorization routine from LAPACK is given in [68]. [69] provides a collection of rules for the forward mode and reverse mode differentiation of matrix operations provided by LAPACK. While deriving the forward mode formulae for BLAS and LAPACK and handcoding them is likely to result in the shortest execution times, doing so is out of the scope of this work. Here, as a first step, we have chosen to differentiate through the source code of the two libraries, which is possible since both libraries are written in Fortran and their source code is available online. We did so by first identifying the subset of library files that contain all the routines required by HFBTHO (see Appendix B for the list of files). Our current workflow, which we plan to improve in the future, requires all files to be concatenated into one file. Because these files are in Fortran77 format, we first preprocessed the Fortran77 files to convert the comments and line continuations in these files to match the Fortran90 format used by HFBTHO. The preprocessed files were then added to the list of files to be concatenated and subsequently differentiated.

All but one of the BLAS and LAPACK routines invoked by HFBTHO were able to be differentiated by Tapenade. The remaining routine, DS-YEVR, computes selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix A [70].

It may not be appropriate for an AD tool to differentiate through the source code of DSYEVR because doing so involves differentiating through a finite number of steps in a process that is only fully resolved asymptotically. In particular, the convergence of function values (e.g., for DSYEVR) within a given tolerance does not imply the convergence of the corresponding deriva-
tive values to the same tolerance, which can lead to more inaccurate derivatives [71]<fnref n="1" />. Instead, the function should be differentiated at a higher level of abstraction.

We follow the higher-level rule for forward mode AD of DSYEVR given by [69]. In this formulation (see Equation 3) – which assumes that A has distinct eigenvalues and real eigenvectors – D is defined to be the diagonal matrix of eigenvalues and U is defined to be the matrix with columns given by the corresponding eigenvectors, so that $AU = UD$. D and U are computed by calling DSYEVR itself.

Denote by $\dot{D}$ the derivative of the eigenvalues, denote by $\dot{U}$ the derivative of the eigenvectors, and denote by $\dot{A}$ the derivative of the associated DSYEVR input. Let $X \odot Y$ denote the Hadamard (i.e., elementwise) product of two matrices of the same size, defined by $(X \odot Y)_{ij} = X_{ij}Y_{ij}$, and let I denote an identity matrix. We may then compute $\dot{D}$ and $\dot{U}$ via

$$ \begin{aligned} \dot{D} &= I \odot (U^{-1}\dot{A}U) \\ \dot{U} &= U(F \odot (U^{-1}\dot{A}U)), \end{aligned} \qquad (3) $$

where F is defined elementwise by $F_{ij} = \begin{cases} (d_j - d_i)^{-1} & \text{for } i \neq j \\ 0 & \text{otherwise.} \end{cases}$

We have thus separately coded the forward mode derivatives for DSYEVR and used that code in place of the code generated by Tapenade.

### 3. Testing Methodology

Our tests serve to 1) verify that the derivatives computed by different configurations are consistent and that we are able to reproduce historical results using the Jacobians produced by HFBTHO-AD, and 2) compare the performance of AD and finite-difference approximation.

#### 3.1. Comparison Parameter Points

For validation purposes, we have chosen three different parameter points from the UNEDF0 study [37]. The first is the SLy4 parameterization [73], which was used as the starting point for the optimization-based parameter

<footnote n="1">¹We note that this issue arises not just for eigensolver-based calculations, similar phenomena have been documented for basic linear system solvers [72].</footnote>
<table><caption>Table 1: Software configurations</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td><strong>Intel/LAPACK</strong></td><td></td></tr><tr><td>Netlib-LAPACK v3.12.1</td><td></td></tr><tr><td><strong>Intel/MKL</strong></td><td></td></tr><tr><td>Intel OneAPI MKL v2023.2.0</td><td>Intel OneAPI v2023.2.1 w/OpenMP<br>HDF5 v1.14.3</td></tr><tr><td><strong>Intel/LAPACK/AD</strong></td><td></td></tr><tr><td>Netlib-LAPACK v3.12.1</td><td>Intel OneAPI MPI v2021.10.0</td></tr><tr><td>OpenJDK v17.0.11_9</td><td></td></tr><tr><td>Tapenade (commit d0c20bda)</td><td></td></tr><tr><td><strong>GCC/LAPACK/AD</strong></td><td></td></tr><tr><td>Netlib-LAPACK v3.9.0</td><td>GCC v13.2.0 w/OpenMP</td></tr><tr><td>OpenJDK v17.0.11_9</td><td>HDF5 v1.14.3</td></tr><tr><td>Tapenade (commit d0c20bda)</td><td>MPICH v4.2.0</td></tr></tbody></table>

calibrations. The second, UNEDFnb, is the UNEDF0 local minimizer determined via an optimization that imposed no bound constraints on parameter values, and was reported in Table VII of [45]. The third point, UNEDFpre, is the UNEDF0 local minimizer reported in Table VIII of [45] that was determined via a bound-constrained optimization designed to ensure that the minimizer had parameter values within physically acceptable ranges. For this solution, both `SMASS_NM` and `K_NM` are at a bound constraint.

## 3.2. Software Architecture

Our software architecture consists of an HFBTHO executable at an *inner level* that computes the model functions $m_i$. Each invocation of the executable provides all values computed for a given nucleus configuration at a given *x*. A single evaluation of *r* requires running the HFBTHO executable for the $n_n = 72$ nucleus configurations associated with the dataset. A Python package was written at an *outer level* to evaluate both *r* and *J* for AD-enabled HFBTHO using MPI parallelization to take advantage of the trivial parallelization across the 72 computational tasks.

## 3.3. Testing Environment

All data acquisition and tests were run on multiple nodes of the Improv cluster managed by the Laboratory Computing Resource Center of Argonne National Laboratory [74]. Each node contains two sockets each holding a sin-
gle AMD EPYC 7713 64-core processor, and each core uses a single hardware thread.

### 3.3.1. *Inner-Level Software Configurations*

We have constructed multiple software configurations for the inner level as shown in Table 1. The Intel/LAPACK configuration is used for data acquisition tasks that only need to evaluate *f* and *r* (e.g., ECNoise and forward differences). The Intel/MKL configuration, which was used for derivative-free optimization and sensitivity analysis in previous studies because it was determined to be resource efficient, is used here only for confirming software robustness of the Intel/LAPACK configuration. The Intel/LAPACK/AD configuration is used to evaluate *f*, *r*, and *J*. This configuration exclusively uses the open-source Netlib LAPACK because Tapenade requires access to the source code in order to differentiate it. The GCC/LAPACK/AD configuration is used for confirming robustness of our builds and verifying the output of the Intel/LAPACK/AD configuration.

The Intel/LAPACK/AD and GCC/LAPACK/AD configurations are designed to give the derivatives with respect to the 14 input parameters used in the UNEDF2 configuration as explained in Section 2. Because our comparison parameter points are from the UNEDF0 study, however, the finite differences compute the derivatives only with respect to the 12 parameters calibrated in that study.

For all non-AD configurations, the OpenMP stack size is set to 16 MB. For all AD configurations, the OpenMP stack size is set to 10 GB to accommodate the memory needs of the parallelized sections of the Tapenade-generated differentiation code. Appendix C presents the full details of generating and executing the Intel/LAPACK/AD configuration.

### 3.3.2. *OpenMP Threads per MPI Process*

Strong-scaling tests were run to determine how many OpenMP threads to use per MPI process and how many nodes to run per job. For setting up the non-AD test environment, the residual function *r* was evaluated by using the full software architecture based on Intel/LAPACK with different numbers of OpenMP threads per MPI process at UNEDFpre. For each thread count setup, the walltime needed to run each of the 72 HFBTHO computations was divided by the number of iterations run by HFBTHO to converge to the solution of the computation. This approach allows for direct comparison of the walltime per iteration across all HFBTHO computations, which is
<figure><img src="image_2.png" /><figcaption>Figure 1: Strong-scaling results for the (top) Intel/LAPACK software configuration and the (bottom) Intel/LAPACK/AD software configuration.</figcaption></figure>

shown in the boxplot visualization in Figure 1. The figure also shows the speedup and parallel efficiency derived from the data. We note that this design assumes that the scaling results at this single point are representative of scaling across the parameter space.

Because AD-based HFBTHO computations involve significantly longer walltimes, a simpler version of the scaling test was performed for setting up the AD test environment. In particular, the single HFBTHO computation for computing observables of the $^{160}$Dy deformed ground-state configuration at the UNEDFpre solution was run directly with the Intel/LAPACK/AD HFBTHO executable with different numbers of OpenMP threads. Therefore, this test assumes that the scaling results, displayed in Figure 1, for this single computation are representative of scaling for the full set of computations across the parameter space.

The scaling results in Figure 1 are typical of an application whose performance does not scale close to ideal scaling. The number of OpenMP threads per MPI process for non-AD evaluations (8 threads/process) and for
AD-based evaluations (4 threads/process) was chosen as a balance between decreasing walltimes and making efficient use of the hardware. In addition, these choices correspond, respectively, to similar parallel efficiencies of 78% and 71%, which should improve the ability to compare performance results cleanly. Our choices also emphasize that there can be additional tradeoffs to consider – based on resource availability, time and energy constraints, etc. – when deploying AD-based solutions.

### 3.3.3. MPI Processes per Job

The walltime needed to perform a single computation can vary strongly from one nucleus configuration to the next as well as from one parameter point to the next. As a result, evaluating **r** at a single point with the number of MPI processes set to the number of underlying HFBTHO computations would likely result in poor resource usage due to MPI processes idling. Therefore, significantly fewer MPI processes are typically used to (hopefully) achieve better load balancing and hence better resource usage. Therefore, all non-AD jobs are run using 32 MPI processes across two nodes; for AD-based jobs, 32 MPI processes on one node. This matching of MPI process count ensures a similar load-balancing strategy, which should also improve the ability to compare cleanly performance results such as CPU usage or core-hour resource utilization.

### 3.4. Noise Levels and Difference Parameters

We employ ECNoise [75] to obtain independent scalar magnitude estimates of the deterministic computational noise, $\epsilon_{ij}$, present in $r_i$, restricted to the $j^{th}$ parameter, at a fixed input $\mathbf{x}_0$. Using a technique outlined in [46] that also coarsely approximates the second derivative of $r_i$ with respect to the $j^{th}$ parameter, we obtain individually for each $(i, j)$ pair the forward difference parameter $h_{ij}$ using $\epsilon_{ij}$. This specific choice of $h_{ij}$ is generally found in practice to yield reliable FD approximations, and theoretical results concerning its resulting approximation error have been established under particular assumptions on the nature of deterministic noise [46]. In general, computing all $h_{ij}$ at a single $\mathbf{x}_0$ via this enumerative approach requires a constant (on the order of 10) times the dimension ($n_x$) number of evaluations of **r** in a neighborhood about $\mathbf{x}_0$. We, therefore, expect to see the AD-based approach perform better than the FD-based approach (which includes the determination of $h_{ij}$), underscoring the former's value, especially as $n_x$ grows.
<figure><img src="image_3.png" /><figcaption>Figure 3: Absolute deviations between the residual vectors **r**, obtained at UNEDFpre as computed using (top) Intel/LAPACK/AD and GCC/LAPACK/AD, (middle) Intel/LAPACK/AD and Intel/LAPACK, and (bottom) Intel/LAPACK/AD and Intel/MKL.</figcaption></figure>

work. We instead set an internal ECNoise tolerance sufficiently high such that all $h_{ij}$ are determined without failures. While some of the $h_{ij}$ are therefore potentially inaccurate, the agreement between AD and FD-based Jacobian matrices is still reasonably good, as shown in subsequent sections.

# 4. Results

In this section we present validation and performance results.

## 4.1. *Output Comparison*

In Figure 3 we compare the deviations in **r** computed by Intel/LAPACK, Intel/LAPACK/AD, and GCC/LAPACK/AD at UNEDFpre. The devia-
tions are all below $10^{-6}$ and indicate that applying Tapenade did not introduce any considerable errors in the residual computation. Similar results are obtained at UNEDFnb and SLy4.

## 4.2. Validation of Derivatives

We now compare the Intel/LAPACK/AD and GCC/LAPACK/AD Jacobian matrices computed at UNEDFpre, UNEDFnb, and SLy4. At each point we compare the relative deviation between the $108 \times 12 = 1296$ elements of $\mathbf{J}$ computed by the two configurations. Figure 4 (left column) shows that computed derivatives of the two configurations match well at the three parameter points. This result indicates that our approach is robust.

Next we compared the values of $\mathbf{J}$ computed by Intel/LAPACK/AD and FD using Intel/LAPACK at the three parameter points. The $h_{ij}$ for FD are determined by the approach explained in Section 3.4. Figure 4 (right column) shows that computed derivatives of the two configurations match well. There are, however, a small number of elements of $\mathbf{J}$ in each case where the relative deviation is greater than $10^{-2}$. The mismatch is likely due to the presence of abnormal computational noise, as highlighted in Section 3.4.

In Figure 5, we plot the absolute deviations in the entries of Intel/LAPACK/AD and Intel/LAPACK/FD ($|J_{ij}^{AD} - J_{ij}^{FD}|$) against the square root of the product of $\epsilon_{ij}$ with a quantity $\mu_{ij}$. This quantity $\mu_{ij}$ is the previously mentioned coarse estimate of the second derivative $\partial_{jj}r_i$ that is computed as part of the ECNoise procedure for determining $h_{ij}$. Although the conditions of [46, Theorem 3.2] are not precisely satisfied, this result suggests that the true Jacobian $\mathbf{J}$ ought to satisfy, up to multiplicative constants, $|J_{ij} - J_{ij}^{FD}| \approx \sqrt{\epsilon_{ij}\mu_{ij}}$. The trend in Figure 5 suggests that $|J_{ij}^{AD} - J_{ij}^{FD}| \approx \sqrt{\epsilon_{ij}\mu_{ij}}$, lending further credence to the claim that $\mathbf{J}^{AD}$ is a good approximation to $\mathbf{J}$.

Lastly, we examine the computed derivatives when used for sensitivity analysis. At each optimization solution $\hat{\mathbf{x}}$, the UNEDF0 study used second-order central differences to approximate the Jacobian matrix $\mathbf{J}(\hat{\mathbf{x}})$. The matrix was used to compute different quantities such as the statistical uncertainties associated with each coordinate in the solution, which are the diagonal elements of the covariance matrix [76],

$$ \frac{f(\hat{\mathbf{x}})}{n_x - d} (\mathbf{J}(\hat{\mathbf{x}})^T \mathbf{J}(\hat{\mathbf{x}}))^{-1} . $$
<figure><img src="image_4.png" /><figcaption>Figure 4: Histograms of relative deviations of all partial derivatives in the Jacobian matrices obtained at the three comparison parameter points between (left) Intel/LAPACK/AD and GCC/LAPACK/AD and (right) between Intel/LAPACK/AD and FD using Intel/LAPACK. All relative deviations were computed as $|1.0 - J_{\text{other }i,j}/J_{\text{reference }i,j}|$ using the Intel/LAPACK/AD as the reference value.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: Comparison of the absolute deviations between the Intel/LAPACK/AD Jacobian and the Intel/LAPACK/FD Jacobian approximation to a theoretical error approximation [46] derived from the ECNoise intermediate quantities $\epsilon_{ij}$ and $\mu_{ij}$.</figcaption></figure>

Table 2 presents a modified version of the uncertainties published in [37] along with the same uncertainties computed using the Jacobian matrices reported here. Despite HFBTHO undergoing several improvements since the original UNEDF0 study and determining good forward and central difference parameters being inherently difficult, the new results reproduce the original results well.

### 4.3. Performance Analysis

Tables 3 and 4 present performance data for evaluating the objective and Jacobian functions at different parameter points using the two principal software configurations, Intel/LAPACK and Intel/LAPACK/AD, used in this study. The number of MPI processes and OpenMP threads used to gather performance data for these configurations is based on the analysis in Sections 3.3.2 and 3.3.3. For approximating Jacobian matrices using FD, we do not include the time taken to use ECNoise to determine appropriate $h_{ij}$ values. Because of the present implementation of AD in HFBTHO, AD-based performance results are penalized for the computation of derivatives with respect to 14 parameters rather than for only the 12 parameters that are free for UNEDF0.

Table 3 presents the walltimes required for function evaluations. The average FD runtime ratio of $T_{FD}/T_r = 15.47$ is larger than the $(n_x + 1 = 13)$ that we would expect, which could be due to HFBTHO computations
<table><caption>Table 3: Walltimes (in minutes) for evaluating the residual function $\mathbf{r}(\mathbf{x})$ using Intel/LAPACK ($T_r$), evaluating $\mathbf{r}(\mathbf{x})$ using Intel/LAPACK, and approximating $\mathbf{J}(\mathbf{x})$ using Intel/LAPACK with forward differences given a set of $h_{ij}$ ($T_{FD}$), and simultaneously evaluating $\mathbf{r}(\mathbf{x})$ while approximating $\mathbf{J}(\mathbf{x})$ with using Intel/LAPACK/AD ($T_{AD}$). The latter two columns provide the ratio of the (approximate) derivative evaluation walltime to the residual evaluation walltime.</caption><thead><tr><th>x</th><th>T<sub>r</sub></th><th>T<sub>FD</sub></th><th>T<sub>AD</sub></th><th>T<sub>FD</sub>/T<sub>r</sub></th><th>T<sub>AD</sub>/T<sub>r</sub></th></tr></thead><tbody><tr><td>SLy4</td><td>0.70</td><td>10.77</td><td>12.20</td><td>15.50</td><td>17.55</td></tr><tr><td>UNEDFnb</td><td>0.85</td><td>13.00</td><td>15.53</td><td>15.24</td><td>18.20</td></tr><tr><td>UNEDFpre</td><td>1.12</td><td>17.60</td><td>20.12</td><td>15.67</td><td>17.91</td></tr></tbody></table>

<table><caption>Table 4: Resource utilization in core-hours needed for evaluating the residual function $\mathbf{r}(\mathbf{x})$ without AD ($U_r$), evaluating $\mathbf{r}(\mathbf{x})$ without AD and approximating $\mathbf{J}(\mathbf{x})$ with forward differences given a set of $h_{ij}$ ($U_{FD}$) with respect to 12 parameters, and simultaneously evaluating $\mathbf{r}(\mathbf{x})$ while approximating $\mathbf{J}(\mathbf{x})$ with AD ($U_{AD}$) with respect to 14 parameters.</caption><thead><tr><th>x</th><th>U<sub>r</sub></th><th>U<sub>FD</sub></th><th>U<sub>AD</sub></th><th>U<sub>FD</sub>/U<sub>r</sub></th><th>U<sub>AD</sub>/U<sub>r</sub></th></tr></thead><tbody><tr><td>SLy4</td><td>2.97</td><td>45.95</td><td>26.03</td><td>15.50</td><td>8.78</td></tr><tr><td>UNEDFnb</td><td>3.64</td><td>55.47</td><td>33.13</td><td>15.24</td><td>9.10</td></tr><tr><td>UNEDFpre</td><td>4.79</td><td>75.09</td><td>42.92</td><td>15.67</td><td>8.95</td></tr></tbody></table>

requires 8.94 times the resources required for evaluating $\mathbf{r}$ with Intel/LAPACK. Note that the FD-based resource requirement ratios are equal to the FD-based runtime ratios because both $T_r$ and $T_{FD}$ were acquired using the same Intel/LAPACK hardware setup, which used two nodes. However, the AD-based resource requirement ratios are half that of its runtime ratios because $T_{AD}$ was acquired using the Intel/LAPACK/AD hardware setup, which used only one node. Therefore, while the FD differentiation technique is superior in terms of walltime for these test setups, the AD-based differentiation technique is superior in terms of resource usage since less hardware is used to achieve the same parallel efficiency as for the FD-based technique.

## 5. Conclusion and Future Work

We have used the AD tool Tapenade to differentiate the EDF solver HFBTHO. We have differentiated through the BLAS and LAPACK routines that are invoked by HFBTHO (except for DSYEVR). We have computed
derivatives at three parameter points from the UNEDF0 study and validated them using finite differences and successfully reproduced statistical uncertainties of the UNEDFpre solution. Our performance comparison shows that within the context of our particular test environments the AD and FD differentiation techniques have similar walltimes but that AD performance is superior in terms of resource usage. However, FD also requires judicious determination of difference parameters, which adds additional overhead to FD-based approximation and is a key reason for favoring AD over FD.

In future work we will employ the computed derivatives to supply Jaco-
bians to a variation on a Levenberg-Marquardt optimization method for non-
linear least squares, and we will compare its performance with that of prior
derivative-free optimization approaches. We will further study HFBTHO's
noisy behavior, with the aim of understanding and mitigating the abnormal
noise discussed in Section 3.4 and shown in Figure 2, to inform the opti-
mization algorithm. We also will differentiate the latest release of HFBTHO.
This will enable us to perform new studies of nuclear energy functionals – for
example, a comparison of Skyrme and Gogny functionals calibrated on the
same dataset. This will also open the way to performing EDF calibration
and uncertainty quantification while restoring broken symmetries, a process
that we presently believe is computationally too expensive for derivative-free
methods.

Acknowledgment

This work was supported in part by the U.S. Department of Energy, Office
of Science, Office of Advanced Scientific Computing Research Applied Mathe-
matics and SciDAC programs, and Office of Nuclear Physics SciDAC program
under Contract Nos. DE-AC02-06CH11357, DE-AC02-05CH11231, and DE-
AC52-07NA27344. This work was partially performed under the auspices
of the U.S. Department of Energy by Lawrence Livermore National Labora-
tory under contract DE-AC52-07NA27344, managed by Lawrence Livermore
National Security, LLC. We gratefully acknowledge the computing resources
provided on Improv and Bebop, high-performance computing clusters oper-
ated by the Laboratory Computing Resource Center at Argonne National
Laboratory. We thank Paul Hovland, Jan Hückelheim, Rodrigo Navarro
Pérez, and Lars Zurek.
# Appendix A. Physics Modules and Options That Are Disabled

The following is a list of modules that are disabled in order to simplify the code for AD. These modules are not relevant in the optimization of the parameters of Skyrme-like energy functionals.

* **Collective inertia package**: calculates the collective inertia mass tensor at the perturbative cranking approximation. This package is relevant only in fission calculations.

* **Finite-range Gogny package**: computes the matrix elements and expectation value of the finite-range Gogny force. This package is inactive when working with Skyrme-like energy functionals.

* **Particle number projection module**: restores particle number by employing projection techniques. This option is not currently used when performing the calibration of energy functionals.

* **Transformed harmonic oscillator basis package**: performs a local transformation of harmonic oscillator basis functions to improve the asymptotic behavior at large distances. This package is relevant only for describing the physics of weakly bound nuclei.

The following is a list of Makefile preprocessor options that are also disabled in order to simplify the code for AD. They activate/deactivate sections of the code that are not relevant for the optimization of Skyrme-like energy functionals.

* **USE MPI**: enables MPI. Both in direct optimization or Bayesian calibration of the energy functional, HFBTHO runs on a single node with shared-memory parallelism only, making this option inactive.

* **DRIP LINES**: when enabled, calculates an entire mass table from the proton to the neutron dripline. This option is used only in production runs and is not relevant for the calibration of energy functionals.

* **DO MASSTABLE**: when enabled, calculates a section of the mass table. Like DRIP LINES, this option is designed only for production runs and is not relevant for the calibration of energy functionals.
* **DO_PES**: when enabled, calculates a potential energy surface for a set of nuclei. This option is mostly used in large-scale calculations of fission and is not relevant in the calibration of energy functionals.

* **READ_FUNCTIONAL**: when enabled, reads the parameters of the functional from a file. Otherwise, the energy functional is defined in the code based on the value of the input keyword "functional." In the optimization of Skyrme-like energy functionals, the parameters of the functional are passed directly to HFBTHO through an external wrapper, making this option irrelevant.

* **USE_LOCALIZATION**: when enabled, computes the nuclear localization functions. This is relevant only for production runs.

* **GOGNY_SYMMETRIES**: assumes several symmetries in the finite-range matrix elements defining the Gogny force. This option is not relevant when working with Skyrme-like energy functionals.

* **GOGNY_HYPER**: when enabled, uses hypergeometric function to calculate the matrix elements of the finite-range component of the Gogny force. It is not relevant for work with Skyrme-like energy functionals.

* **USE_QRPA**: when enabled, produces output for the QRPA-pnFAM code. This option is not relevant for calculations performed at the HFB level.

# Appendix B. Linear Algebra Files

The linear algebra files used by HFBTHO from LAPACK v3.12.1 are given in Table B.5, Table B.6, Table B.7, and Table B.8. The files in Table B.8 cannot be differentiated as is. They require changes to the code because of limitations of Tapenade, because of the way the code is used in HFBTHO, or because differentiating through the function is not appropriate.

# Appendix C. Reproducing Results

Provided here is a brief description of creating the inner-level INTEL/LAPACK/AD configuration and executing it. The GCC/LAPACK/AD configuration requires a similar process.
<table><caption>Table B.5: LAPACK Files from ./SRC/</caption><tr><td>dlansy.f</td><td>dsytrd.f</td><td>dsterf.f</td><td>dstemr.f</td><td>dormtr.f</td><td>dstebz.f</td></tr><tr><td>dstein.f</td><td>dlasyf.f</td><td>dsytf2.f</td><td>dgetrf2.f</td><td>dlaswp.f</td><td>dsyevr.f</td></tr><tr><td>dsyevx.f</td><td>dsytrf.f</td><td>dpotrf.f</td><td>dsytri.f</td><td>dpotri.f</td><td>dgetrf.f</td></tr><tr><td>dgetri.f</td><td>dtrtri.f</td><td>dlagtf.f</td><td>dlagts.f</td><td>disnan.f</td><td>dlascl.f</td></tr><tr><td>dtrti2.f</td><td>dlacpy.f</td><td>dlarrc.f</td><td>dlae2.f</td><td>dlaebz.f</td><td>dlaev2.f</td></tr><tr><td>dlaisnan.f</td><td>dlaneg.f</td><td>dlanst.f</td><td>dlapy2.f</td><td>dlar1v.f</td><td>dlarf.f</td></tr><tr><td>dlarfb.f</td><td>dlarfg.f</td><td>dlarft.f</td><td>dlarnv.f</td><td>dlarrv.f</td><td>dlarra.f</td></tr><tr><td>dlarrb.f</td><td>dlarrd.f</td><td>dlarre.f</td><td>dlarrf.f</td><td>dlarrj.f</td><td>dlarrk.f</td></tr><tr><td>dlarrr.f</td><td>dlaruv.f</td><td>dlaset.f</td><td>dlasq2.f</td><td>dlasq3.f</td><td>dlasq4.f</td></tr><tr><td>dlasq5.f</td><td>dlasq6.f</td><td>dlasrt.f</td><td>dlatrd.f</td><td>dorm2l.f</td><td>dorm2r.f</td></tr><tr><td>dormql.f</td><td>dormqr.f</td><td>dsytd2.f</td><td>ieeeck.f</td><td>iladlc.f</td><td>iladlr.f</td></tr><tr><td>ilaenv.f</td><td>iparmq.f</td><td>xerbla.f</td><td>dcombssq.f</td><td>dlarf1f.f</td><td>dlarf1f.f</td></tr><tr><td>dlarf11.f</td><td>dlassq.f90</td><td></td><td></td><td></td><td></td></tr></table>

<table><caption>Table B.6: BLAS Files from ./SRC/BLAS/</caption><tr><td>ddot.f</td><td>dswap.f</td><td>dscal.f</td><td>dcopy.f</td><td>daxpy.f</td><td>dgemm.f</td></tr><tr><td>dgemv.f</td><td>dtrsm.f</td><td>dsymv.f</td><td>idamax.f</td><td>dsyr.f</td><td>dsyr2.f</td></tr><tr><td>dsyr2k.f</td><td>dtrmm.f</td><td>dtrmv.f</td><td>dger.f</td><td>dgemmtr.f</td><td>dnrm2.f90</td></tr></table>

<table><caption>Table B.7: Install Files from ./INSTALL/</caption><tr><td>lsame.f</td><td>dlamch.f</td></tr></table>

<table><caption>Table B.8: Handcoded Files</caption><tr><td>dlagts.f</td><td>Rewrite to remove use of three arguments MAX</td></tr><tr><td>dlanst.f</td><td></td></tr><tr><td>dlarfg.f</td><td>Add precision from HFBTHO</td></tr><tr><td>dsyevr.f</td><td>Differentiated as explained in Section 2</td></tr></table>

## Appendix C.1. Obtaining HFBTHO UNEDF0

We have created a snapshot of HFBTHO at https://doi.org/10.5281/zenodo.16249941, to ensure reproducibility. Download the file hfbtho.tar.gz and then unpack it.

$ tar zxvf hfbtho.tar.gz
## Appendix C.2. Installing Tapenade

We detail below the instructions for Linux, but the latest instructions for many operating systems can be found at https://tapenade.gitlabpages.inria.fr/tapenade/distrib/README.html.

1. Clone the code

```
$ cd; git clone https://gitlab.inria.fr/tapenade/tapenade.git
$ cd tapenade
```

2. Use the version of Tapenade from Table 1. This step is optional. As underlying dependencies of Tapenade change, it may be easier to build the latest version of Tapenade.

```
$ git checkout
d0c20bdae18ad1c12e2578103cc28ce5e4995570
```

3. On Linux, depending on your distribution, Tapenade may require you to set the shell variable `JAVA_HOME` to your Java installation directory. It is often `JAVA_HOME=/usr/java/default`. You might also need to modify the `PATH` by adding the bin directory from the Tapenade installation. Every time you wish to use the AD capability with Tapenade, you must first ensure that this Java environment is loaded and that all necessary environment variables are properly set. We have used `OpenJDK 17.0.11_9` here to build Tapenade.

4. Build Tapenade

```
$ ./gradlew clean ; ./gradlew
```

You should now have a working copy of Tapenade. For more information on the `tapenade` command and its arguments, type the following.

```
$ tapenade -?
```

## Appendix C.3. Installing LAPACK

We detail below the instructions for building and installing LAPACK v3.12.1 with Intel's ifort compiler. The instructions for using other compilers or LAPACK v3.9.0 are similar.

* Obtain LAPACK
```
$ cd; wget https://github.com/Reference-LAPACK/lapack/archive/
v3.12.1.tar.gz
$ tar zxvf v3.12.1.tar.gz
```

* Build LAPACK

```
$ cd lapack-3.12.1
$ cp INSTALL/make.inc.iftort ./make.inc
$ make
```

* Build the additional target `la_xisnan.mod`

```
$ cd SRC
$ rm la_xisnan.o
$ make la_xisnan.o
```

If the build products are moved to a separate installation location, make certain to include the manually built `la_xisnan.mod` file in the installation.

### Appendix C.4. Other Dependencies

The dependencies listed below must be installed.

1. Python 3.X

2. HDF5

In case the HFBTHO build system does not find `hdf5.mod` using `pkg-config`, locate `hdf5.mod` and add `-I/location/of/hdf5.mod` to the `OPTIONS_FC` flag HFBTHO-AD's Makefile.

### Appendix C.5. Building and Running HFBTHO-AD

Running a derivative simulation with HFBTHO using Tapenade is not very different from running any other typical HFBTHO simulation and broadly involves the same steps.

* Set up environment variables

```
$ export TAPENADE_HOME=/Location/of/Tapenade
$ export TAPENADE_LAPACK_SRC=/Location/of/LAPACK
$ export TAPENADE_LAPACK_INSTALLATION=\
${TAPENADE_LAPACK_SRC}/SRC/
```
* Go to the `src` directory of HFBHTO. From the root of the HFBHTO repository, the path is as follows.

```
$ cd src
```

* Edit `Makefile` to set variables.

```
COMPILER     = IFORT
FORTRANICI = ifort
USE_OPENMP   = 1
AD           = 1
AD VECTOR MODE = 1
```

* Clean any remnant files in your build subdirectory from a previous simulation.

```
$ make clean
```

* Run the make commands. The results in Section 4 are derived by using the forward mode with independent and dependent variables that are specified in the file `src/hfbtho_ad/MakefileTap`.

```
$ make
```

These commands will generate the executable `hfbtho_ad/hftbho_main`.

* Run the executable, and examine the output.

If `OpenMP` is used, do the following first.

```
$ export OMP_STACKSIZE="10 G"; export OMP_NUM_THREADS=32;
```

Then run the code using the provided input files.

```
$ cd hfbtho_ad;
$ ./hfbtho_main -f hfbtho_UNEDFO_Z20_N20_sphGS.dat \
  -u unedf_UNEDFO_Z20_N20_sphGS.dat > hfbtho_out.txt 2>&1
```

* Analyze the output written to `hfbtho_out.txt` as well as final `r` values written in `eresu.csv` and `ala2.csv`. The derivatives are written to `eresu_d.csv` and `ala2_d.csv`.
References

[1] N. Schunck, Energy Density Functional Methods for Atomic Nuclei, IOP Expanding Physics, IOP Publishing, Bristol, UK, 2019. doi:10.1088/2053-2563/aae0ed.

[2] R. Parr, W. Yang, Density Functional Theory of Atoms and Molecules, International Series of Monographs on Chemistry, Oxford University Press, New York, 1989. doi:10.1093/oso/9780195092769.001.0001.

[3] R. Dreizler, E. Gross, Density Functional Theory: An Approach to the Quantum Many-Body Problem, Springer-Verlag, 1990. doi:10.1007/978-3-642-86105-5.

[4] R. Eschrig, Fundamentals of Density Functional Theory, Teubner, Leipzig, 1996. doi:10.1007/978-3-322-97620-8.

[5] J. R. Stone, P.-G. Reinhard, The Skyrme interaction in finite nuclei and nuclear matter, Prog. Part. Nucl. Phys. 58 (2007) 587. doi:10.1016/j.ppnp.2006.07.001.

[6] L. M. Robledo, T. R. Rodríguez, R. R. Rodríguez-Guzmán, Mean field and beyond description of nuclear structure with the Gogny force: A review, J. Phys. G: Nucl. Part. Phys. 46 (2019) 013001. doi:10.1088/1361-6471/aadebd.

[7] T. Nikšić, D. Vretenar, P. Ring, Relativistic nuclear energy density functionals: Mean-field and beyond, Prog. Part. Nucl. Phys. 66 (2011) 519. doi:10.1016/j.ppnp.2011.01.055.

[8] B. G. Carlsson, J. Dobaczewski, J. Toivanen, P. Veselý, Solution of self-consistent equations for the N3LO nuclear energy density functional in spherical symmetry. The program HOSPHE (v1.02), Comput. Phys. Commun. 181 (2010) 1641. doi:10.1016/j.cpc.2010.05.022.

[9] M. V. Stoitsov, J. Dobaczewski, W. Nazarewicz, P. Ring, Axially deformed solution of the Skyrme-Hartree-Fock-Bogolyubov equations using the transformed harmonic oscillator basis. The program HFBTHO (v1.66p), Comput. Phys. Commun. 167 (2005) 43. doi:10.1016/j.cpc.2005.01.001.
[10] M. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, S. Wild, Axially deformed solution of the Skyrme-Hartree-Fock-Bogoliubov equations using the transformed harmonic oscillator basis (II) HFBTHO v2.00d: A new version of the program, Comput. Phys. Commun. 184 (2013) 1592. doi:10.1016/j.cpc.2013.01.013.

[11] R. N. Perez, N. Schunck, R.-D. Lasseri, C. Zhang, J. Sarich, Axially deformed solution of the Skyrme-Hartree-Fock-Bogolyubov equations using the transformed harmonic oscillator basis (III) HFBTHO (v3.00): A new version of the program, Comput. Phys. Commun. 220 (2017) 363. doi:10.1016/j.cpc.2017.06.022.

[12] P. Marević, N. Schunck, E. M. Ney, R. Navarro Pérez, M. Verriere, J. O'Neal, Axially-deformed solution of the Skyrme-Hartree-Fock-Bogoliubov equations using the transformed harmonic oscillator basis (IV) HFBTHO (v4.0): A new version of the program, Comput. Phys. Commun. 276 (2022) 108367. doi:10.1016/j.cpc.2022.108367.

[13] J. Dobaczewski, J. Dudek, Solution of the Skyrme-Hartree-Fock equations in the Cartesian deformed harmonic oscillator basis I. The method, Comput. Phys. Commun. 102 (1997) 166. doi:10.1016/S0010-4655(97)00004-0.

[14] J. Dobaczewski, J. Dudek, Solution of the Skyrme-Hartree-Fock equations in the Cartesian deformed harmonic oscillator basis II. The program HFODD, Comput. Phys. Commun. 102 (1997) 183. doi:10.1016/S0010-4655(97)00005-2.

[15] J. Dobaczewski, J. Dudek, Solution of the Skyrme-Hartree-Fock equations in the Cartesian deformed harmonic-oscillator basis. (III) HFODD (v1.75r): A new version of the program, Comput. Phys. Commun. 131 (2000) 164. doi:10.1016/S0010-4655(00)00121-1.

[16] J. Dobaczewski, P. Olbratowski, Solution of the Skyrme-Hartree-Fock-Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis. (IV) HFODD (v2.08i): A new version of the program, Comput. Phys. Commun. 158 (2004) 158. doi:10.1016/j.cpc.2004.02.003.

[17] J. Dobaczewski, W. Satuła, B. G. Carlsson, J. Engel, P. Olbratowski, P. Powałowski, M. Sadziak, J. Sarich, N. Schunck,
A. Staszczak, M. Stoitsov, M. Zalewski, H. Zduńczuk, Solution of the Skyrme-Hartree-Fock-Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis. (VI) HFODD (v.240h): A new version of the program, Comput. Phys. Commun. 180 (2009) 2361. doi:10.1016/j.cpc.2009.08.009.

[18] N. Schunck, J. Dobaczewski, J. McDonnell, W. Satuła, J. Sheikh, A. Staszczak, M. Stoitsov, P. Toivanen, Solution of the Skyrme-Hartree-Fock-Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis. (VII) HFODD (v2.49t): A new version of the program, Comput. Phys. Commun. 183 (2012) 166. doi:10.1016/j.cpc.2011.08.013.

[19] N. Schunck, J. Dobaczewski, W. Satuła, P. Bączyk, J. Dudek, Y. Gao, M. Konieczka, K. Sato, Y. Shi, X. B. Wang, T. R. Werner, Solution of the Skyrme-Hartree-Fock-Bogolyubov equations in the Cartesian deformed harmonic-oscillator basis. (VIII) HFODD (v2.73y): A new version of the program, Comput. Phys. Commun. 216 (2017) 145. doi:10.1016/j.cpc.2017.03.007.

[20] J. Dobaczewski, P. Bączyk, P. Becker, M. Bender, K. Bennaceur, J. Bonnard, Y. Gao, A. Idini, M. Konieczka, M. Kortelainen, L. Próchniak, A. M. Romero, W. Satuła, Y. Shi, T. R. Werner, L. F. Yu, Solution of universal nonrelativistic nuclear DFT equations in the Cartesian deformed harmonic-oscillator basis. (IX) HFODD (v3.06h): A new version of the program, J. Phys. G: Nucl. Part. Phys. 48 (2021) 102001. doi:10.1088/1361-6471/ac0a82.

[21] T. Nikśic, N. Paar, D. Vretenar, P. Ring, DIRHB - A relativistic self-consistent mean-field framework for atomic nuclei, Comput. Phys. Commun. 185 (2014) 1808. doi:10.1016/j.cpc.2014.02.027.

[22] K. Bennaceur, J. Dobaczewski, Coordinate-space solution of the Skyrme-Hartree-Fock-Bogolyubov equations within spherical symmetry. The program HFBRAD (v1.00), Comput. Phys. Commun. 168 (2005) 96. doi:10.1016/j.cpc.2005.02.002.

[23] P. Bonche, H. Flocard, P.-H. Heenen, Solution of the Skyrme HF+BCS equation on a 3D mesh, Comput. Phys. Commun. 171 (2005) 49. doi:10.1016/j.cpc.2005.05.001.
[24] W. Ryssens, V. Hellemans, M. Bender, P.-H. Heenen, Solution of the Skyrme-HF+BCS equation on a 3D mesh, II: A new version of the Ev8 code, Comput. Phys. Commun. 187 (2015) 175. doi:10.1016/j.cpc.2014.10.001.

[25] P. G. Reinhard, B. Schuetrumpf, J. A. Maruhn, The Axial Hartree-Fock + BCS Code SkyAx, Comput. Phys. Commun. 258 (2021) 107603. doi:10.1016/j.cpc.2020.107603.

[26] J. C. Pei, G. I. Fann, R. J. Harrison, W. Nazarewicz, Y. Shi, S. Thornton, Adaptive multi-resolution 3D Hartree-Fock-Bogoliubov solver for nuclear structure, Phys. Rev. C 90 (2014) 024317. doi:10.1103/PhysRevC.90.024317.

[27] M. Chen, T. Li, B. Schuetrumpf, P.-G. Reinhard, W. Nazarewicz, Three-dimensional Skyrme Hartree-Fock-Bogoliubov solver in coordinate-space representation, Comput. Phys. Commun. 276 (2022) 108344. doi:10.1016/j.cpc.2022.108344.

[28] J. Engel, Intrinsic-density functionals, Phys. Rev. C 75 (2007) 014306. doi:10.1103/PhysRevC.75.014306.

[29] J. Messud, M. Bender, E. Suraud, Density functional theory and Kohn-Sham scheme for self-bound systems, Phys. Rev. C 80 (2009) 054314. doi:10.1103/PhysRevC.80.054314.

[30] T. Duguet, Symmetry broken and restored coupled-cluster theory: I. Rotational symmetry and angular momentum, J. Phys. G: Nucl. Part. Phys. 42 (2015) 025107. doi:10.1088/0954-3899/42/2/025107.

[31] T. Duguet, A. Signoracci, Symmetry broken and restored coupled-cluster theory: II. Global gauge symmetry and particle number, J. Phys. G: Nucl. Part. Phys. 44 (2016) 015103. doi:10.1088/0954-3899/44/1/015103.

[32] J. M. Yao, J. Engel, L. J. Wang, C. F. Jiao, H. Hergert, Generator-coordinate reference states for spectra and $0\nu\beta\beta$ decay in the in-medium similarity renormalization group, Phys. Rev. C 98 (2018) 054311. doi:10.1103/PhysRevC.98.054311.
[33] G. Hagen, S. J. Novario, Z. H. Sun, T. Papenbrock, G. R. Jansen, J. G. Lietz, T. Duguet, A. Tichai, Angular-momentum projection in coupled-cluster theory: Structure of $^{34}$Mg, Phys. Rev. C 105 (2022) 064311. doi:10.1103/PhysRevC.105.064311.

[34] B. S. Hu, Z. H. Sun, G. Hagen, T. Papenbrock, Ab initio computations of strongly deformed nuclei near $^{80}$Zr, Phys. Rev. C 110 (2024) L011302. doi:10.1103/PhysRevC.110.L011302.

[35] Z. H. Sun, A. Ekström, C. Forssén, G. Hagen, G. R. Jansen, T. Papenbrock, Multiscale physics of atomic nuclei from first principles, Phys. Rev. X 15 (2025) 011028. doi:10.1103/PhysRevX.15.011028.

[36] S. Bogner, A. Bulgac, J. Carlson, J. Engel, G. Fann, R. J. Furnstahl, S. Gandolfi, G. Hagen, M. Horoi, C. Johnson, Computational nuclear quantum many-body problem: The UNEDF project, Comput. Phys. Commun. 184 (2013) 2235. doi:10.1016/j.cpc.2013.05.020.

[37] M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, S. Wild, Nuclear energy density optimization, Phys. Rev. C 82 (2010) 024313. doi:10.1103/PhysRevC.82.024313.

[38] M. Kortelainen, J. McDonnell, W. Nazarewicz, P.-G. Reinhard, J. Sarich, N. Schunck, M. V. Stoitsov, S. M. Wild, Nuclear energy density optimization: Large deformations, Phys. Rev. C 85 (2012) 024304. doi:10.1103/PhysRevC.85.024304.

[39] M. Kortelainen, J. McDonnell, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sarich, N. Schunck, S. M. Wild, D. Davesne, J. Erler, A. Pastore, Nuclear energy density optimization: Shell structure, Phys. Rev. C 89 (2014) 054314. doi:10.1103/PhysRevC.89.054314.

[40] S. M. Wild, J. Sarich, N. Schunck, Derivative-free optimization for parameter estimation in computational nuclear physics, J. Phys. G: Nucl. Part. Phys. 42 (2015) 034031. doi:10.1088/0954-3899/42/3/034031.

[41] J. Larson, M. Menickelly, S. M. Wild, Derivative-free optimization methods, Acta Numerica 28 (2019) 287-404. doi:10.1017/S0962492919000060.
[42] J. Nocedal, S. J. Wright, Numerical Optimization, second ed., Springer, 2006. doi:10.1007/978-0-387-40065-5.

[43] R. C. Smith, Uncertainty Quantification: Theory, Implementation, and Applications, SIAM, 2014. doi:10.1137/1.9781611973228.

[44] B. D. Carlsson, A. Ekström, C. Forssén, D. F. Strömberg, G. R. Jansen, O. Lilja, M. Lindby, B. A. Mattsson, K. A. Wendt, Uncertainty analysis and order-by-order optimization of chiral nuclear interactions, Phys. Rev. X 6 (2016) 011019. doi:10.1103/PhysRevX.6.011019.

[45] M. Kortelainen, T. Lesinski, J. J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, S. M. Wild, Nuclear energy density optimization, Physical Review C 82 (2010) 024313. doi:10.1103/PhysRevC.82.024313.

[46] J. J. Moré, S. M. Wild, Estimating derivatives of noisy simulations, ACM Transactions on Mathematical Software 38 (2012) 19:1–19:21. doi:10.1145/2168773.2168777.

[47] A. Griewank, A. Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, SIAM, 2008. doi:10.1137/1.9780898717761.

[48] U. Naumann, The Art of Differentiating Computer Programs, SIAM, Philadelphia, PA, 2011. doi:10.1137/1.9781611972078.

[49] S. H. K. Narayanan, B. Norris, B. Winnicka, ADIC2: Development of a component source transformation system for differentiating C and C++, Procedia Computer Science 1 (2010) 1845–1853. doi:10.1016/j.procs.2010.04.206, ICCS 2010.

[50] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: A survey, Journal of Machine Learning Research 18 (2018) 1–43. URL: http://jmlr.org/papers/v18/17-468.html.

[51] L. Hascoët, V. Pascual, The Tapenade automatic differentiation tool: Principles, model, and specification, ACM Transactions on Mathematical Software 39 (2013) 20:1–20:43. doi:10.1145/2450153.2450158.
[52] R. Giering, T. Kaminski, Recipes for adjoint code construction, ACM Trans. Math. Softw. 24 (1998) 437–474. doi:10.1145/293686.293695.

[53] A. Griewank, D. Juedes, J. Utke, Algorithm 755: ADOL-C: A package for the automatic differentiation of algorithms written in C/C++, ACM Trans. Math. Softw. 22 (1996) 131–167. doi:10.1145/229473.229474.

[54] N. G. M. Sagebaum, T. Albring, High-performance derivative computations using CoDiPack, ACM Trans. Math. Softw. 45 (2019). doi:10.1145/3356900.

[55] W. Moses, V. Churavy, Instead of rewriting foreign code for machine learning, automatically synthesize fast gradients, in: H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, H. Lin (Eds.), Advances in Neural Information Processing Systems, volume 33, Curran Associates, Inc., 2020, pp. 12472–12485. URL: https://proceedings.neurips.cc/paper/2020/file/9332c513ef44b682e9347822c2e457ac-Paper.pdf.

[56] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, Q. Zhang, JAX: composable transformations of Python+NumPy programs, 2018. URL: http://github.com/google/jax.

[57] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala, PyTorch: An imperative style, high-performance deep learning library, in: H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, R. Garnett (Eds.), Advances in Neural Information Processing Systems, volume 32, Curran Associates, Inc., 2019. URL: https://proceedings.neurips.cc/paper_files/paper/2019/file/bdbca288fee7f92f2bfa9f7012727740-Paper.pdf.

[58] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Mur ray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng, TensorFlow:
Large-scale machine learning on heterogeneous systems, 2015. URL: https://www.tensorflow.org/.

[59] W. S. Moses, V. Churavy, L. Paehler, J. Hückelheim, S. H. K. Narayanan, M. Schanen, J. Doerfert, Reverse-mode automatic differentiation and optimization of GPU kernels via Enzyme, in: Proc. Int. Conf. High Perf. Comp., Network., Storage and Analysis, SC ’21, ACM, New York, 2021. doi:10.1145/3458817.3476165.

[60] S. S. Gaikwad, L. Hascoet, S. H. K. Narayanan, L. Curry-Logan, R. Greve, P. Heimbach, SICOPOLIS-AD v2: tangent linear and adjoint modeling framework for ice sheet modeling enabled by automatic differentiation tool Tapenade, Journal of Open Source Software 8 (2023) 4679. doi:10.21105/joss.04679.

[61] B. Maugars, S. Bourasseau, C. Content, B. Michel, B. Berthoul, J. N. Ramirez, P. Raud, L. Hascoët, Algorithmic differentiation for an efficient CFD solver, in: ECCOMAS 2022 - 8th European Congress on Computational Methods in Applied Sciences and Engineering, Oslo, Norway, 2022. URL: https://hal.science/hal-03759125.

[62] S. Carli, L. Hascoët, W. Dekeyser, M. Blommaert, Algorithmic differentiation for adjoint sensitivity calculation in plasma edge codes, Journal of Computational Physics 491 (2023) 112403. doi:10.1016/j.jcp.2023.112403.

[63] N. Schunck, J. O’Neal, M. Grosskopf, E. Lawrence, S. M. Wild, Calibration of energy density functionals with deformed nuclei, J. Phys. G: Nucl. Part. Phys. 47 (2020) 074001. doi:10.1088/1361-6471/ab8745.

[64] J. Hückelheim, L. Hascoët, Source-to-source automatic differentiation of OpenMP parallel loops, ACM Trans. Math. Softw. 48 (2022). doi:10.1145/3472796.

[65] E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users’ Guide, third ed., Society for Industrial and Applied Mathematics, 1999. doi:10.1137/1.9780898719604.
[66] L. S. Blackford, et al., An updated set of basic linear algebra subprograms (BLAS), ACM Trans. Math. Softw. 28 (2002) 135-151. doi:10.1145/567806.567807.

[67] C. 11/2024, LAPACK – Linear Algebra PACKage, 2024. URL: https://www.netlib.org/lapack/index.html.

[68] K. Jonasson, S. Sigurdsson, H. F. Yngvason, P. O. Ragnarsson, P. Meltsted, Algorithm 1005: Fortran subroutines for reverse mode algorithmic differentiation of blas matrix operations, ACM Trans. Math. Softw. 46 (2020). URL: https://doi.org/10.1145/3382191. doi:10.1145/3382191.

[69] M. B. Giles, An extended collection of matrix derivative results for forward and reverse mode algorithmic differentiation, 08/01, 2008. URL: https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf.

[70] C. 11/2024, LAPACK documentation, 2024. URL: https://www.netlib.org/lapack/explore-html/.

[71] B. Christianson, Reverse accumulation and attractive fixed points, Optimization Methods and Software 3 (1994) 311-326. doi:10.1080/10556789408805572.

[72] J. J. Moré, S. M. Wild, Do you trust derivatives or differences?, Journal of Computational Physics 273 (2014) 268-277. doi:10.1016/j.jcp.2014.04.056.

[73] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, R. Schaeffer, New Skyrme effective forces for supernovae and neutron rich nuclei, Physica Scripta 1995 (1995) 231. doi:10.1088/0031-8949/1995/T56/034.

[74] C. 07/2025, Improv, 2025. URL: https://www.lcrc.anl.gov/systems/improv.

[75] J. J. Moré, S. M. Wild, Estimating computational noise, SIAM Journal on Scientific Computing 33 (2011) 1292-1314. doi:10.1137/100786125.

[76] G. A. F. Seber, C. J. Wild, Nonlinear Regression, John Wiley & Sons, Inc., 1989. doi:10.1002/0471725315.
The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory ('Argonne'). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doi-public-access-plan.