Method to manage integration error in the Green-Kubo method

The Green-Kubo method is a commonly used approach for predicting transport properties in a system from equilibrium molecular dynamics simulations. The approach is founded on the fluctuation dissipation theorem and relates the property of interest to the lifetime of fluctuations in its thermodynamic driving potential. For heat transport, the lattice thermal conductivity is related to the integral of the autocorrelation of the instantaneous heat flux. A principal source of error in these calculations is that the autocorrelation function requires a long averaging time to reduce remnant noise. Integrating the noise in the tail of the autocorrelation function becomes conflated with physically important slow relaxation processes. In this paper we present a method to quantify the uncertainty on transport properties computed using the Green-Kubo formulation based on recognizing that the integrated noise is a random walk, with a growing envelope of uncertainty. By characterizing the noise we can choose integration conditions to best trade off systematic truncation error with unbiased integration noise, to minimize uncertainty for a given allocation of computational resources.

Figure 1

Panel (a) shows the HCACFs (the decaying functions) plotted along side their integrals (the curves that rise to a plateau) computed from nine separate simulations of a 10 648-atom, perfectly crystalline, and periodically contiguous block of graphite. The data were taken from a study to determine the influence of Wigner defects on thermal transport in graphite [22]. The dashed lines correspond to the heat flux along the $\left[2\overline{1}\overline{1}0\right]$ direction and the solid lines correspond to the heat flux along $\left[01\overline{1}0\right]$. The system was found to be converged for size, and $\kappa$ is expected to be the same in both directions along the basal plane. This plot illustrates the increasingly diverging noise of the HCACF integrals, present even after 50 ps. To the eye, the ACFs look nicely converged after 10–15 ps. Plot (b) shows the gradual convergence of the HCACF with increasing averaging time during a single simulation. The amplitude of the fluctuations in the tail of the HCACF decays over time, but it is notable that continued averaging does not remove the pattern of the fluctuations.

Figure 2

Panel (a) corresponds to the step or velocity fluctuations that give rise to a random walk; in panel (b) a set of 10 random walks is shown in black and the expected root mean square translation distance at time $t$ is plotted in red; and panel (c) is the distribution of the random walks shown in panel (b). Panel (d) corresponds to the tail of a HCACF, depicting the noise fluctuations that integrate to a large error akin to a random walk, shown in panel (e) for all HCACF tails. The black lines correspond to heat-flux measurements along the $x$ direction ($\left[2\overline{1}\overline{1}0\right]$), and the blue ones are along the $y$ direction ($\left[01\overline{1}0\right]$). Both values were measured along the basal plane, and this distinction should not matter. The data set is explained in the Methods section. For the selected $20–50\mathrm{ps}$ interval, the distribution of all data points across the multiple simulation tails is shown in panel (f). A 1-ps moving average was used along with a peak find algorithm to plot major peaks in the HCACF tails, as shown in panel (d), in magenta. The peak distribution for all data is offered in panel (g). The dashed red lines in panels (f) and (g) correspond to a normal distribution with the standard deviation of each of the distributions and mean zero. A normal distribution with the mean for each of the data sets is shown in the solid lines for each case.

Figure 3

The noise of a HCACF tail in the $30–50\text{ps}$ interval is shown (a) decomposed into high-frequency (blue) and low-frequency (red) noise. The autocorrelations of the noise (black) and the high-frequency (blue) and low-frequency (red) components of the noise are shown in panel (b), along with fits through the high-frequency (cyan) and the low-frequency (magenta) autocorrelations. In panel (c) the integrated tail appears in black and the uncertainty envelope for $\delta t$ equal to the interval of the HCACFs is shown in dashed red; the uncertainty envelopes corresponding to the high-frequency and low-frequency noise are in cyan and magenta, respectively. The dashed black line that follows along the magenta is the combined uncertainty envelope of the high- and low-frequency noises, i.e., the square root of the sum of their squares.

Figure 4

