{'fraction_non_alphanumeric': 0.037213032581453635, 'fraction_numerical': 0.009363408521303258, 'mean_word_length': 4.57709940735771, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
QEHeat: An open-source energy flux calculator for the computation of heat-transport coefficients from first principles

We give a detailed presentation of the theory and numerical implementation of an expression for the adiabatic energy flux in extended systems, derived from density-functional theory. This expression can be used to estimate the heat conductivity from equilibrium ab initio molecular dynamics, using the Green-Kubo linear response theory of transport coefficients. Our expression is implemented in an open-source component of the QE suite of computer codes for quantum mechanical materials modelling, which is being made publicly available.

first-principles methods, a definition of the heat-flux compatible with density-functional theory is required [1]. The evaluation of such a heat flux requires an extension of stateof-the-art atomic simulation codes. Solution method: This work describes in detail the numerical implementation of the adiabatic energy current derived in Refs. [1,2] and makes it available to the users of the Quantum ESPRESSO suite of computer codes [3]. Used in conjunction with the cp.x code, to perform Car-Parrinello ab initio molecular dynamics, and the SPORTRAN post-processing tool for data analysis [4][5][6], the program allows to estimate heat transport coefficients in extended systems entirely from first principles. The new code provides as well to developers a modular and easily extendable framework to evaluate time derivatives of electronic properties (e.g. electronic densities or potentials) via a finite difference approach.

Introduction
The Green-Kubo (GK) linear-response theory of heat transport [1][2][3][4] has long been deemed incompatible with modern quantum-simulation methods based on Density-Functional Theory (DFT) [5,6], essentially because the concepts of energy density and current, upon which that theory stands, are ill-defined at the molecular scale. This predicament has been reversed by a recent paper [7] by Marcolongo, Umari, and Baroni (MUB) where it was shown that, due to a general gauge invariance principle of transport coefficients [8], the heat conductivity is largely independent of the details of the microscopic definition of the energy densities and currents from which it is derived. Leveraging this remarkable finding, an explicit expression for the energy flux, based on DFT, was given in that same paper. This expression was implemented in a private branch of the Quantum ESPRESSO project [9][10][11] and, due to its considerable complexity, it has not been publicly available so far to the scientific community at large. The purpose of the present paper is to provide a detailed derivation of the MUB expression for the DFT energy flux, a description of its implementation in Quantum ESPRESSO, and to release and document an open-source code for its computation, named QEHeat. QEHeat can be easily interfaced to read a dynamical trajectory generated with a code of choice and compute the MUB flux for the corresponding steps. QEHeat is already delivered with a user-friendly interface for the cp.x program of the Quantum ESPRESSO distribution, which is routinely used to perform Car-Parrinello molecular dynamics simulations [9][10][11]. The combination of QEHeat, cp.x and the post-processing tool SPORTRAN [12], designed to perform the statistical analysis needed to evaluate transport coefficients, provides a convenient framework to compute the heat conductivity of extended insulating systems-be they crystalline, amorphous, or liquid-entirely from first principles. We warn the user that the MUB flux is an adiabatic energy flux. This refers to the fact that electrons are supposed to populate the ground state during time evolution, which is the case for insulators with a finite band gap.
We note that other computational efforts have been proposed to evaluate microscopic expressions for the energy flux based on DFT. For example, Kang et. al [13] evaluate an atomic decomposition of the total energy and then use a two-step procedure to take into account periodic boundary conditions when evaluating the energy flux. Carbogno et. al [14] neglect the convective contribution to the energy flux and use a DFT-based expression for the virial component. This approach is developed as an approximation suited to the description of thermal transport in solids. The open-source distribution of QEHeat will make it easier to compare the computational advantages of the various expressions which, after the assessment of the principle of gauge invariance, are being developed by the community. Finally, the evaluation of the energy flux provided by QEHeat involves a finite difference evaluation of electronic properties, e.g. electronic densities or potentials, which are performed entirely in-memory, storing the results of different total energy calculations into ad-hoc data structures. Therefore, QEHeat provides as well to developers a modular and easily extendable framework to evaluate time derivatives of electronic properties via a finite difference approach.

Overview
According to the GK theory of linear response [1][2][3][4], the heat conductivity of an isotropic system is given by: where Ω is the system volume, J (t) = Ω j(r, t)dr, i.e. the volume integral of the energy current density, is the energy flux (note that with this definition the flux is extensive), T the temperature, and · indicates an equilibrium average over the initial conditions of a molecular trajectory. QEHeat serves the purpose of evaluating J at the DFT level of theory. A formal expression for the energy flux can be obtained by integrating by parts the continuity equation˙ = −∇·j as [5,6]: where (r) is the energy density of the system, whose integral is its total energy. Strictly speaking, Eq. (2) is ill-defined in periodic boundary conditions (PBC), which are commonly adopted in molecular simulations. In order to compute it explicitly within PBC, one has to first recast it in a boundary-insensitive form, and Ω can then be replaced with the volume of the simulation cell in Eqs. (1)(2). Once this is done, the DFT energy flux can be cast into the MUB form, called here J MU B , and discussed in detail in later sections.
QEHeat computes the MUB energy flux as a function of the atomic positions, {R s } and velocities, {V s }, i.e. for any selected snapshot of an ab initio molecular dynamics (AIMD) trajectory. Despite the complexity of the resulting formula for the energy current, from a practical point of view the use of QEHeat relies on a limited number of additional input parameters with respect to a standard Quantum ESPRESSO DFT computation. These are reported in the energy_current input namelist, which is shown in figure 1. The meaning of all the keywords is explained in more detail in section Section 5. The only additional parameters are eta and n_max, controlling the Ewald summations, which appear only in classical contributions to the energy current, and delta_t, a time-discretization parameter, used to perform numerical derivatives. The default values should work for most systems. &energy_current delta_t = 1.000, file_output = 'current_hz', eta = 0.100, n_max = 5, trajdir = 'traj/cp', first_step = 1, vel_input_units = 'CP' / Figure 1: Example of the energy current namelist. delta t is the time used for numerical derivatives. eta and n max are the parameters used to converge Ewald sums. trajdir is the prefix of the trajectory files. In this example, the program reads the files "traj/cp.pos" and "traj/cp.vel". first step tell the program the first step id to compute. The step ids are part of the trajectory file format. After this namelist the full pw.x input is required. A full example of the input and the documentation of the keywords can be found at https://gitlab.com/QEF/q-e/-/tree/master/QEHeat/Doc. See also Section 5.
The work is organized as follows. In Section 3 we show the formulas that QEHeat implements for the abiabatic energy flux, and Appendix 3 is dedicated to the exact numerical schemes. Section 4 contains an overview of the code, and Section 5 describes the input parameters in more details and shows some stability and implementation checks. Finally, in Section 6 we showcase how to use the energy flux time series provided QEHeat to evaluate thermal conductivity, using the SporTran data analysis post-processing tool (the data analysis is fully explained in the Jupyter notebook data analysis.ipynb in the Supplementary Materials).

DFT energy flux
In this section we recall the expression of MUB energy flux [7] and some of the notation used throughout the paper. For a more extensive and detailed study of the implementations of the many components of the MUB current the reader is referred to Appendix B. J MU B is expressed as a sum of five components: where In the following, these components are referred to as the Kohn-Sham, Zero, Ionic, Hartree and Exchange-Correlation fluxes, respectively. For the insulating systems of interest in this work, the ionic degrees of freedom completely define the state of the system and the electrons populate the ground state, according to the adiabatic approximation. Each time derivative, indicated with the usual dot operator, has then to be understood from the implicit dependence on the atomic positions. We note that a complete understanding of the different components of the MUB energy flux is not needed to perform, as a user, a thermal-conductivity calculation. Here and in the following, we indicate with L the lattice vector. We stress that periodic boundary conditions (PBC) are assumed here over the unit (simulation) cell. This implies that Kohn-Sham orbitals and energies are sampled at the Γ point of the Brillouin cell. The notation ∇ sL is a shorthand for the gradient with respect to displacement of the atom at location R s + L. A summation over s runs over all atoms belonging to the simulation cell. Unless otherwise specified, carets indicate quantum-mechanical operators, as inĤ KS orr. Following is a brief report of the definition of the most important terms according to their physical meaning. For a more extensive summary of the notation used here and throughout the text the reader is referred to Appendix A ("Notation"). The ionic energy e 0 s is the sum of the kinetic energy, 1 2 M s V 2 s , and w s the classical electrostatic interaction between the s-atom and all other atoms in the system. The electronic degrees of freedom are instead described by the instantaneous Kohn-Sham Hamiltonian,Ĥ KS , and its eigenvalues and eigenvectors, the Kohn-Sham energies, , and orbitals, ϕ . A summation over runs over all occupied orbitals. Several quantities, like the electronic number-density distribution, n(r), and Hartree potential, H (r), are implicit functions of the wave-functions.ˆ 0 is the total external atomic pseudo-potential, describing the interaction between electrons and nuclei. For a more detailed analysis see Appendix B.2. The symbol GGA stands for the generalized gradient approximation (GGA) exchange-correlation local energy per particle and its derivative with respect to density gradients is indicated with ∂ GGA , which is a vector whose component along direction i ∈ {x, y, z} is given by ∂ GGA (n, ∇n)/∂(∇ i n). In the present version of QEHeat only the local density approximation (LDA) and the generalized gradient approximation given by the PBE [15] functional are implemented. The formulas reported are compatible with PBC. Thus, they can be implemented for periodic systems, where some of the summations need to be extended to all the periodic replicas of the atoms. Nevertheless, the computation of the various contributions to the energy current, Eqs. (4)(5)(6)(7)(8), is plagued by the occurrence of several divergences, arising from the long range character of the Coulomb interaction. As it is the case for the total energies, atomic forces, and stress, the individual electronic, ionic, and electron-ion contributions diverge and it is only their sum that is regular in the thermodynamic limit. In order to regularize the individual components of the MUB flux, we compute all the relevant terms by screening the Coulomb interaction with a Yukawa cutoff, 1 x → e −µx x . In Appendix B, we check explicitly that the singular contributions to the various terms cancel each other in the µ → 0 limit, so that they can be consistently and safely neglected and do not appear in the final formulas reported in Appendix B.

Code structure
We start by describing the strategy implemented to compute numerical derivatives of quantities appearing in Eqs. (4), (7), and (8), likeṅ and˙ H , since they require a special treatment. Quite generally, one needs to evaluate terms of the typeḟ ({R s (t)}), where the function f can be a scalar function, which depends on time only through the set of the instantaneous ionic positions {R s (t)}, evolving according to Hamilton's equations of motion. QEHeat implements a finite-difference scheme, using by default a symmetric numerical differentiation formula: The small parameter dt is an input of the computation. In such a scheme quantities that are not differentiated are evaluated at time t, so three wave-functions are required to be kept in memory at the same time. QEHeat performs therefore for each step two additional self-consistent-field (SCF) DFT calculations, using the same DFT solver of the Quantum ESPRESSO distribution, at slightly displaced positions, i.e. {R s − V s dt/2} and {R s + V s dt/2} along the AIMD trajectory. The wave-functions of the previous calculation are used as a starting point for the next one, which require much less iterations to converge. We note that QEHeat gives the user also the possibility to use a non-symmetric differentiation scheme, which is shown in Appendix E. This scheme is computationally cheaper. Nevertheless, for differentiable functions, the order of convergence of the symmetric scheme is quadratic in dt, whereas the non-symmetric one is linear. Accordingly, the stability is improved with the default symmetric scheme. We recommend therefore the latter and use it for all calculations here presented. See also the dedicated Appendix E.
The trajectory data is managed by the Fortran derived data type cpv_trajectory defined in the file cpv_traj.f90, while the orbitals and the associated atomic position are managed by the derived type scf_result, implemented in scf_result.f90. The most relevant subroutine that acts on this object is scf_result_set_from_global_variable, that copies the eigenfunctions, the eigenvalues, the potential and the atomic positions from the Quantum ESPRESSO's global variables to the instance of scf_result. The results for each of the three (or two) wave-functions that are required by the computation routines are stored in the variable scf_all, defined in the main program routine. Global variables are avoided as mush as possible.
The code starts by reading the input "namelists": first the energy_current namelist, then all the pw.x namelists. Then it calls all the pw.x-related initialization routines. After eventually reading the previously generated output file that allows the program to set the correct starting timestep, it enters the main loop over the input trajectory timesteps. The trajectory files have the same format of Quantum ESPRESSO's cp.x code output files.
The most important routines where the above mentioned data structures are used are the following: Carries out the computation of Eq. (5). This routine is called in the middle of the computation using the same timestep t of the positions stored in the input trajectory, so that the result does not depend on dt.
• SUBROUTINE current_ionic (module ionic_mod) Computes all parts of Eq. (6), and it is called as current_zero at the same timestep of the input trajectory • SUBROUTINE current_hartree_xc (module hartree_xc_mod) Computes (7) and (8). Since a numerical derivative is needed, this routine reads the wave-functions from the global type scf_all and it is run at the end of all necessary run_pwscf calls.
• SUBROUTINE current_kohn_sham (module kohn_sham_mod) Computes (4). As current_hartree_xc, it needs all the wave-functions calculated by the DFT solver for this step.
• SUBROUTINE run_pwscf Uses Quantum ESPRESSO's routines to solve the DFT problem for the atomic positions stored in the global array tau. Equivalent (but the starting wave-function and potential, that can be the last computed one) to a standard call to the pw.x program with the input stripped of the ENERGY_CURRENT namelist. The result is stored in the Quantum Espresso's global arrays (evc) • SUBROUTINE prepare_next_step This routine is used to change the global array tau to tau + vel·dt·ipm, where ipm is the argument of the subroutine that can be -1,0,1. After doing that it calls the necessary routines to prepare the potential for run_pwscf.
The 4 modules, one for each part of the MUB current, are completely independent of one another. The structure of the main loop over the trajectory's time steps is summarized in Algorithm 1. As every big computational code an extended test suite is needed to safeguard the correctness of the calculation after every source code modification. We implemented small tests that are able to run on a single core of a cheap laptop that check against changes in the numerical output of many parts of the code, using the standard Quantum ESPRESSO's test suite framework.
To conclude the section we want to do some remarks on the code and its interactivity with others typical ab initio simulations tools. In principle the wavefunctions computed on-the fly by cp.x during the AIMD run could be used, but we preferred to implement a workflow where the computation of the currents is completely decoupled from the AIMD engine, thus the wave-functions are always recomputed by pw.x. The chosen approach allows the user to run the calculation in post-processing mode, thus using the preferred code to generate the dynamics, not to be limited to those in the QE packages. This way it allows, also, a trivial and powerful per-snaphot parallelization. call current_zero, evaluate currents derived from the pseudo-potential 7: call current_ionic, evaluate the electrostatic and kinetic Ionic current 8: call run_pwscf with positions displaced at t + dt/2 9: call current_hartree_xc, evaluate Excange and Hartree currents 10: call current_kohn_sham, evaluate Kohn-Sham current Steps 6 and 7 do not require any finite differences, while steps 9 and 10 do.
Step 10 is the most expensive.

Input description
The input is organized in a traditional fortran namelist input file, similar to the input files of many Quantum Espresso's programs, and an optional trajectory file (that is a file for the atomic velocities and a file for atomic positions) if the user wants to compute the energy current for more than one snapshot with a single run. A full example of the input can be found at https://gitlab.com/QEF/q-e/ -/tree/master/QEHeat/examples. Before running QEHeat it is necessary to obtain velocities and positions from a different code, for a complete description of the units of measure see Table 1. If the Quantum Espresso's cp.x program is used for this purpose, its output trajectory files can be recycled as input trajectory files without any modification. The program's mandatory input is organized into an ENERGY_CURRENT namelist and all the usual pw.x namelists. We remind the user that, up to the present version, only norm conserving pseupotentials and the PBE exchange correlation functional are supported. At the end of the input file the ATOMIC_VELOCITIES card is required. In the IONS namelists the value ion_velocities = 'from_input' is required, since the program must read the atomic velocities to compute the energy current. An extensive input description can be found in the file INPUT_ALL_CURRENTS.html in the Doc folder of the code repository. Here we remark the most important parameters of the ENERGY_CURRENT namelist: • delta_t : time in PW's atomic unit used to compute all the numerical derivatives like the one in Eq. (9); • trajdir : prefix of the cp-formatted trajectory. Optional: if not setted, only the positions and the velocities of the input file are read; • n_max : number of periodic images along the directions of each basis cell vector to converge Ewald sums. This fixes the range of L in Eq. (6) ; • eta : convergence parameter of the Ewald sums needed in the computation of J n , for more details see Appendix B. 3.
An example of the namelist is provided in fig. 1. An additional output file is written and updated at the end of each step in the folder where the program is run. All the currents are printed in a column format, ready to be analyzed by an external post-processing tool.
As discussed in Appendix E, CONTROL's conv_thr and ENERGY_CURRENT's delta_t have a profound link and influence heavily each other, and despite we think the default value of delta_t=1.0 is safe enough, they must be carefully tested, veryfing that the standard deviation of the result is low enough.
The standard deviation of the output energy current can be estimated by repeating the same calculation for every step, many different times, setting for each repetition a random starting potential and a random starting wave-function. The input options re_init_wfc_1 = .true. together with n_repeat_every_step = 20, for example, do 20 repetition of every timestep, resetting the starting wave-functions/potential before the first scf calculation. The pw.x's input option startingwfc = 'random' is suggested, to obtain a faithful error estimation. If more reinitializations are desired, the options re_init_wfc_2 and re_init_wfc_3 can control the randomness of the starting wfc and potential of every of the 3 (or 2) wave-functions needed to perform the numerical derivatives, as explained in section 4. Note that when the wave-function is reinitialized from scratch, the computation time raises since more scf cycles are required to reach the target convergence threshold. When n_repeat_every_step is greater than 1, an additional column formatted output file with the averages and the standard deviations is produced.

Implementation checks: Finite systems translating at constant speed
The Green-Kubo current associated with a localized energy density (r, t) rigidly translating with constant velocity , is equal to E tot × . One possible way to show this is to consider (r, t) = (r − t, 0) ≡ 0 (r − t),thus J a = ˙ r a dr = Note that we used the fact that 0 can be taken identically equal to zero at the boundary of the integration volume, to remove boundary contributions from the integraton by parts. The identity requires therefore the energy density to be localized and this condition can be mimicked in PBC considering a large enough supercell. We used this property to check the correctness of our implementation for each individual current in Eq. (3). We simulate a single Argon atom and a water molecule at equilibrium, both translating at constant speed. We then compare J MU B output from QEHeat and E tot , where E tot is evaluated using an independent computation from the QE code. As discussed, the resulting currents need to be equal only in the limit of large cells, where boundary effects can be neglected, i.e. the energy density is truly localized, and under tight convergence criteria. In Fig. 2, we report the ratio between the computed and theoretical values as a function of the cell parameter, showing that the correct limit behavior is recovered. For this calculation we used a cutoff of 120Ry and econv = 10 −14 Ry. Additionally, in Appendix C we perform the same test using a large cell parameter but removing individual current components (i.e. the {XC,IONIC,ZERO,KOHN} components) from the total energy flux. Since removing each component changes the difference with respect to the limiting theoretical value, this proves the correct implementation of each individual current.
In the same figure, using the same approach, we tested the electronic density current defined in eq. (B.58) as well. In the infinite cell limit the electron density current of a system translating at a constant speed is J el = N e where N e is   Proof that the behaviour of the code is correct for systems translating at equilibrium and constant velocity. In this setting, the output of QEHeat can be exactly compared with the known theoretical value, indicated with J theory , in the large cell limit. Indeed the plot shows that the ratio J/J theory goes to one increasing the simulation cell. This test can be performed for J equal to the electronic density current (dotted), that is used to calculate a part of the energy current, and for the energy current itself (not dotted). Tests are performed for a single relaxed water molecule (H 2 O) and a single argon atom (Ar). The electronic current should be, in the infinite cell limit, J el theory = N el × , where N el is the number of electrons and the energy flux should be, in the infinite cell limit, equal to J ene theory = E tot × , E tot being the total energy. the number of electrons. It is possible to see in figure 2 that the correct limit is obtained, validating the implementation of the code.

Example of a thermal conductivity calculation
As an example of a complete QEHeat calculation, coupled with the SPORTRAN tool for signal analysis, we report the results of a thermal conductivity calculation for PBE ( [15] ) water. It is known in the literature [16][17][18] that DFT with PBE functional does not provide a correctly structured water. Nevertheless, PBE water is a simple example of a single component molecular fluid that can be used to show all aspects of the theory. For a complete description of the signal-analysis theory used here to post-process the simulation, based on cepstral analysis, we refer to the works [5,6,8,19] and to [12] for example notebooks and tutorials. The details of the present analysis are reported for reproducibility in the Jupyter notebook "data analysis.ipynb", provided as supplementary material, and we report here the results. We warn the reader that the heat flux returned by QEHeat almost always contains non-diffusive signals, i.e. signals not contributing to thermal conductivity, with a high amplitude, that may be a problem for a standard numerical Green-Kubo integration. The high formation energies present in abinitio computations can generate such components in the energy flux. By using the fact that these components are highly correlated with the mass and electronic density flux, renormalization or multicomponent techniques, effectively decorrelating the energy flux from these spurious signals, solve this issue. For a longer discussion we refer to [19] and references therein. Here we took advantage of the multicomponent filters implemented in SPORTRAN. We consider a 64 molecules water system at 600K. First, the molecular dynamics trajectory was computed with the cp.x code using a cubic box with size of 12.43Å and an integration timestep ∆t = 3 in Hartree atomic units, with a total length of 240ps. QEHeat shows excellence stability properties as a function of the time-discretization parameter dt, which are shown and thoroughly discussed in Appendix E. A value of dt equal to twice the Car-Parinello MD simulation timestep, dt = 2∆t, is of particular interest because it would be beneficial in an on-the-fly computation, allowing to reuse the same wavefunctions computed in the MD simulation, neglecting the need for the recomputation of the scf cycles.
In the Appendix we also show that, at a fixed value of the self consistency threshold for convergence econv, the variance of the MUB estimator decreases with dt. Thus, it is favourable, in order to control numerical noise, to choose dt as large as possible, as long as the bias introduced by nonlinear effects is negligi-ble. For these calculations, we first investigated two values of econv, namely 10 −8 and 10 −11 Ry with a rather small dt = 0.5 Rydberg atomic units. The final result of the computation, i.e. the thermal conductivity coefficient, is affected very weakly by this parameter, that on the contrary changes in a noticeable way the computational cost. In this case both computations gave (1.07 ± 0.09)W/mK. Moreover we checked that using a dt that is equal to the integration timestep of the Car-Parrinello simulation and econv = 10 −8 produces the same spectrum with no noticeable differences and the same thermal conductivity coefficient: in this case we had (0.98 ± 0.09)W/mK. According to our experience, the overall cost of a QEHeat calculation is often of the same order of magnitude of the full Car Parrinello molecular dynamics simulation even if the QEHeat computation can be trivially parallelized. The data analysis with SPORTRAN adds a negligible computational cost. A more detailed study of the computational costs and how to reduce it can be found in Appendix F. Here we just stress that the total computational time can be also reduced choosing an appropriate length of the simulation, according to the desired precision. Fig. 3 show the dependence of κ from the simulation length both for our water system (upper panel) and a solid ice X structure from ref. [20] (lower panel), where the oxygen are packed in a bbc lattice and tetrahedrally coordinated to hydrogen atoms located exactly midway between two neighboring oxygen atoms. In the present didactic work we choose a very long simulation time of 240 ps but from the figure is it clear that a simulation length around 100 ps can provide a reasonable estimate. Depending on the system, even shorter simulations lengths could suffice: the lower panel of Fig. 3 presents the same data of ice X from [20], where a simulation of 20 − 30 ps proved to be sufficient. We also warn the reader that for strongly harmonic systems, e.g. solids at ambient temperature and pressure, also large cell sizes may be required to remove boundary effects.

Conclusions
We have presented a robust implementation of the MUB energy current, allowing to compute the heat conductivity from an explicit expression of the energy flux based on DFT, something that is made possible only by recent findings The dependence of κ with respect to the trajectory length for 64 molecules of water system at 600K and for the ice X system in [20]. The orange line correspond to the values with the longest trajectory possible, κ = 1.07 ± 0.09 W/(mK) and κ = 16.1 ± 1.1 W/(mK) for water and ice X respectively. The data for the ice X system are taken from the Material Cloud repository [21] matic framework to evaluate time derivatives of electronic properties, e.g. electronic densities or potentials.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
[21] F. Grasselli, L. Stixrude, S. Baroni Notations and definitions used throughout the text : • e : electron charge ; • e 0 s : ionic energy, equal to 1 2 M s 2 s + w s , where w s is the electrostatic energy.
• XC : local XC energy per particle, defined by the relation: E XC = XC [n](r)n(r)dr. "LDA" and "GGA" in Eq. (8) indicate the local-density and generalizedgradient approximations to the XC energy functional ; • : electronic eigenvalues ; •f (G) : Fourier transform of periodic functions. Given a function f (r) periodic with respect to the unit cell, the Fourier transformf (G), evaluated at the reciprocal lattice vector G is defined by the convention: • F [ f ](G) : Fourier transform of non periodic and localized functions. Given a function f (r) smooth enough and localized (hence not periodic), its Fourier transform is defined for every value of G (hence not only for reciprocal lattice vectors), by the convention: where a convenient normalization factor Ω has been introduced ; •Ĥ KS : instantaneous Kohn-Sham (KS) Hamiltonian ; • L : lattice vector, as defined by the unit cell ; • M s : atomic mass of atomic s ; • n(r) : ground-state electron-density distribution, defined as n(r) = |φ | 2 ; • H : Hartree potential, defined as: Note that if n(r) is periodic, also the Hartree potential is periodic; •r : multiplicative position operator or versor along direction pointed by vector r. The meaning should be clear from the context. • φ, |φ : a generic normalized wave-function. With brackets, the same wavefunction is considered as a vector of a Hilbert space with a scalar product ; • Ω : volume of the unit cell ; • ∂ GGA : derivative of the GGA XC local energy per particle with respect to density gradients. It is a vector whose component along direction i is explicitly given by ∂ GGA /∂(∇n i ) , • ∇ : gradient with respect to the spatial coordinate r ; • : standard scalar product between wave-functions ; •[ ] : derivative with respect to time ;

Appendix B. Numerical implementation of each current component
In the following appendix we will give an exhaustive and detailed description of the implemented components of the MUB current.

Appendix B.1. Kohn-Sham current
Starting from Eq. (4) and after some simple algebra [22] we end up with the following expression for the Kohn-Sham current. For every Cartesian component, i, one has: andP v andP c are the projectors over the occupied-(valence-) and empty-(conduction-) state manifolds respectively.φ c is calculated by solving the linear system: where [·, ·] indicates the commutator between quantum mechanical operators and α is a positive constant that removes the singularity of the linear system and forces the solution to be orthogonal to the valence manyfold, using standard techniques from density-functional perturbation theory (DFPT) [23]. In order to avoid alignment problems between wave-functions at different time steps,φ c is calculated by moving the derivative to the projector [10,22], using the relation: which can be evaluated using a finite-difference scheme, as explained in detail in Sec. 4, reading: where dt is a time-discretization parameter, input of a QEHeatcomputation, which is discussed in more detail in Section 4 and Appendix E. After evaluating Eq. (B.8) and solving Eq. (B.6), the results can be inserted into (B.1). We note that the computational cost to evaluate all the components of the MUB current is dominated by the solution of the linear system, Eq. (B.6).

Appendix B.2. Zero current
The Zero current J 0 , Eq. (5), can be better treated by separating the local and non-local contributions from the ionic pseudo-potential,ˆ 0 . Each contribution acts on a generic wave-function φ in the following way: where the total pseudopotential has been separated into atomic contributions and into its local (LOC), long-tailed, and non-local (NL), short-range, components. The subscripts in the expressionsˆ LOC sL andˆ NL sL indicate that the corresponding ionic pseudo-potential is centered at the atomic position R s + L. Note that the atomic contributionsˆ LOC/NL sL , when differentiated with respect to ionic positions, only depend on the position of the atom located at R s + L.
The local and nonlocal potentials provide two contributions to the Zero current, which we discuss individually in the next sections. For the local contribution, f LOC s (r) is the radial local pseudo-potential provided in the pseudo-potential datasets for each species. The centered beta functions β s lm (r) (also called projectors) define the non-local component of the pseudo-potentials and are defined as β s lm (r) ≡ β s l (r)Y lm (r), where Y lm are the real spherical harmonics with quantum numbers l, m andr is the unit versor, not to be confused with the multiplicative position operator. We denote the radial components of the beta functions with the similar notation β s l (r). These are the ones provided for each atomic species in the pseudo-potential datasets, alongside the constant D matrix. At every time step, the beta functions need to be centered on the instantaneous ionic positions and in Eq. (B.11) we indicated the translated beta function for atom at position R s − L with an apex. More explicitly : We fix some handy notation and define from a (real) localized function γ(r) (e.g. a β function or the local pseudo-potential) its periodic counterpart as: whose Fourier components can be computed as: where we introduced the symbol F [γ](G) to identify a standard Fourier transform over all the three dimensional space, here defined for localized functions. In the following we use the notation γ a to indicate the translated and periodic function built from its localized counterpart, analogously to the notation introduced for the projector. Note that the following relation is used in the code, γ a (G) = e −iG·a F [γ](G), thanks to standard properties of the Fourier transform.
The remaining Fourier transforms are evaluated numerically by QEHeat whenever needed, as explained below.

Appendix B.2.1. Zero current: local contribution
The local pseudo-potential of the s-atom f LOC s (r) behaves as ∼ −Z s e 2 /r for large r. We call this local long-range contribution to the flux J LR 0 : We use the chain rule and the definitions: where f LOC s (r) is the derivative of the local pseudo-potential. Note that both h and u are periodic functions. The current can then be rewritten as: The reciprocal Fourier components of h s i j (r) are computed through the following procedure, which avoids an explicit numerical differentiation of the pseudopotential. One writes: The resulting expressions ofh 1 (G) andh 2 (G) can be evaluted using Eq. (B.14) and the standard expansion of e −iG·r into spherical harmonics and Bessel functions where the G = 0 component vanishes thanks to the presence of a derivative in the definition of h 1 . For h 2 the situation is different: These expressions need still to be a bit elaborated before being evaluated by QEHeat. In order to evaluate integrals of the localized functions, one needs to add and subtract the asymptotic long-range tail of the local pseudo-potential. The long-range part can be integrated analytically for any finite value of the Yukawa screening parameter, µ, after plugging in the exact form of the spherical Bessel functions. For h2, this permits to extract the divergent part in the G = 0 contribution as well. The final results read: (B.24) We note here again that only from h2 we get a Coulombian divergence when µ → 0. Overall, the divergent part of the Zero current is equal to where Z tot = s Z s .

Appendix B.2.2. Zero current: nonlocal contribution
The non-local part is inherently short range and we call it J S R 0 : For a pair of two localized functions (γ 1 (r), γ 2 (r)) we introduce the notation, given a generic translation a: where we recall that γ a 1 and γ a 2 are the translated and periodic counterpart of localized functions. We note that, once the Fourier components of φ, γ 1 and γ 2 are known, evaluating A[γ 1 , γ 2 ](a) involves just scalar products between periodic functions and is straightforward. Using these expressions and expanding the projector operator in their integral form, the current can be rewritten as: (B.28) In order to evaluate Eq. (B.28) via (B.27) the Fourier transforms of the following four localized functions are needed : β lm (r), where we dropped the atomic index s. We need only to evaluate the first two expressions thanks to the identities: which are valid for every localized function γ(r). We show how to evaluate the Fourier transform of the second function, r i β lm (r), for i = x, which is the most complex one. The procedure is similar for the remaining expressions. One starts by replacing the factor x using that Y 11 (r) = − 3 4π x r (according to the convention followed by Quantum ESPRESSO for the sign of the spherical harmonics). By expanding e −iG·r as well into spherical harmonics, one gets: where dr indicates an integral over the solid angle, such that dr = r 2 drdr. In Quantum ESPRESSO the Clebsch-Gordan coefficients are stored in a array ap defined as: Application of this relation to the product Y lm Y 11 allows us to simplify Eq. (B.31) to the final form: M, l, m, 1, 1) .
(B.33) The remaining radial integral is performed numerically on a grid.

