OpenBTE: a Solver for ab-initio Phonon Transport in Multidimensional Structures

The increasing accuracy of first-principles calculations of heat transport, along with the growing availability of computing resources, has enabled the systematic prediction of lattice thermal conductivity in semiconductors and two-dimensional (2D) materials [1, 2, 3]. If key theoretical developments have solidified our understanding of phonon dynamics at the nanoscales, the fast progress of this field has been unarguably aided by the release of open-source packages. For example, ShengBTE [4] implements the iterative solution of the Boltzmann transport equation (BTE) [5], with force constants computed by first-principles [6, 7]. The tool can be primarily used for either bulk materials or nanowires [8]. Later iterations include  AlmaBTE [9], which, among several improvements, features a Monte Carlo solver for superlattices, and FourPhonons [10, 11], a tool that, as the name suggests, can handle four-phonon scattering. Another notable package is Phono3Py [12], where the BTE is solved both within the relaxation time approximation (RTA) and including the full scattering operator [13]. Albeit not strictly required, these tools are mostly used in tandem with super-cell approaches. When a density functional perturbation theory (DFPT) is available for a given material, it is often convenient to compute directly the Fourier-transformed force constants, without resorting to supercells; a prominent package in this space is d3q [14], implemented on top of QUANTUM ESPRESSO [15]. This software also includes the possibility of simulating thin films [16] and coherence effects [17]. In a similar effort, the recently released Kaldo[18] combines the Green-Kubo model and the BTE to take into account material disorder [19].

While the dissemination of these software manifests the increasing excitement around phonon related applications, their focus is primarily on bulk materials or simple geometries. In fact, most phonon dynamics simulations in complex geometries are currently performed by MonteCarlo solvers [20, 21], with a recent implementation provided by MCBTE [22]. Thus, an open-source tool that solves space-dependent heat transport deterministically, i.e. in the same spirit as in bulk calculations, is still pending. We fill this gap by introducing OpenBTE, a solver for the first-principles, multidimensional phonon BTE. Key capabilities include the calculation of the temperature and heat flux maps, as well as the mode-resolve effective thermal conductivity. Ab-initio harmonic and anharmonic properties are delegated to external software. While flexible and straightforward to expand, OpenBTE currently implements the anisotropic-mean-free-path BTE (aMFP-BTE) [23], a method based on the interpolation of the phonon populations onto the vectorial MFP space. This method has shown a 50x speed up with respect to the mode-resolve BTE for thick Si membranes, while not compromising on the accuracy. Furthermore, it scales with constant time with respect to the number of phonon branches, therefore enabling fast simulations of size effects in complex-unit-cell materials, such as bismuth telluride.

The paper is organized as follows. First, we review the model and the assumptions therein, then we provide details on the implemented algorithm and the software’s main modules. An example of a porous Si membrane will follow, and final remarks conclude the paper. We also provides technical implementation details in the appendices. Enabling fast simulations of heat transport in nanostructures, OpenBTE sets out to accelerate the development of thermal-related applications, such as thermal management and thermoelectrics. OpenBTE is released under the GPL-3.0 license and currently hosted on GitHub [24].

The steady-state, linearized phonon BTE in the temperature formulation reads [23]

where $\mu,\nu$ collectively indicate phonon branch and wave vector, running up to $N_{p}$ and $N_{q}$, respectively. The unknowns $\Delta T_{\mu}(\mathbf{r})$ are the phonon pseudo-temperatures (referred to hereafter simply as temperatures), connected to the deviational non-equilibrium phonon population via $\Delta T_{\mu}(\mathbf{r})=\Delta n_{\mu}(\mathbf{r})\hbar\omega_{\mu}C_{\mu}^{-1}$. The mode-specific heat capacity is $C_{\mu}=k_{B}\left(\eta_{\mu}/\sinh{\eta_{\mu}}\right)^{2}$, with $\eta_{\mu}=\hbar\omega_{\mu}/k_{B}/T/2$. The terms $\hbar\omega_{\mu}$ and $\mathbf{v}_{\mu}$ are the phonon frequencies and the group velocities, respectively. In this formulation, the thermal flux is given by $\mathbf{J}(\mathbf{r})=\left(N_{q}\mathcal{V}\right)^{-1}\sum_{\mu}\Delta T_{\mu}(\mathbf{r})\mathbf{S}_{\mu}$, where $\mathcal{V}$ is the volume of the unit cell, and $\mathbf{S}_{\mu}=C_{\mu}\mathbf{v}_{\mu}$. The term $W_{\mu\nu}$ is the energy form of the scattering operator [25]. While we plan to make available the implementation of Eq. 1, currently only the relaxation-time-approximation (RTA) is supported; within the RTA, the RHS of Eq. 1 becomes $C_{\mu}/\tau_{\mu}\left[\Delta T_{\mu}(\mathbf{r})-\Delta T(\mathbf{r})\right]$, leading to