This graph shows the application of multiple pass filters to isolate existing frequencies in the HCACF noise. The first filter applied selects out data below a $0.04-\text{ps}$ interval (the blue high-frequency line at the bottom of the graph) and leaves the remaining frequencies. The next filter has a $0.08-\text{ps}$ window and is used to filter the low-frequency data remnant from the first pass. This procedure is performed for $0.04-\text{ps}$ intervals up to a filter with a $0.56-\text{ps}$ window.

Figure 5

In panel (a) the tail of the HCACFs, their integral, the uncertainty envelope (cyan) calculated as described in the text, and its error (blue) are all plotted. In the inset in panel (b), instead of only considering the noisy tails of the HCACFs, the whole HCACFs are represented. In both panel (a) and the inset in panel (b) the solid black lines correspond to results along the $y$ direction, and the dashed black lines correspond to results along the $x$ direction. The bold red line in the inset in panel (b) is the integral of the average of the HCACFs; the solid green line is the standard error computed for the 18 HCACF integrals; and the dashed green line is the standard error of the 216 $50-\mathrm{ps}$ HCACF integrals that can be obtained from the 18 sets of data with $600\mathrm{ps}$ each. These lines are shown in the inset in panel (b) for perspective, but also in the larger plot in panel (b) for a clearer distinction between them and the cyan line, which shows the uncertainty calculated as described in the text, using the random walk approach.

Figure 6

Panel (a) is the normal distribution over all $\mathbf{J}$. Panel (b) is the distribution of the noise from the tails in the $30–50\mathrm{ps}$ interval. Panel (c) is the distribution of the peaks fit to the noise from the tails in the $30–50\mathrm{ps}$ interval, as shown in Fig. 2.

Figure 7

Panel (a) shows two extremes both in terms of their total integrated value and the interval, ${\tau }_L$, of their low-frequency oscillations. The uncertainty envelope for the integrated HCACF in purple is slightly above the maximum standard error (blue), whereas that of the HCACF integral in brown is below. The corresponding noise and noise integrals for these extrema are shown in panel (b).

Figure 8

Panel (a) shows the averaged HCACFs for all simulations along $x$ (cyan) and $y$ (magenta), the HCAFCs for $x$ (blue) and $y$ (red) for the large, $8-\mathrm{ns}$, simulation and the corresponding integrals in the same color. To observe the effect of a single outlier, all HCACFs except the purple one (see Fig. 7) are averaged. The resulting HCACF and integral are plotted in dashed yellow. Panel (b) shows the integrals [using the same color scheme as in panel (a)] with the corresponding uncertainty envelope around them.

Figure 9

This shows the integrated HCACF average for all simulations along $x$ (cyan) and $y$ (magenta) for the subset of $800-\mathrm{ps}$ simulations resulting from the $8-\mathrm{ns}$ simulation, the integrated HCAFCs for $x$ (blue) and $y$ (red) for the large, $8-\mathrm{ns}$, simulation and the corresponding uncertainty envelope around them.

Figure 10

In panel (a), in addition to the HCACF, the moving average of the uncertainty envelope computed using the random walk approach is also propagated through the simulation time. In panel (b) the percentage error is computed as the uncertainty envelope over the total integral.

Figure 11

The heat flux (black), $\mathbf{J}$, a $0.4\text{-ps}$ moving average of $\mathbf{J}$ (red), and a gradual cutoff of the higher peaks of $\mathbf{J}$ (green) are shown in panel (a). The HCACF and integral for each of the above cases is shown in panel (b) as is, and is normalized in panel (c). Panels (d) and (e) are close-ups of panels (b) and (c), respectively. The color coding is maintained throughout the figures.

Figure 12

Panel (a) shows $\mathbf{J}$ (black), a transform on $\mathbf{J}$ that keeps its higher peaks and replaces data between the peaks with a zero value (blue), and a line at $550\text{ps}$ representing a cutoff of the $\mathbf{J}$ data above it. Panel (b) shows the normalized HCACF for the above cases, including those depicted in Fig. 11. The HCACF as is shown in panel (c). The color code is kept constant between Figs. 11 and 12.