Appendix B.3. Ionic current
The so-called Ionic flux, J n , is the contribution to the total energy flux that depends only on the ionic positions, R s , and velocities, V s . First of all let's take, as reported in Appendix A, the definition of the ionic energy e 0 s = 1 2 M s V 2 s + w s and the electrostatic energy: where L = |L|. We can separate the expression of the Ionic flux, Eq. (6), in two terms: one depending only on the mass and velocity of the ions, an other depending on w s and its gradient.
where ∇ tL is the shorthand notation for the gradient with respect to displacement of the atom with position R t + L, introduced in Section 2. We, also, introduced f (x) = 1 x to keep track of the Coulombian contributions when applying the Yukawa screening.
For the sake of simplicity, let's introduce the following four quantities: where i, j ∈ {x, y, z} represent the Cartesian coordinates. Then, consider the following properties for the function f (x): The following relations be- Introducing the matrix T i j = S B δ i j + S A i j , which only depends on the cell and not on the individual atoms within it, we can rewrite Eq. (B.45) as a sum of terms that depends only on the atomic species and the cell: j∈{x,y,z} where S is a sum over the atomic species. Recent theoretical developments in the statistical analysis of fluxes [5,7,8,19,24], and in particular the so called gauge [7] and convective [24] invariances tell us that fluxes written as sum of terms that depends only on the atomic species do not contribute to the value of the thermal transport coefficient. Thus, J nB can be neglected in the computation of J n , reducing the implemented formula to: We still need to deal with the long-range Coulomb contributions in J nC , J nD and J nB . Following the scheme explained in Section 2 we introduce a Yukawa screened interaction, 1 x → e −µx x , and we will study the limit µ → 0. We can, then, straightforwardly apply the Ewald techniques [25], with a convergence parameter η, for the computation of S A i j , S B , S C (R s − R t ) and S D i j (R s − R t ). Further details of the computation can be found in Appendix C: tive of the electron number density, to obtain for every Cartesian component i: All quantities needed to evaluate the Electronic density current have already been discussed in the section dedicated to the Kohn-Sham current. The electronic flux thus evaluated is interesting per se, e.g. to compute the electric conductivity in ionic conductors, and also as an ingredient to facilitate the statistical analysis of the energy-flux time series, using multi-component [24] or decorrelation [19] techniques.

