*2.3. Error Analysis*

The output of the MD simulations is the stress tensor as a function of time for a given thermodynamic state. Generally, the observations in this sequence are correlated for time differences that are shorter than some intrinsic correlation time. Therefore, we cannot use the sample standard deviation for error analysis. Instead, we need to apply a more sophisticated analysis applicable to correlated data. We follow the approach of Flyvbjerg and Petersen [28], which requires calculation of *c*(*t*), the time correlation function for each dataset (stress component in our case), for data obtained after *n* time steps (*n* + 1 data points). The correlation function is given by

$$\mathcal{L}(t) = \frac{1}{n-t} \sum\_{k=1}^{n-t} (\sigma\_i(t) - \langle \sigma\_i \rangle)(\sigma\_i(k+t) - \langle \sigma\_i \rangle),\tag{7}$$

where angular brackets denote the time average. Then, the error *s* is calculated using the following equation:

$$s^2 = \frac{c(0) + 2\sum\_{t=1}^{t'} (1 - \frac{t}{n})c(t)}{n - 2t' - 1 + \frac{t'(t'+1)}{n}}.\tag{8}$$

Here, *t* is a cut-off time chosen to be much larger than the correlation time of *c*(*t*) but much smaller than *n*. We observed typical relaxation times for *c*(*t*) of about 5000 time steps (1 ps) so we used *t* = 50,000 (10 ps). The errors for the elastic coefficient values reported in the next section are reported as ±*s*.