where $\tau_{\mu}$ is the relaxation time, and $\Delta T(\mathbf{r})$ is the pseudo lattice temperature; this term is evaluated so that energy conservation is satisfied, i.e. $\nabla\cdot\mathbf{J}(\mathbf{r})=0$, yielding $\Delta T(\mathbf{r})=\left[\sum_{\nu}C_{\nu}/\tau_{\nu}\right]^{-1}\sum_{\mu}C_{\mu}/\tau_{\mu}\Delta T_{\mu}(\mathbf{r})$ [26].

The simulation domain is a cuboid of size $L_{x}$, $L_{y}$ and $L_{z}$ with periodic boundary conditions applied throughout by default. A difference of temperature $\Delta T_{\mathrm{ext}}$ is also enforced along a chosen Cartesian axis. At the walls of the pores, energy conservation leads to $\sum_{\nu}\Delta T_{\nu}\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}}=0$, which can be split into incoming and outgoing flux, i.e. $\sum_{\mu}\Delta T_{\mu}\mathrm{ReLu}(\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}})=-\sum_{\mu}\Delta T_{\mu}\mathrm{ReLu}(-\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}})$; note that we define $\mathrm{ReLu}(x)=x\Theta(x)$, with $\Theta(x)$ being the Heaviside function. The boundary condition is imposed to phonons leaving the surface; generally, we have $\Delta T_{\mu}^{-}=\sum_{\nu}R_{\mu\nu}\Delta T_{\nu}^{+}$, where $-$ ($+$) are for outgoing (incoming) phonons, and $R_{\mu\nu}$ is a mode-resolved bouncing matrix. However, currently only a simplified model is available, with which we assume that phonons thermalize to a single temperature [20]; this choice corresponds to $R_{\mu\nu}=\left[\sum_{k}\mathrm{ReLu}(\mathbf{S}_{k}\cdot\mathbf{\hat{n}})\right]^{-1}\mathrm{ReLu}(\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}})$. Lastly, the effective thermal conductivity, assuming we have imposed the temperature gradient along $x$, reads

Equation 2, when solved iteratively, requires inverting as many linear systems as the number of phonons modes, an endeavour that may easily become prohibitive. To ameliorate the computational effort, we instead implement the recently introduced anisotropic-MFP-BTE (aMFP-BTE), a model based on interpolating the phonon temperatures onto the vectorial MFP space [23]. In practice, instead of solving Eq. 2 $N_{p}N_{q}$ times at each iteration, we solve it for a set of vectorial MFPs, $\mathbf{F}_{ml}$, located on a spherical grid. The indices $m$ and $l$ label the MFP $\Lambda_{m}$ (up to $N_{\Lambda}$) and phonon directions $\mathbf{\hat{s}}_{l}$ (up to $N_{\Omega}$), respectively. As $N_{\Lambda}N_{\Omega}<<N_{p}N_{q}$, while not compromising on accuracy, this scheme results in a much faster runtime with respect to the case with mode-resolved resolution, i.e. no interpolation.

As derived in A, the discretization of aMFP-BTE model in real and momentum space yields the following iterative linear system

where $\mathbf{A}_{ml}$ is the stiffness matrix for a given MFP and direction, $\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}$ is associated to the lattice temperature and adiabatic boundary conditions, and $\mathbf{P}_{ml}$ represents the perturbation. Once Eq. 4 converges, the effective thermal conductivity is computed as $\kappa_{\mathrm{eff}}=\kappa_{\mathrm{bal}}/2+\sum_{ml}\mathbf{K}_{ml}^{T}\cdot\mathbf{\Delta T}_{ml}$, where $\kappa_{\mathrm{bal}}$ is the ballistic thermal conductivity, and $\mathbf{K}_{ml}$ is associated to the thermal flux.