Appendix C. Reciprocal space computation of the S A i j and S B
In the following appendix we show explicitly the computation of S A i j , S B . S C and S D i j can be computed applying more straightforwardly standard Ewald techniques [25], then the full derivation is left to reader.
Recalling the definition of S A i j from Eq. (B.41): The term containing erfc( √ η|x − L|) can be simply computed in real space: L , 3) The first term of Eq. (C.3) can be easily computed in reciprocal space, while it can be demonstrated that lim Summing all these pieces together we get the expression in Eq. (B.50): S B is, instead, defined as: introducing a functional dependence on x: (C.6) The second expression can be computed in direct space as it is, while the first requires some further work: (C.7) Now the first term of Eq. (C.7) can be computed in reciprocal space, thus, summing all the contributions, we can regain the expression in Eq. (B.51) for S B : Finally we remark that the two terms computed in this Appendix are general and work for generic cells. However, since J nB does not contribute to the value of thermal transport coefficient, up to the present version of QEHeat we only implemented the simpler expression for cubic cells: where we defined: Tr [·] indicating the trace of a matrix. The following identity can be used to recover the formula for cubic systems from the general one:

Appendix D. Implementation check of individual currents
As discussed in the text, in the case of a finite system at equilibrium translating at constant speed , the current J MU B must be equal to E tot × . In Fig. D.4 we report the modulus of J MU B − E tot , normalized by the module of E tot and indicated with REF ERROR. In order to check the correct implementation of each individual current, we also report values of ∆ X , which represent the same quantity after the substitution J MU B → J MU B − J X , X ∈ {XC, IONIC, ZERO, KOHN}. If the current X is correctly computed, the error with respect to the reference value E tot should increase, as is indeed observed. For Argon, we used the QE parameters ecut = 160 Ry and econv = 10 −16 Ry. For water ecut = 120Ry, econv = 10 −14 Ry. We used a cubic simulation cell of 20 and 30Å for Argon and water respectively.

