The effects of the finite size of the simulation box in equilibrium molecular dynamics simulations are investigated for prototypical superionic conductors of different types, namely, the fluorite-structure materials PbF2, CaF2, and UO2 (type II), and the α phase of AgI (type I). Largely validated empirical force-fields are employed to run ns-long simulations and extract general trends for several properties, at increasing size and in a wide temperature range. This work shows that, for the considered type-II superionic conductors, the diffusivity dramatically depends on the system size and that the superionic regime is shifted to larger temperatures in smaller cells. Furthermore, only simulations of several hundred atoms are able to capture the experimentally observed, characteristic change in the activation energy of the diffusion process, occurring at the order–disorder transition to the superionic regime. Finite-size effects on ion diffusion are instead much weaker in α-AgI. The thermal conductivity is found generally smaller for smaller cells, where the temperature-independent (Allen-Feldman) regime is also reached at significantly lower temperatures. The finite-size effects on the thermal motion of the non-mobile ions composing the solid matrix follow the simple law that holds for solids.
I. INTRODUCTION
Superionic (SI) materials are characterized by a matrix of atoms arranged in a (crystalline or amorphous) solid, and of one (or more) mobile species, which diffuses above some critical temperature. The interest for SI materials has largely increased during the past decades, along with the quest for good candidates in the realization of solid-state batteries, where charge carriers like lithium ions move through a solid-state electrolytic matrix.1–4 Moreover, the superionic phases of water and ammonia5 have been predicted to compose a large fraction of the outer cores of ice giant planets, such as Uranus and Neptune,6,7 and many recent theoretical8,9 and experimental10,11 studies have focused on transport properties of materials becoming SI at planetary conditions, to study the evolution of these celestial bodies.12,13
The complexity and the variety of new SI materials naturally imply that, from the computational material-science standpoint, much effort is devoted to the prediction of static and dynamical properties and their microscopic description by means of atomistic simulations. In particular, molecular dynamics (MD) simulations are needed whenever dynamical properties (such as transport coefficients and correlation functions) are investigated.14
Due to the chemical complexity of many of these materials, ab initio MD simulations are often performed, whose computational cost currently limits the simulation box to a few hundred atoms at most. This limitation poses serious questions on the role of finite-size effects (FSE) in the characterization of the physical properties of a SI material. For instance, recent tests were run on Li10GeP2S12-type SI conductors,15 where machine-learning interatomic potentials trained on ab initio calculations allowed for simulations that are inaccessible to ab initio MD.16 These calculations showed that simulation boxes containing even some hundred atoms overestimate the Li-ion diffusivity by one order of magnitude with respect to the largest size considered (1600 atoms). This may have dramatic consequences in calculations aiming to find the best SI conductors for realistic devices.17–19
In liquids, FSE affecting particle diffusion have been extensively investigated (see, e.g., the recent review of Ref. 20): The hydrodynamics arguments by Yeh and Hummer21,22 suggest that, for a cubic simulation box under periodic boundary conditions (PBC) and for a given particle density, the diffusivity of the liquid can be corrected by a factor proportional to the inverse of the box side, and that the proportionality coefficient only depends on geometric factors, the temperature, and the viscosity of the liquid, which is usually largely independent of the size.22,23 Recent works evidenced that the application of the Yeh and Hummer correction is justified also for multicomponent liquids and ionic melts.24–26 Nevertheless, the hydrodynamics equations of an SI material are different from (and more complicated than) those of simple liquids.27,28 For instance, transverse modes of the lattice survive, the static shear modulus is non-vanishing, and the atoms of the mobile species diffuse via hopping mechanisms that are qualitatively different from the motion of particles in a simple fluid:29 the Yeh–Hummer arguments are, therefore, inappropriate to account for FSE in the diffusion of charge carriers in SI materials. Furthermore, while general trends for the FSE on heat transport in solids and liquids have been reported in the literature,31–32 such an analysis is currently missing for equilibrium MD simulations of thermal conduction in SI materials.
This article aims at investigating of the FSE in the calculation of relevant static and dynamical properties of SI materials via equilibrium MD simulations. I focus on simple yet paradigmatic examples of type-I and type-II SI conductors, which can be effectively described in terms of largely validated empirical potentials. I restrict my analysis to systems with perfect stoichiometry. Furthermore, I consider cubic simulation boxes where the unit cell, of lattice parameter a, is equally replicated ℓ times in all the three spatial directions, to avoid additional and non-trivial effects that arise, even for simple liquids, in the case of anisotropic replications.33
Section II discusses the SI materials I selected to investigate. Section III gives methodological details on the equilibrium MD simulations that I performed. Section IV provides the main results of the calculations, by analyzing the size dependence of the specific heat capacity, the mobile-ion diffusivity, the thermal conductivity, and the Debye–Waller B-factors, for each of the selected SI materials. Finally, I draw general conclusions in Sec. V.
II. DISCUSSION
I choose the fluorite-structure materials PbF2, CaF2, and UO2, as simple, yet prototypical examples of SI materials where FSE should be particularly relevant: all the energetically equivalent regular sites of the mobile species—the anions—are occupied, and the hopping of one diffusing anion can eventually occur only with a net hopping of other anions, since anion diffusion “occurs by discrete hops between regular sites,” and anions “do not reside in a well-defined manner on the cube-centre sites” (verbatim from Ref. 29. See also Ref. 34 for an insightful analysis of anion distribution in fluorites, and Ref. 35 for a recent, comprehensive study of cooperative F dynamics in PbF2). Such a concerted hopping mechanism may easily extend to more than one lattice constant, leading to a size dependence. Furthermore, these materials (see Fig. 1, left panels) are characterized by a continuous order–disorder transition to the SI phase with no structural change in the crystalline structure of the non-diffusive species (type II superionic materials). The diffusion mechanism depends on the specific temperature regime, and, in particular, whether the system is in the SI phase or not.36
Structure and behavior of the diffusivity, D, of the mobile species against inverse temperature for fluorite-structure materials (left) and α-AgI (right). The (non-)mobile ions are represented in (blue) red. The empty red circles indicate degenerate tetrahedral positions. Rear-faces’ atoms are not displayed. The shaded yellow area indicates the superionic regime.
Structure and behavior of the diffusivity, D, of the mobile species against inverse temperature for fluorite-structure materials (left) and α-AgI (right). The (non-)mobile ions are represented in (blue) red. The empty red circles indicate degenerate tetrahedral positions. Rear-faces’ atoms are not displayed. The shaded yellow area indicates the superionic regime.
The finite size of the sample is known to cause a shift in the critical temperature of second-order phase transitions.37 Therefore, for a given temperature and particle density, a small simulation box may be in a different thermodynamic phase with respect to a larger one, with different diffusive mechanisms and a dramatic effect on diffusion. This is likely to be the case in all those materials where the diffusion mechanisms are strongly dependent on T. As we shall see in Sec. IV, this tangling between the diffusion mechanisms (hydrodynamics) and the phase of the system (thermodynamics) is responsible for dramatic FSE on ionic transport in these materials.
I also investigate a different material, the cubic phase of silver iodide (α-AgI), as a typical example of a system where FSE on the diffusion coefficient should be less relevant (Fig. 1, right panels): in contrast with fluorite-structure materials, in α-AgI the large degeneracy of equivalent positions that one Ag ion—the diffusive species—may take inside a unit cell results in a large freedom in the choice of empty sites (empty red circles in Fig. 1) that a selected Ag ion can hop to: the temperature affects the probability that hopping occurs, but not the general mechanism of diffusion. Moreover, due to the large degeneracy of empty regular sites, there is no need for the hopping of one Ag cation to be accompanied by the hopping of other neighbor Ag cations. The α-AgI phase is superionic and is reached after a first-order phase transition (at K at ambient pressure38) from the hexagonal, non-conducting β-phase. The sudden, discontinuous change in the ionic diffusion at the phase transition makes AgI a type-I SI conductor.
III. METHODOLOGY
A. Empirical interatomic potentials
The choice of these materials is also motivated by the availability, in the scientific literature, of largely validated empirical potentials that proved to qualitatively describe the ion diffusion mechanisms as well as the static properties of these systems.39 These potentials make it possible to run, at a feasible computational cost, reliably long simulations (∼ns) at different temperatures and sizes, to extract general behaviors.
For PbF2 and CaF2, I employ the following two-body potential, combination of Coulomb and Buckingham potentials:
with the parameters optimized in Ref. 40 (PbF2) and in Ref. 29 (CaF2), reported in Table I. The success of this potential in the microscopic study of ionic diffusion in these materials dates back to the 1980s. For UO2, I employ a recently developed potential, described in Refs. 41 and 42, which combines Buckingham, Coulomb, and Morse potentials to treat two-body interactions, as well as the embedded atom method (EAM) to account for many-body interactions. I refer the interested reader to the original literature and to the Materials Cloud repository of this work for the explicit parametrization values (see “Data Availability”).
Parameters employed in the potential of Eq. (1). For cation–cation interaction, only the Coulomb interaction is considered (i.e., A++ = C++ = 0, with + indicating Pb or Ca). The integer charges zPb = zCa = +2 and zF = −1 are used.
. | Aij (eV) . | ρij (Å) . | Cij (eV Å6) . |
---|---|---|---|
Pb–F | 122.7 | 0.516 | 0.0 |
F–F | 10 225.0 | 0.225 | 107.3 |
Ca–F | 674.3 | 0.336 | 0.0 |
F–F | 1808.0 | 0.293 | 109.1 |
. | Aij (eV) . | ρij (Å) . | Cij (eV Å6) . |
---|---|---|---|
Pb–F | 122.7 | 0.516 | 0.0 |
F–F | 10 225.0 | 0.225 | 107.3 |
Ca–F | 674.3 | 0.336 | 0.0 |
F–F | 1808.0 | 0.293 | 109.1 |
Finally, for AgI, I use the following combination of Coulomb and Morse potentials:
with the parameters of Ref. 43, reported in Table II, derived via the Chen–Möbius lattice inversion method from ab initio calculations of cohesive energies. This potential displays good agreement with experiments44 concerning static properties of different phases, as well as Ag diffusivity in α-AgI, and it is also consistent, in a wide temperature range, with the widely used Parrinello–Vashishta–Rahman empirical potential.45
Parameters of the Morse component of the potential in Eq. (2). For cation–cation interaction, only the Coulomb interaction is considered (i.e., DAg,Ag = 0). The fractional charges zAg = −zI = 0.3181 are used for the Coulomb term.
. | Dij (eV) . | αij (Å−1) . | Rij (Å) . |
---|---|---|---|
I–Ag | 0.55 | 1.600 | 2.6 |
I–I | 0.16 | 0.684 | 5.7 |
. | Dij (eV) . | αij (Å−1) . | Rij (Å) . |
---|---|---|---|
I–Ag | 0.55 | 1.600 | 2.6 |
I–I | 0.16 | 0.684 | 5.7 |
B. Details on MD simulations
All simulations are performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software.46 It has the great advantage, with respect to other MD codes, that force computation is not subject to minimum image conventions, and one can use cutoffs larger than half the simulation domain size, thanks to the inclusion of “ghost” atoms. This is particularly important for the purpose of this work, where small boxes are needed for the FSE analysis, yet the cutoff radius should be the same to avoid changes in the form of the potential. The long-range interactions are included by means of the Ewald-summation technique in MD simulations of PbF2, CaF2, and AgI, and with the particle-particle particle-mesh (PPPM) method for the MD simulations of UO2.47 The simulations of PbF2, CaF2, and α-AgI are run with a MD timestep of 4 fs. For UO2, the MD time step is set to 2 fs. The trajectories ( ps), from which the mean square displacements of the atoms are computed, are sampled each 10 MD time steps. For a given material, all constant-NVT (canonical) and constant-NVE (microcanonical) simulations (N is the number of particles, V is the cell volume, T is the temperature, and E is the total energy) are run at fixed lattice constant, a, irrespective of the temperature, i.e., no thermal expansion is considered for simplicity. I employ the following lattice constants: aPbF2 = 6.056 Å (value at T = 792 K in Ref. 40), aCaF2 = 5.712 Å,29 aUO2 = 5.65 Å (online material of Ref. 42), and aαAgI = 5.37Å.43 The NVT simulations are run with a Bussi–Donadio–Parrinello stochastic-velocity-rescaling (SVR) thermostat,48 as implemented in LAMMPS. The temperature damping parameter of the SVR thermostat is set to 100 MD timesteps. Further details on the equilibration procedures and on the cutoffs employed are reported, for the sake of reproducibility, in the input scripts of the simulations, which are all available in the Materials Cloud repository of this work (see “Data Availability”).
IV. RESULTS
I proceed by investigating FSE for physical quantities important for superionics, namely, the specific isochoric heat capacity (Sec. IV A), the diffusivity of the mobile species (Sec. IV B), the thermal conductivity (Sec. IV C), and the Debye–Waller B-factor of the non-diffusive ions of the solid matrix (Sec. IV D).
A. Specific heat capacity and critical temperature
The isochoric molar specific heat capacity, cV, is obtained from the finite-difference derivative of the average energy with respect to the temperature and displayed vs T in Fig. 2 for the four materials considered. Let us first focus on the fluorite-structure materials (first three panels of Fig. 2). For sufficiently large cells (ℓ ≥ 2, i.e., N ≥ 96), the heat capacity clearly displays a peak. Experimentally, a peak in the heat capacity—at some high critical temperature, yet below the melting point—has been observed via heat-content measurements of fluorite-structure materials, and associated with a transition that is not of the first order; such anomaly accompanies a sensible onset of electrical conduction, indicating a transition to the superionic phase.49,50 At a large enough size, the calculated cV is in fairly good agreement with experiments. As mentioned in Sec. II, it is known since the late 1960s that the effect of a finite size is to broaden a second-order transition and to shift the (pseudo)critical temperature Tc(ℓ) with respect to its thermodynamic-limit value Tc(∞).51 Whether the shift is positive or negative depends, among other factors, on the boundary conditions: usually, in PBC, Tc(ℓ) > Tc(∞), as a result of extra “communication” via paths that encircle the torus (verbatim from Ref. 51). This is, in fact, the behavior observed in Fig. 2, where the peak shifts toward lower temperatures as the size of the simulation is increased, in agreement with the existing literature.52 In the smallest cell, ℓ = 1 and N = 12, no peak is observed and, in the temperature range that I consider, the heat capacity is always sensibly lower than the one obtained with larger simulation boxes. I remark that, in agreement with Ref. 42, the melting of UO2, predicted for this potential at K41 via a moving interface method, is not observed in these simulations, where the lattice parameter is kept fixed to its value at K53 (see Sec. III B).
Isochoric specific heat capacity, cV, as a function of temperature and system size for the fluorite-structure materials investigated. The markers indicate the calculated quantities, while the line is a three point running-average window filter. The shaded area indicates the statistical uncertainty, obtained from the standard block analysis.
Isochoric specific heat capacity, cV, as a function of temperature and system size for the fluorite-structure materials investigated. The markers indicate the calculated quantities, while the line is a three point running-average window filter. The shaded area indicates the statistical uncertainty, obtained from the standard block analysis.
In striking contrast with fluorite structure materials, for α-AgI, no significant size effect is observed, in line with its intrinsic superionic structure. The link between FSE and the onset of a superionic phase transition is clearly highlighted from the analysis of the diffusivity of the mobile species, as discussed below.
B. Diffusion coefficient
The diffusion coefficient of the mobile species is the most characterizing quantity of superionic conductors: its behavior at the superionic transition dictates the classification of SI materials (see Fig. 1). Ion diffusion is usually described as an Arrhenius-like process,
where k is Boltzmann’s constant, A is a prefactor with the dimensions of a diffusivity, and Ea is the activation energy of the hopping mechanism leading to particle diffusion. In SI materials, in general, for a specific regime of diffusion, A and Ea weakly depend on the temperature: it makes sense, therefore, to plot the logarithm of the diffusion coefficient against the inverse temperature (in the so-called Arrhenius plot) to highlight significant changes in the diffusion mechanism. Ea represents the slope of the Arrhenius plot, and A is its intercept.
Figure 3 displays the Arrhenius plots of the diffusivity of the mobile species for the materials considered in this work, at different system sizes. The reported values for D (markers) are obtained from the slope of the mean square displacement (MSD) of the mobile species at large enough time,
where Nd is the number of atoms of the diffusive species and ri(t) is the position of the ith atom of the mobile species at time t. A linear fit on the MSD at large t was used to obtain D. A block analysis, with four blocks of ps each, is performed to extract the uncertainty on the MSD. This uncertainty propagates to the asymptotic-time slope of the MSD and eventually to the diffusion coefficient. The uncertainty on D is reported as the shaded area in the Arrhenius plots. The reported D is computed in the “laboratory” reference frame, where the center of mass of the non-mobile species is fixed, although the MD simulations are performed in the barycentric reference frame, where the total momentum vanishes. The relations between the different, reference-frame dependent definitions of the diffusion coefficients are described in the supplementary material, Sec. S2 A. To check whether the SVR thermostat affects the estimate of D, I also run NVE simulations, previously equilibrated at the target temperature. Diffusion coefficients extracted from NVE and NVT simulations display fully compatible Arrhenius plots for each simulation cell size, as reported in the supplementary material, Fig. S9.
Arrhenius plot of the diffusion coefficient of the mobile species for the fluorite-structure materials considered in this work, as a function of the cell size. The markers indicate the calculated D and the shaded area indicates the uncertainty from a block analysis on four blocks of ps each. In the plots for fluorite-structure materials, the solid (dashed) lines are fits to Eq. (3) above (below) Tc(ℓ). The dotted vertical lines indicate the size-dependent critical temperature, Tc(ℓ), obtained from the position of the maximum of the heat capacity (same color code).
Arrhenius plot of the diffusion coefficient of the mobile species for the fluorite-structure materials considered in this work, as a function of the cell size. The markers indicate the calculated D and the shaded area indicates the uncertainty from a block analysis on four blocks of ps each. In the plots for fluorite-structure materials, the solid (dashed) lines are fits to Eq. (3) above (below) Tc(ℓ). The dotted vertical lines indicate the size-dependent critical temperature, Tc(ℓ), obtained from the position of the maximum of the heat capacity (same color code).
In PbF2, CaF2, and UO2, I find that D(T) strongly depends on the system size. For instance, the diffusivity for N = 96 is greater (smaller) than the large-N values at low (high) T, in agreement with previous observations.29 Furthermore, the kink in the Arrhenius plot at Tc, characterizing type II materials, is observable only for N ≥ 324 (i.e., ℓ ≥ 3). This is quantified by the values of Ea(ℓ) (the slope of the Arrhenius plot) extracted from the fit of D to Eq. (3) below and above the size-dependent critical temperature, Tc(ℓ), and reported in Table III: only for N ≥ 324, a sensible difference between Ea(T < Tc) and Ea(T > Tc) is observed. Notice, as well, that the kink in the Arrhenius plot, associated with the SI transition, moves toward lower temperatures as the size is increased, in accordance to the existing literature,54 and with the shift in Tc extracted from the maximum of the heat capacity (Fig. 2).
Activation energies, Ea (in meV), of the considered fluorite-structure materials, as a function of the system size, from the weighted fit to Eq. (3). For each system, the first row corresponds to the regime T < Tc(ℓ), and the second row to the regime T > Tc(ℓ). For N = 12 (i.e., ℓ = 1), where no phase transition is observed, only one value is provided. See the Materials Cloud Repository for data and details of the fit.
N = . | 12 . | 96 . | 324 . | 768 . |
---|---|---|---|---|
PbF2 | 160 ± 9 | 320 ± 20 | 590 ± 20 | 620 ± 20 |
390 ± 40 | 324 ± 11 | 308 ± 8 | ||
CaF2 | 560 ± 50 | 1300 ± 60 | 1930 ± 40 | 1830 ± 60 |
1070 ± 330 | 940 ± 70 | 750 ± 40 | ||
UO2 | 1020 ± 80 | 2490 ± 70 | 4110 ± 80 | 4160 ± 50 |
2410 ± 150 | 1490 ± 70 | 1240 ± 40 |
N = . | 12 . | 96 . | 324 . | 768 . |
---|---|---|---|---|
PbF2 | 160 ± 9 | 320 ± 20 | 590 ± 20 | 620 ± 20 |
390 ± 40 | 324 ± 11 | 308 ± 8 | ||
CaF2 | 560 ± 50 | 1300 ± 60 | 1930 ± 40 | 1830 ± 60 |
1070 ± 330 | 940 ± 70 | 750 ± 40 | ||
UO2 | 1020 ± 80 | 2490 ± 70 | 4110 ± 80 | 4160 ± 50 |
2410 ± 150 | 1490 ± 70 | 1240 ± 40 |
As shown in the supplementary material, Figs. S4 and S5, all these results are in qualitative agreement also with frozen-matrix simulations,55 where the non-diffusive ions of the solid matrix are kept fixed to their equilibrium position, and with simulations employing a short-range version of the Coulomb interaction.56 This agreement confirms the proposed picture whereby the FSE on the diffusivity of mobile ions are mainly imputable to geometric factors (such as the degeneracy of empty sites and the extra interatomic communication paths occurring in PBC), rather than to the details of the potential or the vibrations of the solid matrix. This is also justified by a set of simulations run on non-stoichiometric lead fluoride, where a size-independent concentration of empty sites is generated by removing a set of randomly selected F− ions accordingly: the presence of available empty sites favors hopping on a more local scale than in systems with perfect stoichiometry, where, below the SI transition, all the regular sites are occupied and the hopping of one F− anion can only occur with a net hopping of other F− anions, and results in a drastic reduction of FSE on fluorine-ion diffusion, as shown quantitatively in Appendix B.
In striking contrast with fluorite-structure materials, in α-AgI, the values of the diffusivity of Ag ions at different cell sizes (even for very small cells) are all consistent with each other. The activation energies range from 93 ± 3 to 100 ± 7 meV (without any particular trend connected to the box size) and agree with the experimental value Ea = 94.97 meV of Ref. 44. As suggested above, I ascribe the size independence of the diffusivity to the large degeneracy of equivalent Ag ion sites within the conventional unit cell of α-AgI: in this superionic phase, Ag-ion hopping occurs on a smaller length scale than in fluorite-structure materials, where, instead, the hopping of one ion can occur only if accompanied by the concerted hopping of other ions, since—apart from the short transient of the jump—all the regular sites are occupied, as confirmed also in the literature (see, e.g., Refs. 40 and 57).
C. Thermal conductivity
Figure 4 shows the temperature and size dependence of the thermal conductivity, κ, for the materials considered, extracted from NVT MD simulations according to the Green–Kubo (GK) theory of linear response for multicomponent, diffusive systems, as described in Appendix A. Even if a mode-based analysis58,59 of heat conduction would be required for a quantitative assessment of the role of disorder on phonon propagation, here, I employ the GK formula and MD simulations to include the role of diffusing ions, which, by construction, is not considered in mode-based calculations.
Temperature behavior of thermal conductivity (markers) and associated uncertainty (shaded area) from cepstral analysis. The vertical dotted lines indicate the size-dependent critical temperature to the SI phase, obtained from the position of the maximum of the specific heat capacity, see Fig. 2.
Temperature behavior of thermal conductivity (markers) and associated uncertainty (shaded area) from cepstral analysis. The vertical dotted lines indicate the size-dependent critical temperature to the SI phase, obtained from the position of the maximum of the specific heat capacity, see Fig. 2.
Let us first focus, again, on fluorite-structure materials. At large T, κ is almost independent of T due to the increasing disorder that accompanies the diffusion of mobile ions, suppressing phonon propagation. This is a manifestation of the Allen-Feldman regime,60 and a typical feature also observed in other SI materials, such as solid-state electrolytes.61 The FSE of κ are system dependent, the values for PbF2 being almost converged for N = 96, in contrast with CaF2 and UO2. Nonetheless, a general trend can be observed. The N = 12 (ℓ = 1) cell dramatically fails in describing heat transport and displays a non-physical, slight increase of κ with temperature. For ℓ > 1, the Allen-Feldman limit is achieved at lower temperatures for smaller simulation boxes, in agreement with the larger diffusivity at low T observed for small size, and the higher degree of disorder in smaller cells.62 Furthermore, the presence of a single defect is not supposed to strongly perturb a large system, while it would dramatically affect a few-atom cell, suppressing phonon propagation in favor of the Allen-Feldman limit. The activation of diffusion with the onset of significant disorder is much more relevant for heat transport than the actual SI phase transition: at T = Tc(ℓ) (vertical dotted lines), no particular feature of κ(T) is, in fact, observed.
Overall, the large-cell values are in good agreement with existing literature: κPbF2 = 1.4 Wm−1 K−1 at 300 K, Ref. 63 (experimental); κCaF2 = 1.46 ± 0.29 Wm−1 K−1 at T = 1694 K, Ref. 65 (numerical simulation, same density, and potential); κUO2 = 1.5 Wm−1 K−1 at T = 3000 K, Ref. 64 (numerical simulation, close density but different potential). I remark that, in UO2, for T ≳ 2000 K, the growing contribution of electrons to heat conduction must be added for a comparison with experiments.
The thermal conductivity of α-AgI is almost constant in the considered temperature range. This indicates that phonon propagation is always suppressed by disorder in favor of the AF regime in this intrinsically superionic phase. All simulations with N ≥ 108 are fully compatible. For the small cell, N = 32, κ is lower, though only by 10%–15%, than the converged value. The results are in good agreement with the experimental value of 0.17 Wm−1 K−1 at T ≈ 500 K.65
D. Dynamics of the non-diffusive species
FSE also affect the dynamics and thermal vibrations of the solid matrix, which I investigate in this section in terms of the mean-square-displacement of the non-diffusive (n.d.) ions with respect to their equilibrium position , which is a time-independent quantity. Furthermore, following a standard convention, I shall employ the so-called B-factor, , entering the Debye–Waller factor that dictates the attenuation of x ray or neutron scattering in experiments.
For SI phases, MD simulations are needed to compute the B-factor, since normal-mode-based approaches66 cannot be applied. Nonetheless, the values obtained from MD simulations are known to slowly converge with size. Early calculations for cubic hard-sphere solids show FSE corrections that follow a 1/ℓ law.67 This law stems from the minimum frequency that can be sampled in a finite box size. An extensive derivation of this FSE is provided in the supplementary material, with an example on FCC solid argon. Interestingly, the same trend is observed in my simulations on SI materials. Figure 5 shows the B-factor of PbF2, as a function of 1/ℓ, at different temperatures. The linear behavior is evident, and a linear fit of the data can be used to extrapolate the B-factor for infinite size, B(∞) (blue crosses). Note that the B-factor is here re-scaled by T and other geometric factors so that the slope represents the inverse of an effective elastic modulus of the material (see the supplementary material, Sec. V A). Only the smallest box considered (i.e., the conventional unit cell, ℓ = 1, with N = 12 atoms) deviates from the 1/ℓ behavior, not surprisingly, as it fails in describing most of the quantities analyzed so far, and where higher order FSE may enter.
Finite-size effects affecting the Debye–Waller B-factor of Pb atoms in PbF2, at different temperatures. Values obtained from simulations (markers and error bars) and the fit B(ℓ) = m/ℓ + B(∞) (for ℓ ≥ 2) are displayed. The blue crosses indicate the extrapolated asymptotic value B(∞). The absolute value of the slope, |m|, is reported in the inset. The linear behavior in 1/ℓ is evident. For ℓ = 1, higher order FSE may occur.
Finite-size effects affecting the Debye–Waller B-factor of Pb atoms in PbF2, at different temperatures. Values obtained from simulations (markers and error bars) and the fit B(ℓ) = m/ℓ + B(∞) (for ℓ ≥ 2) are displayed. The blue crosses indicate the extrapolated asymptotic value B(∞). The absolute value of the slope, |m|, is reported in the inset. The linear behavior in 1/ℓ is evident. For ℓ = 1, higher order FSE may occur.
The 1/ℓ behavior is observed in all the systems considered. An exception is the temperature range between 2600 and 2900 K for UO2: in fact, the drastic change in the elastic properties of UO2, and, therefore, of B, at the SI phase transition, occurs at different critical temperatures for different sizes (see Sec. IV A). Figure 6 displays the B-factor as a function of temperature for different systems and sizes. The blue, shaded area represents B(∞) with its uncertainty. Note that significant FSE are here observed not only in fluorite-structure materials but also in α-AgI. In fact, the very reason why the B-factor exhibits FSE—the minimum mode frequency that can be probed for a given cell—is very general and not system dependent. Once again, these calculations show that wrong results may be obtained even from an accurate description of interatomic interactions in a MD simulation (e.g., by computing forces ab initio), whenever FSE are not correctly accounted for.
B-factors of the non-diffusive ions for the materials considered in this work, as a function of the system size. The vertical dotted lines indicate the size-dependent critical temperature to the SI phase, obtained from the position of the maximum of the specific heat capacity. The blue, shaded area represents the extrapolated value, accounting for the finite-size correction.
B-factors of the non-diffusive ions for the materials considered in this work, as a function of the system size. The vertical dotted lines indicate the size-dependent critical temperature to the SI phase, obtained from the position of the maximum of the specific heat capacity. The blue, shaded area represents the extrapolated value, accounting for the finite-size correction.
V. CONCLUSIONS
Finite size effects (FSE) in superionic materials are, in general, not negligible and strongly depend on the specific system and property examined. In materials that possess intrinsically a high availability of degenerate hopping sites, such as α-AgI, the FSE on the diffusivity—the key quantity of superionic conductors—are weak. Things change, instead, whenever the mechanism of diffusion strongly depends on the temperature range, like in fluorites. These systems display an order–disorder transition to the superionic regime without a net change neither in the lattice structure nor in the available hopping sites, which coincide with those occupied by the mobile ions in the non-superionic phase. This work shows that, in these materials, the diffusivity is strongly affected by the system size, even qualitatively: the change of activation energy of the diffusion process at the SI critical temperature, resulting in the typical kink in the Arrhenius plot of the diffusivity, can be observed only above a certain threshold size. The order–disorder critical temperature is also largely affected by FSE, being larger for smaller sizes, as indicated by the shift of the maximum of the specific heat capacity. Therefore, in these materials, tangling between hydrodynamics (diffusion processes) and thermodynamics (the specific phase of a material) is responsible for changes in the diffusivity of even some orders of magnitude, at a given temperature, depending on the size. The thermal conductivity is also affected by FSE, mainly due to the role of disorder in hindering phonon propagation: this effect is larger at smaller simulation boxes, where the Allen-Feldman regime is reached at lower temperatures. The transition to the SI phase does not seem to strongly affect κ, instead. In general, an a priori determination on the minimum size that is sufficient to obtain satisfactory results is hardly feasible, due to strong dependence of FSE on the specific superionic material. My analysis suggests, a posteriori, that, for fluorite-structure materials, a minimum of N = 324 atoms is needed to correctly capture at least the main qualitative features of particle diffusion and heat transport, while for the α-AgI phase, even a relatively small cell of N = 108 atoms is sufficient to obtain a quantitative convergence of all the analyzed properties.
FSE also affect the Debye–Waller B-factor incorporating the thermal motion of the non-diffusive ions of the solid matrix. The 1/ℓ (or 1/N1/3) behavior, predicted for simple solids, is observed also for the SI materials considered. This is motivated by the minimum frequency of oscillation, which can be captured by a simulation of a given size, irrespective of the specific system or phase considered.
Increasingly reliable interatomic potentials have provided, during the past few years, a systematically more accurate description of SI materials. However, this work clearly shows that accurate potentials alone are not sufficient to converge the calculation of many key properties, if not complemented with a full treatment of FSE. From the application point of view, particular attention must be paid in calculations aiming to compare different candidates for realistic devices, such as solid-state batteries. A given simulation size may be sufficient for some superionic materials but not for others.
SUPPLEMENTARY MATERIAL
See the supplementary material for more details about (i) the calculation of the heat capacity, (ii) the calculation of the diffusivity (role of the reference frame, of modified Coulomb interactions and of vibrations of the solid matrix), (iii) a comparison of the results obtained with NVT and NVE simulations, (iv) the calculation of the thermal conductivity, (v) the derivation of the 1/ℓ law for FSE of the B-factor and a simple application to solid argon, and (vi) the heat capacity of the defected structure of Appendix B.
ACKNOWLEDGMENTS
I thank Michele Ceriotti, Loris Ercole, Alfredo Fiorentino, Lorenzo Gigli, and Paolo Pegolo for insightful discussions and fruitful comments on the manuscript. I acknowledge funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Action IF-EF-ST (Grant Agreement No. 101018557; TRANQUIL).
AUTHOR DECLARATIONS
Conflict of Interest
The author has no conflicts to disclose.
DATA AVAILABILITY
The data that support the plots and relevant results within this paper are available on the Materials Cloud platform.68 See DOI: https://doi.org/10.24435/materialscloud:jy-tw.
APPENDIX A: THERMAL CONDUCTIVITY FOR DIFFUSIVE, MULTICOMPONENT SYSTEMS
The thermal conductivity, κ, is the proportionality coefficient between the energy flux and the (negative of the) temperature gradient in the absence of any convection. For a two-component system, like the SI materials studied in this work, characterized by the presence of one diffusive species,61,69 the Green–Kubo theory of linear response allows to extract κ from MD simulations as the zero frequency component of
where
is the power spectrum of the fluxes JA(t) and JB(t). The flux Je(t) is the total energy flux, here computed via the compute heat/flux command of LAMMPS, while is the convective flux of the diffusive species. Both Je(t) and Jd(t) have been computed in the barycentric reference frame. A multivariate technique72,73 for the cepstral analysis of time-series of the energy flux obtained from MD simulations has been employed in this work,70,71 allowing one to compute κ, efficiently and rigorously, for multicomponent superionic materials. As for the diffusivity, I checked any possible influence of the thermostat on the value of κ, by repeating several simulations in the NVE ensemble, previously equilibrated at the desired temperature. The results from NVT and NVE are fully compatible (see the supplementary material, Fig. S10).
The second term between square brackets in Eq. (A1) represents the contribution to heat flux due to convection. It must be removed from the first term, See, to correctly calculate κ, i.e., the coefficient of thermal conduction.64,69,74,75. Its effects on the thermal conductivity of CaF2 at different temperatures and system sizes are shown in Fig. 7: the full calculation, employing Eq. (A1) (solid lines and error bars), coincides with the single component calculation (shaded areas), that is, , only at low temperatures, where diffusion is negligible. Note that the departure from the multicomponent value occurs at lower temperature for small sizes, in agreement with overestimation of D in small boxes at low T.
Comparison between multicomponent (lines and error bars) and single-component (shaded areas) calculation of the thermal conductivity, κ, as a function of temperature and system size, for CaF2.
Comparison between multicomponent (lines and error bars) and single-component (shaded areas) calculation of the thermal conductivity, κ, as a function of temperature and system size, for CaF2.
APPENDIX B: ROLE OF ADDITIONAL F− EMPTY SITES ON DIFFUSION
I investigated the role of a higher number of accessible empty sites, to which a F− ion can jump, by considering a strongly defected lead fluoride system where, for each simulation, ℓ3 fluorine ions are randomly selected and removed from the ℓ × ℓ × ℓ supercell with perfect stoichiometry, ℓ being, as usual, the number of replicas of the conventional cell. In this way, the concentration of F− vacancies is set to 1/8, and it is independent of ℓ. Therefore, the simulation cells possess 4ℓ3 lead ions and 7ℓ3 fluorine ions. The selection of the F− ions to remove, and the generation of resulting defected cell, is performed with the code ATOMSK.76 A different seed for the random selection is used for each temperature and size. After a long initial equilibration (400 ps), a production run in the NVT ensemble is performed with LAMMPS. A uniform, neutralizing, background charge distribution is implicitly applied in the simulations employing Ewald’s method.
Although such a high concentration of F− vacancies is hardly attainable in actual experimental samples, this model serves as a test bench to validate the general picture whereby the presence of empty sites tends to reduce FSE, as extensively discussed in the main manuscript. In fact, as shown in Fig. 8, in these simulations, convergence in the Arrhenius plots is reached already for the ℓ = 2 cell (N = 88 atoms), in striking contrast to cells with perfect stoichiometry. Another interesting difference with the perfect stoichiometry case is that the activation energy below the transition to the superionic phase is significantly lower than that above the SI transition (whose value is itself quite close to the one for perfect stoichiometry): the presence of empty sites facilitates F− ion diffusion below the transition to the SI phase, while it becomes much less relevant above it, where the system is globally disordered.
Arrhenius plot of the diffusivity of F− ions at different sizes, for systems where a number ℓ3 of F− ions are randomly selected and removed from the stoichiometric system. The activation energies for the diffusion process below and above the transition temperature to the SI regime are extracted from the fits to Eq. (3) and reported in the legend.
Arrhenius plot of the diffusivity of F− ions at different sizes, for systems where a number ℓ3 of F− ions are randomly selected and removed from the stoichiometric system. The activation energies for the diffusion process below and above the transition temperature to the SI regime are extracted from the fits to Eq. (3) and reported in the legend.
REFERENCES
Notice that, even in the case of liquids, the specific functional dependence of the thermal conductivity on the system size seems to be qualitatively affected by the particular pressure and temperature conditions of the simulation.31
Below the SI transition, the transient hopping mechanism is dominated by vacancy motion, while, in the SI phase, it can be “attributed in roughly equal measure to vacancy and interstitial motion” (verbatim from Ref. 57).
This is also confirmed by the temperature, lower for smaller systems, at which the multicomponent analysis departs from the single-component one, which assumes no atomic diffusion (see Fig. 7).