The first guess to Eq. 4 is given by the standard heat conduction equation. While it is straightforward to discretize the diffusive equation in orthogonal grids, it becomes convoluted for unstructured meshes. In OpenBTE, we use the approach described in [27], where the non-orthogonal contribution of the flux is treated explicitly, i.e. it depends on previous solutions. The resulting system, analogously to Eq. 4, is an iterative linear system,

where $\mathbf{S}$ arises from the non-orthogonality of the mesh. Note that $\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}$ in Eq. 4 and $\mathbf{S}$ in Eq. 5 are not related. The same applies with the terms sharing the same basename in these two equations. The derivation of Eq. 5 is given in B.

Equation 4 entails, for each iteration $n$, the inversion of $N_{\Lambda}N_{\Omega}$ matrices, whose size is the number of volumes in the computational domain. To enhance the computational efficiency, three distinct strategies have been put in place. First, by noting that the stiffness matrix does not depend on $n$, we are able to compute LU factorizations for the first iteration, and reuse them until convergence. Another major speed up is achieved by vectorizing the assembly of the stiffness matrices, while also exploiting their sparsity. Lastly, for each $n$, Eq. 4 consists of a set of independent problems, thus it can be trivially parallelized [28]; to this end, we assign a set of phonon directions to each core ($\Omega_{c}$), while the indices $m$ run serially. The implementations of these three methods rely on SCIPY, NUMPY [29] and MPI4PY [30], respectively.

When vectorization and parallelization are both used, we may save memory by noting that, as shown in A, the stiffness matrices have the form $\mathbf{A}_{ml}=\mathbf{I}+\Lambda_{m}\mathbf{A}_{l}$; thus vectorization can be done only for $\mathbf{A}_{l}$, while $\mathbf{A}_{ml}$ is assembled on the fly right before LU factorization. Similar arguments hold for $\mathbf{P}_{ml}$ and $\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}$. The overall methodology is summerized in Algorithm 1, where, for readability, only $\mathbf{A}_{ml}$ has been expanded upon.

The easiest way to run OpenBTE is to use it as a module in Python. A sketch for the workflow is illustrated in Fig. 1. The core module is Solver, which runs both the heat conduction equation as a first guess to Eq. 4 and all subsequent iterations, outlined in Algorithm 1. Before invoking it, though, it is necessary to define the material’s internal boundaries using the Geometry module, and bulk-related data via the Material module. The latter is the interface to external packages. Once the simulation is finalized, results can be plotted using the Plot module. Note that OpenBTE follows a pure functional approach, where data is simply stored as dictionaries. Each module is expanded upon in the next sections.

The Material module converts mode-resolved bulk data, $C_{\mu}$, $\tau_{\mu}$ and $\mathbf{v}_{\mu}$, into the spherically-resolved quantities $\mathbf{G}_{ml}$ and $C_{ml}$. The input files, rta.npz, are created by postprocessing data obtained by external packages. Currently, OpenBTE provides interfaces to AlmbaBTE [9] and Phono3Py [12]. A small set of precomputed rta.npz is also provided. The file created by Material is stored into material.npz and ready to be used by Solver. Depending on the geometry, this module has two models:

• 1.

  rta3D: this model must be used for geometries with $l_{z}>0$. In this case, the angular interpolation of $\mathbf{F}_{\mu}$ is performed in both polar and azimuthal angles.

• 2.

  rta2DSym: when $l_{z}$ is not specified, the actual simulation domain is 2D, and this model must be used. Note that when the base material is 3D, using this option is equivalent to considering an infinite thickness. In fact, in this case, the mode-resolved $\mathbf{F}_{\mu}$ are mapped into $\mathbf{F}_{tk}=\Lambda_{t}\cos(\phi_{k})$, where $\Lambda_{t}$ implicitly includes the MFP and the azimuthal angle, and $\phi_{k}$ is the polar angle on the x-y plane [23].

The model as well as the discretization grid can be chosen with

where, by default, the material file is taken from the database, whose entries are updated online. In this case, we choose to discretize the polar angle into 48 slices.