Appendix E. Numerical stability of QEHeat
In the following appendix we show the stability and convergence properties of a QEHeat calculation on a snapshot of 125 water molecules.  substantially up to dt = 2∆t. At this value of dt, even reducing it by a factor of 10 would change the component of the current by less than 0.001%, a negligible error given that a typical value of the thermal transport coefficient has an accuracy of 10% [20,22,24]. Moreover dt = 2∆t would be beneficial in an on-the-fly computation, allowing to reuse the same wavefunctions computed in the MD simulation, neglecting the need for the recomputation of the scf cycles. The errorbars in the figure are computed using the testing feature, provided with QEHeat, presented in Section 5 and averaging 20 fluxes obtained from different initial random wave-functions. The previous calculations were performed with the self-consistent threshold econv equal to 10 −14 Ry. The lower, the better the quality of the wavefunctions calculated. In order to further test and prove the stability of the MUB current, we used again the aforementioned testing feature to compute the dt dependence of the statistical uncertainty of the MUB current, at different values of econv. Fig. E.6 shows the percentage error of the x-component of the flux for a specific snapshot of 125 molecule. Even though the dependence on dt is similar for the two values of econv, it is clear that increasing econv, at a fixed dt, increases the statistical error by orders of magnitude, three orders when econv is changed from 10 −14 Ry to 10 −08 Ry and at any fixed dt. Note also that the variance decreases when increasing dt, in an exponential way. It must be stated that the errors in In Section 4 we showed how QEHeat implements numerical derivatives with the symmetric approach. The code allows also, simply changing the three_point_derivative keyword to false, to compute the derivative within a non-symmetric approach. In this case the implemented expression is: The errorbars are obtained repeating the computation of the current with a different initial wavefunction. The picture also shows that for the non-symmetric derivative there are only slightly stronger non-linear behavior on the current for large dt.
• Simulation cell sizes. One has to take into account that strongly harmonic systems, like crystalline solids at ambient temperature, may require large simulation cells to remove boundary effects. For high temperature or high pressure solids this should be a smaller issue. Anharhmonic effects in disordered systems or glasses can also reduce the typical path lengths and reduce therefore boundary effects.
• Required accuracy. All settings should be tuned according to the desired accuracy on the thermal conductivity coefficient, which may vary according to the application.
We provide here the overall computational cost of the computation presented in section 6 for liquid water, as a guideline, even if system specific setups are suggested. For this test, we chose conservative parameters and the default symmetric derivative scheme. We ran the MD calculations with 4 nodes, 192 processors, whereas the QEHeat calculations, being trivially parallelizable, where each run on 1 node and 48 processors. All times here reported are multiplied by the number of processors and identify therefore the total cost of the computation. For the 64-molecules system, a general MD step costs ≈ 1.2 × 10 2 s of cpu time (average over ≈ 2.5 × 10 5 steps, leading to a trajectory 240ps long). The cost of evaluating J MU B for a single snapshot was ≈ 1.5 × 10 3 s of cpu time, in which 84% of the time is reserved to the solution of the linear system in Eq. (4) during the evaluation of J KS . The 3 minimizations needed for the symmetric derivative scheme cost all together ≈ 2.2 × 10 3 s (mean over all PW calculations). We evaluated J MU B every 20 MD timesteps, thus the overall overhead was around 1.5 times the cost of the whole ab initio molecular dynamics simulation.
To conclude, at the present state our code recomputes the self-consistent cycles for each step, adding some extra time that we considered in this analysis. We are currently working on using directly the wafefunctions provided during the Car-Parrinello molecular dynamics (CP-MD) simulation, a promising feature considering that our analysis on the dt dependence implies that the discretization step can be chosen equal to the CP-MD time step, and it will be available in a future release. All the computations are done on the Tier-0 system called Marconi CINECA, which have 3188 nodes equipped with 2 Intel Xeon 8160 (SkyLake) at 2.10 GHz with 24-cores each [28].