The Geometry module allows the user to build an arbitrarily patterned nanomaterial. It interfaces with GMSH [31] to create an unstructured Delaney mesh, and creates the file geometry.npz. While we plan to handle user-provided meshes, OpenBTE currently supports 2D membranes with finite and infinite thickness. The size of the unit cell is specified by lx and ly, while the thickness is defined by lz. The direction of the applied gradient is indicated with direction (default is x). Periodic boundary conditions are applied along all directions; however, adiabatic surfaces can be imposed along specific directions with Periodic=[True,False,True], where each component declares whether a surface is periodic. Currently, there are two geometry models, as expanded upon in the next sections.

Simulations results can be plotted via the Plot module. Currently, the following features are implemented:

• 1.

  vtu: the temperature and flux maps are stored in a vtu file and ready to be opened by Paraview [32].

• 2.

  maps: the maps are shown as web-ready, interactive structures.

• 3.

  kappa_mode: this model computes the mode-resolved thermal conductivity from Eq 21, and stores in a file.

• 4.

  line: this option, which currently is implemented only for 2D systems, allows for line plots, using linear interpolation.

In some cases, it is convenient to use Fourier’s law to estimate macroscopic geometric effects. In this first part, we thus aim to assess the reliability of the diffusive solver, implemented in OpenBTE. To this end, we consider a membrane with circular pores and infinite thickness, for which $\kappa_{\mathrm{fourier}}/\kappa_{\mathrm{bulk}}=r$ can be evaluated analytically. In fact, in this case, the thermal transport suppression is given by the Eucken-Maxwell model $r=(1-\phi)/(1+\phi)$ [33], $\phi$ being the porosity. In Fig. LABEL:fig1a, we compare this model with the results from OpenBTE for different porosities, obtaining excellent agreement ($<1\%$). We note that by using the option only_fourier=True it is possible to run only the Fourier’s model. Next, we compute $r$ for a more complicated structure, comprising two concatenated sub-lattices of pores with different user-provided shapes; this structure is generated using the shape=custom option from the lattice model, a two-pore base, and different options for each pore. As shown in Fig. LABEL:fig1a, we obtain a much larger heat transport suppression due to the longest path heat must travel through the material.

Let us now compute $\kappa^{\mathrm{eff}}$ using the BTE, encoded in Eq. 4. For the chosen porosity of $0.1$ and $L$ = 50nm, we obtain  17.4 Wm^-1K^-1, well below the corresponding macroscopic value 99.4 Wm^-1K^-1. In Fig. LABEL:fig1b, we plot a cut of the temperature maps, both for the BTE and Fourier cases, along the $x$-axis. As expected, they have different trends, the BTE one being flatter due to ballistic transport. When a finite thickness is added, phonon undergo stronger suppression; in fact, for thickness $t$=10 nm and porosity 0.1, we obtain $\kappa^{\mathrm{eff}}=8.9$ W m^-1K^-1, roughly a half of the value with $t=\infty$. In Fig. LABEL:fig1c, $\kappa^{\mathrm{eff}}$ for different porosities is plotted. From the insets, we note that, as the porosity becomes larger, the heat flux peaks around the phonon bottleneck, e.g. the areas between the pores. Lastly, the mode-resolved thermal conductivities for $t$ = 10 nm and $\infty$ are plotted in Fig. LABEL:fig1d; analogously to Ref. [23], we can appreciate the stronger suppression for large MFPs for the case with finite thickness.

We have presented OpenBTE, a software for computing deterministically space-dependent heat transport in 2D and 3D structures. The efficient and modern implementation of the underlying algorithms enable accurate simulations both in the cloud and on common laptops, thus democratizing this type of simulations. Future work includes going beyond the RTA of the BTE and adding transient dynamics. GPU support and differentiability are also on the roadmap. Furthermore, we plan to develop a hybrid Fourier/BTE solver, which has the potential to dramatically reduce runtimes for the cases with 3D structures. Filling an important gap in the current offer of publicly available phonon transport solvers, OpenBTE sets out to guide experiments on nanoscale heat transport and facilitate high-throughput prediction of thermal properties of nanostructures.

Research was partially supported by the Solid-State SolarThermal Energy Conversion Center (S3TEC), an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award No. DESC0001. The author thanks all the early users, who helped improving OpenBTE.