Polytypism in the ground state structure of the Lennard-Jonesium

TL;DR

This paper systematically investigates the stability of various polytypic structures in Lennard-Jones materials, revealing that the ground state can differ from traditional fcc or hcp structures depending on pressure and potential truncation.

Contribution

It provides a comprehensive analysis of nineteen structures, highlighting the influence of pressure, cutoff, and potential parameters on the ground state configurations.

Findings

Polytypic sequences can be more stable than fcc or hcp under certain conditions.

Ground state structures depend strongly on potential truncation methods.

Different polytypes emerge as stable configurations based on pressure and cutoff choices.

Abstract

We present a systematic study of the stability of nineteen different periodic structures using the finite range Lennard-Jones potential model discussing the effects of pressure, potential truncation, cutoff distance and Lennard-Jones exponents. The structures considered are the hexagonal close packed (hcp), face centred cubic (fcc) and seventeen other polytype stacking sequences, such as dhcp and $9R$. We found that at certain pressure and cutoff distance values, neither fcc nor hcp is the ground state structure as previously documented, but different polytypic sequences. This behaviour shows a strong dependence on the way the tail of the potential is truncated.

Tables

Table 1. Table 1: Different stacking variants studied, listed with both their hc and ABC notation. Alternative names and Ramsdell notations are shown for specific stacking sequences in parenthesis.

stacking	min. number	physical stacking of the layers
variants	of layers
c (fcc, $3C$)	3	${\left[ABC\right]}_n$
h (hcp, $2H$)	2	${\left[AB\right]}_n$
hc (dhcp,$4H$)	4	${\left[ABCB\right]}_n$
hcc (thcp,$6H_1$)	6	${\left[ABCACB\right]}_n$
hccc	8	${\left[ABCABACB\right]}_n$
hcccc	10	${\left[ABCABCBACB\right]}_n$
hccccc	12	${\left[ABCABCACBACB\right]}_n$
hhc ($9R$)	9	${\left[ABACACBCB\right]}_n$
hhcc	12	${\left[ABACBCBACACB\right]}_n$
hhccc	5	${\left[ABACB\right]}_n$
hhcccc	18	${\left[ABACBACACBACBCBACB\right]}_n$
hhhc	8	${\left[ABABCBCB\right]}_n$
hhhcc	10	${\left[ABABCACACB\right]}_n$
hhhccc	12	${\left[ABABCABABACB\right]}_n$
hhhhc	15	${\left[ABABACACACBCBCB\right]}_n$
hhhhcc	18	${\left[ABABACBCBCBACACACB\right]}_n$
hhhhhc	12	${\left[ABABABCBCBCB\right]}_n$
hchcc ($15R$)	15	${\left[ABCBACABACBCACB\right]}_n$
hchhc	10	${\left[ABCBCACBCB\right]}_n$

Equations

$$U_{\mathrm{LJ}}(r)=\frac{p}{(p-q)(\frac{q}{p})^{q/(p-q)}}\epsilon\bigg{[}\Big{(}\frac{\sigma}{r}\Big{)}^{p}-\Big{(}\frac{\sigma}{r}\Big{)}^{q}\bigg{]}.$$

$$U_{\mathrm{LJ}}(r)=4\epsilon\bigg{[}\Big{(}\frac{\sigma}{r}\Big{)}^{12}-\Big{(}\frac{\sigma}{r}\Big{)}^{6}\bigg{]},$$

$$U_{\mathrm{LJ}-C^0}(r)=\{\begin{array}{ll}{U}_{\mathrm{LJ}}(r)-U_{\mathrm{LJ}}(r_c)& r\le r_c\\ 0& r>r_c.\end{array}$$

$$U_{\mathrm{LJ}-C^1}(r)=\{\begin{array}{ll}\begin{array}{rl}{U}_{\mathrm{LJ}}(r)& -U_{\mathrm{LJ}}(r_c)-\\ & -(r-r_c)U_{\mathrm{LJ}}^{\prime }(r_c)\end{array}& r\le r_c\\ 0& r>r_c.\end{array}$$

$$U_{\mathrm{LJ}-C^2}(r)=\{\begin{array}{ll}\begin{array}{rl}{U}_{\mathrm{LJ}}(r)-U_{\mathrm{LJ}}(r_c)& \\ -(r-r_c)U_{\mathrm{LJ}}^{\prime }(r_c)& \\ -\frac{1}{2}(r-r_c{)}^2{U}_{\mathrm{LJ}}^{\prime \prime }(r_c)& \end{array}& r\le r_c\\ 0& r>r_c.\end{array}$$

$$U_{\mathrm{LJ}-C^1-\mathrm{smoothed}}(r)=f_s(r)U_{\mathrm{LJ}-C^1}(r).$$

$$f_s(r)\{\begin{array}{ll}1& r\le (r_c-r_s)\\ 1-\frac{e^{(-1/x)}}{e^{(-1/x)}+e^{(-1/(1-x))}}& (r_c-r_s)<r\le r_c\\ 0& r>r_c,\end{array}$$

$$U_{\mathrm{pl}}(r)=4\epsilon\Big{(}\frac{\sigma}{r}\Big{)}^{12},$$

$$U_{\mathrm{M}}(r)=D_{e}\Big{(}e^{-2a(r-r_{e})}-2e^{-a(r-r_{e})}\Big{)},$$

$$H=-\frac{1}{2}{J}_0\sum _{ij{j}^{\prime }}{S}_{i,j}{S}_{i,j^{\prime }}-J_1\sum _{ij}{S}_{i,j}{S}_{i+1,j}-J_2\sum _{ij}{S}_{i,j}{S}_{i+2,j}$$

Full text

I Introduction

Polytypism is a special form of polymorphism, occurring in layered materials, in which the polymorphs are derived simply by varying the way in which the layers are arranged relative to each other. This means that the various stacking arrangements do not affect the chemistry of the phase as a whole, but some of the physical properties (e.g. density, Young modulus, band gap or electron mobility) can be significantly different. A large variety of materials have several different stable polytype phases Trigunayat (1991), one of the most extensively studied being SiC Jepps and Page (1983); Nakashima and Hangyo (1991). SiC has more than 200 identified polytypes, a few being more favoured in applications than the rest due to their superior electronic properties. Many materials with similar structural properties also form polytypes, such as metal sulphides and halogenides, e.g. ZnS Mardix et al. (1969); Boutaiba et al. (2014) and CdI2 Katahama et al. (1984). The physical properties of such materials can be tuned by changing the stacking sequences, e.g. in ZnO Zagorac et al. (2015). Polytypism also occurs in the case of diamond. The common cubic form of diamond has a hexagonal polytype called Lonsdaleite, which is suggested to be a complex mixture of different stacking sequences Wen et al. (2008); Salzmann et al. (2015). Similarly, hexagonal (Ih) and cubic (Ic) ice are also polytypic structures. Some elements are also known to form polytype phases, such as lanthanum, which exists in the dhcp form Bağci et al. (2010), samarium and lithium having the $9R$ stacking sequence as the ground state structure Ellinger and Zachariasen (1953); Overhauser (1984), erbium which is stable in both dhcp and 9R stacking sequences at different pressures Samudrala et al. (2011), and bismuth, long suspected to exist in several polytypic forms Shu et al. (2016). It has been speculated that the transformation from fcc to hcp structure with increasing pressure might occur through a series of different stacking fault structures, e.g. as in the case of noble gases xenon and krypton, suggested by some experimental results Jephcoat et al. (1987); Errandonea et al. (2002), or in the case of iron at high pressure and high temperature Yoo et al. (1996); Mikhaylushkin et al. (2007). Finally, if a wider definition of polytypism is used such that structurally compatible modules are also considered, a range of minerals which include the pyroxenes, perovskites, spinelloids, chlorites and oxides form polytypic structures as well Price and Yeomans (1984).

In order to model the polytypic behaviour, the axial next-nearest-neighbour Ising (ANNNI) model, was used in the 1980s Price and Yeomans (1984); J Yeomans and G D Price (1986); J Yeomans (1988). (An overview of the ANNNI and A3NNI ground state structures are given in the Appendix.) However, the two possible layer types in ANNNI, usually marked by $\downarrow$ and $\uparrow$ are interchangeable, a phase is only defined by the number of consecutive layers of the same orientation, but not the orientation of the layer itself, thus phases $\downarrow\downarrow\uparrow$ and $\uparrow\uparrow\downarrow$ are identical. In contrast, close packed stacking structures are built up by layers in three different possible positions, usually denoted by A, B and C, forming either a hexagonal (in ABA stacking) or a cubic (in ABC stacking) layer. These two are not identical nor interchangeable, meaning that the ability of ANNNI to describe the behaviour of close packed materials is limited.

One of the most widely used models to study close-packed materials is the Lennard-Jones pair potential. It has been long known that its low temperature dominant structures are the hexagonal close packed (hcp) and the face centred cubic (fcc) Kihara and Koba (1952). Interestingly, although other structures (bcc, simple cubic, diamond) have been studied Travesset (2014), to the best of our knowledge no other polytype sequences have ever been investigated from the point of view of phase stability, only as stacking faults in relation to crystal growth defects or nucleation de Souza and Wales (2016). It is also notable that the customary finite range truncation of the potential has a significant effect on the liquid-vapour equilibrium Smit (1992); Panagiotopoulos (1994); Shi and Johnson (2001) and on the melting temperature Mastny and de Pablo (2007); Ahmed and Sadusa (2010), yet it is rarely mentioned and almost never taken into account in the discussion of the low temperature solid phases, causing an apparent inconsistency in the literature with regard to the lowest energy structure: some works refer to the fcc de Souza and Wales (2016); van der Hoef (2000), others to hcp van de Waal (1991); Travesset (2014) as the global minimum of the Lennard-Jones model. An exception is an article by Jackson et. al Jackson et al. (2002) showing that the ground state can be either fcc or hcp depending on the cutoff distance and method of truncation.

Our aim in this work is to provide a systematic study of the ground state structure of the Lennard-Jones potential considering different polytypic stacking sequences, and fill the gap in the literature regarding its dependence on pressure, potential truncation and potential parameters. The rich diversity of structures we find serve as a reminder that complex material behaviour can result from comparatively simple models, and that implementation details can have a strong effect on phase stability when modelling materials.

II Computational details

A generalised form of the Lennard-Jones potential can be given by

[TABLE]

The values $p=12$ and $q=6$ are most common, in which case one obtains

[TABLE]

where $\epsilon$ is the depth of the potential well, $\sigma$ is the size of the repulsive core and $r$ is the distance between two particles. This potential is usually truncated at a cut-off distance, $r_{c}$, and to avoid the discontinuity at this point, the potential can be shifted. The energy-shifted LJ potential is a continuous ($C^{0}$) function,

[TABLE]

In order to obtain continuous forces at the cut-off distance, and thus make the potential function differentiable, it can be force-shifted, which leads to

[TABLE]

In order to make also the second derivatives continuous at $r_{c}$ (i.e. create a $C^{2}$ function), the potential can be further shifted by a third term,

[TABLE]

Alternatively, a sigmoidal shaped function, $f_{s}(r)$, can be used to ”smooth out” the force shifted potential within a distance $r_{s}$ of the cutoff,

[TABLE]

For $f_{s}(r)$ we used an infinitely differentiable ($C^{\infty}$) function,

[TABLE]

where $x=[r-(r_{c}-r_{s})]/r_{s}$).

Nineteen different stacking sequences were considered (shown in Table 1), all variations up to five stacking layers, and some of the possible six layer arrangements. The notation used to mark the sequences is that layers with fcc surroundings (the two neighbouring layers occupy different positions) is marked “c” as cubic, while the layers with hcp surrounding (sandwiched between two layers occupying the same position) is marked “h” as hexagonal. An example structure, with associated notation, is shown in Figure 2.

The geometry optimisations were performed with the QUIP package qui , with the conjugate gradient method and double-checked with the steepest descent method for several cases. The minimisation tolerance was set to $10^{-6}\epsilon/\sigma$ for the norm of the forces, which corresponds to $10^{-15}\epsilon/\mathrm{atom}$ accuracy in the energy calculation. Hydrostatic pressure was applied. The minimisations were started from configurations where atoms were placed $1.0\sigma$ distance from each other, and during the minimisation the atomic positions and all the lattice parameters were allowed to relax. The calculations were done with different shifted and smoothed versions of the potential to compare their effects within the truncation range $2.0\sigma\leq r_{c}\leq 6.0\sigma$, in $0.05\sigma$ intervals.

III Results: ground state phase diagrams of the Lennard-Jones potential

During the minimisation process configurations retained their stacking order, and the atoms forming the stacking plane also stayed perfectly in the plane. Lattices remained orthorhombic, but the lattice height corresponding to the stacking direction changed with respect to the other two, as expected. To be able to identify the ground state structures and draw the phase diagrams, the enthalpies of the minimised configurations were calculated and compared. These enthalpy curves were individually checked and more calculations were performed with a finer pressure scale whenever it was necessary, thus we believe that no phases have been missed. A set of example enthalpy curves can be seen on Figure 3.

Truncation distance vs. pressure phase diagrams of the Lennard-Jones type potentials with different cutoff schemes are shown on Figures 4, 5 and 6. The coloured regions show the series of phases found to be the most stable at a given truncation length and in a given pressure range. Dark grey colour corresponds to the fcc and light grey to hcp structures, with other colours representing different stacking variants (yellow and red shades represent stackings with a single “h” layer in the repeated subunit, green shades are polytypes with two consecutive “h” layers, while blues and purples correspond to three and four consecutive “h” layers, respectively).

It is clear from all the phase diagrams, that the hcp structure tends to be the most favourable stacking variant at lower pressure values, while for every value of the cutoff there is a pressure above which the fcc is the most stable polytype. In order to see whether one of the studied polytypes becomes the ground state again at even higher pressures, the structures were minimised up to $p=6\times 10^{5}\epsilon/\sigma^{3}$ for a few randomly chosen cutoffs: fcc remained the lowest enthalpy structure in every case.

However, the most striking result is that at the boundary between the ground state regions of fcc and hcp structures, several other stacking variants are found to be more stable. This means that, in contrast with the common belief, the (truncated) Lennard-Jones potential exhibits a wide range of different global minima depending on the fine details of the potential.

The boundary between the fcc and hcp regions appear to have a “wave”-like pattern for all the potential function variants we used. The shape of these waves reflect how the distance of atomic shells decreases as the density increases with increasing pressure. (The white dashed graph in Figure 4 represents the curve along which the number of atoms within the cutoff radius jumps from 177 to 201 in the fcc crystal.) As different polytypes have different numbers of neighbours in each shell, their relative energy will be different depending on which shells lie within the cutoff radius. As the pressure increases the atoms get closer, farther shells appear within the smaller cutoff radii, causing the phase boundary to be shifted towards smaller cutoffs.

At small cutoff values, fewer polytypes appear along the fcc-hcp boundary and these remain the same as the pressure increases. As the number of shells are increased using larger cutoffs, this is no longer true, due to the fact that the distance between the layers of fcc and hcp can be different depending on the pressure, thus the neighbour shells will no longer be isotropic. Finally, with increasing cutoff, the energy contribution of the outmost shells gets smaller, the “waves” gradually flatten out.

Figure 4 shows the phase diagram of the simple energy shifted Lennard-Jones function (see eq. 3). At smaller cutoff values, only one phase is found to be stable other than fcc and hcp, the $\langle\mathrm{hc}\rangle$ phase. As the cutoff increases, first the $\langle\mathrm{hhc}\rangle$ and $\langle\mathrm{hhcc}\rangle$ phases appear, than other sequences with longer repeated subunits as well.

Applying additional shifts to the potential, such as force-shift and second-derivative shift, the “wave” like pattern of the phase diagram becomes significantly less pronounced (see the first two panels of Figure 5), but the order in which the more complex polytypes appear on the phase diagram remains similar. For example $\langle\mathrm{hhc}\rangle$ and $\langle\mathrm{hhcc}\rangle$ phases appear on the second “wave”, the same two and $\langle\mathrm{hc}\rangle$ on the two sides of the third “wave” and then $\langle\mathrm{hcccc}\rangle$ first appears on the tip of the fourth “wave” in all three phase diagrams. Although we are unable to offer a rigorous explanation for the flattening trend of the “waves”, we speculate that it is due to the fact that a non-smooth cut-off mechanism leads to large variations in energy as new neighbour shells cross the interaction range. Since different polytypes have different neighbour shells, significant changes in energy differences can therefore occur.

Using a smoothing function to obtain a completely smooth potential, $U_{{\rm LJ}-C^{\infty}}$ has seemingly the opposite effect, while the width of the stability region of the $\langle\mathrm{hc}\rangle$ phase becomes significantly narrower, especially at shorter cutoffs, the magnitude of the “waves” increases (see Figure 5). The explanation for this effect is that $U_{\rm LJ-C^{\infty}}$ is only qualitatively smooth but not quantitively. Indeed we observe in Figure 1 that its second derivative has rapid variations in the interval where the cut-off is applied. It also has to be noted that the exact effect of the smoothing will also depend on the widths of the smoothing region.

In order to study the effect of the Lennard-Jones exponents, thus the shape of the pair potential on the ground state phase diagram, we repeated our calculations on the force shifted potential with the following different $p$ and $q$ exponents: 12 and 8, 12 and 4, 14 and 6, 10 and 6. The phase diagrams are shown in Figure 6. It is clear from these figures that changing the exponents does not notably change the order in which the polytypic phases appear to be stable, thus the same stacking variants appear at the same cutoff values, but the corresponding pressure of the phase transitions are significantly different. As a general rule, if either of the exponents are increased, the pressure above which the fcc phase is the most stable increases as well.

IV Other interatomic potentials

IV.1 Power law potential

A simple repulsive power law potential with the exponent set to 12,

[TABLE]

was also tested to see whether polytypic phases are stable in this case. A force-shift was applied here as well. The results indicate that the ground state structure is fcc at every pressure and cutoff studied.

IV.2 Morse potential

In order to see how another simple pair potential with an attractive term behave we also tested the Morse potential,

[TABLE]

with parameters $D_{e}=1.0$, $r_{e}=1.0$ controlling the depth and location of the minimum, respectively, and $a=4.0$. (We chose the value for parameter $a$ so that the pair potential is similar in shape to the Lennard-Jones potential.) The potential was force shifted. The Morse phase diagram (see Figure 7) shows some stacking fault structures too: $\langle\mathrm{hc}\rangle$, $\langle\mathrm{hhc}\rangle$ and $\langle\mathrm{hhcc}\rangle$, similarly to the LJ potential, but only at small cutoffs. Above $r_{cut}=3.6r_{e}$ the fcc is the only stable phase. This fast decay of the “waves” can be explained by the exponential decay of the potential function, but the fact that polytypes other than fcc and hcp are found to be ground states also in case of the Morse potential indicates that the results we obtained for Lennard-Jones potentials in the previous sections might be generic for truncated pair potentials.

V Conclusion

We have systematically studied the global minimum structure of the bulk Lennard-Jones model as a function of pressure and the details of potential truncation. Our results demonstrate that its ground state structures are far more complex than previously reported, the stable phases including not only fcc and hcp but a wide range of more complex stacking sequences. Most notably we obtained $\langle\mathrm{hc}\rangle$ and $\langle\mathrm{hhc}\rangle$ phases, the two polytypes most often observed in real materials (as dhcp and 9R) other then hcp and fcc. This suggests that well-known pair potentials might be useful models of polytypism and can help us to understand and predict polytypic behaviour.

The relative stability of polytypes was found to be especially sensitive to the degree of smoothness of the potential around the cutoff. This shows that the effect of truncation and the way the derivative of the potential is treated should not be underestimated when using pair potentials. Further work is still needed, however, to obtain a clear theoretical explanation for the polytypism that we observed, e.g. to construct potentials with prescribed polytypes, as well as to confirm analogous effects in case of more complex model systems.

Acknowledgements.

LBP acknowledges support from the Royal Society through a Dorothy Hodgkin Research Fellowship. CO was supported by ERC Starting Grant 335120. ABP was supported by a Leverhulme Early Career Fellowship and the Isaac Newton Trust until 2016. CJP is also supported by the Royal Society through a Royal Society Wolfson Research Merit award.

Appendix A Stable structures of the ANNNI and A3NNI models

One of the simplest nontrivial models exhibiting periodically ordered phases, is the axial next-nearest-neighbour Ising model, ANNNI Elliott (1961); Fisher and Selke (1980), which is a variant of the Ising model with a two-state spin on each lattice site. Interactions are between nearest neighbours, together with a second-neighbour interaction along one lattice direction (this is the axial direction, $z$). The model is defined by the Hamiltonian

[TABLE]

where $S_{ij}$ is the two-state spin on each lattice site, and $i$ denotes the layers perpendicular to the axial direction and $j$ and $j^{\prime}$ are nearest neighbour spins within the layer.

The ANNNI model is considered a prototype for polytypism Price and Yeomans (1984); J Yeomans and G D Price (1986), since its phase diagram contains sequences of long-wavelength-modulated phases, hence showing that short-range competing interactions are sufficient to stabilise long periodic structures as ground states.

Within ANNNI, the layers building up the polytypes can be characterised by two signs, $\downarrow$ and $\uparrow$. In these models $\cdots\uparrow\uparrow\uparrow\cdots$ and $\cdots\downarrow\downarrow\downarrow\cdots$ are identical, and often simply marked as $\langle\infty\rangle$ (this is called the Zhdanov notation where the numbers in the brackets show the band widths, i.e. the number of layers with the same spin, e.g. $\langle 2\rangle$=(2,2) means $\cdots\uparrow\uparrow\downarrow\downarrow\cdots$).

The ground state phase diagram of ANNNI can be easily determined  Price and Yeomans (1984), see Figure 8. There are three main stable phases at 0 K, $\langle 1\rangle$, $\langle 2\rangle$ and $\langle\infty\rangle$, but the two dashed lines mark regions where the ground state is highly degenerate: along the line between $\langle 1\rangle$ and $\langle 2\rangle$ all phases containing only 1 and 2 bands have equal energy (e.g. $\langle 12\rangle$, $\langle 12^{2}\rangle$,…etc.), and along the line between $\langle 2\rangle$ and $\langle\infty\rangle$ all phases which contain no 1-bands have the same energy (e.g. $\langle 23\rangle$, $\langle 2^{2}4\rangle$, etc.). Note, that the boundary between $\langle 1\rangle$ and $\langle\infty\rangle$ is not degenerate! This means in particular that there are several phases missing from this phase diagram, e.g. $\langle 13\rangle$, $\langle 14\rangle$, and so forth, are not ground states at any value of $J_{1}$ or $J_{2}$.

In order to see how longer range interactions effect the stability of phases, the third neighbour Ising model, A3NNI has been studied too and discussed as a model for polytypism J Yeomans (1988); Muraoka et al. (1998); Salje et al. (1987). If the third neighbour interactions are also taken into account, the ground state phase diagram becomes more complicated and two additional structures appear as possible ground states phases compared to ANNNI, the $\langle 12\rangle$ and $\langle 3\rangle$ with several additional phases along the multiphase lines J Yeomans (1988); Salje et al. (1987); see Figure 9.

Many of the polytypic structures found in different materials resemble the phases seen in the ANNNI and A3NNI models. For example, PbI2 has a reversible phase transition between the phases $\langle\infty\rangle$ and $\langle 1\rangle$, with the phases $\langle 2\rangle$ and $\langle 12\rangle$ being observed under different growth conditions only Salje et al. (1987). Same is true for ZnS and AgI. Phase transitions between $\langle 1\rangle-\langle 2\rangle-\langle\infty\rangle$ are found in MgSiO3, and $\langle 1\rangle-\langle\infty\rangle-\langle 3\rangle-\langle 23\rangle$ are seen in case of SiC Jepps and Page (1983). However, there are also polytypic structures seen in experiments, e.g. $\langle 13\rangle$ in spinelloids, which do not occur in the ANNNI model.

At temperatures slightly above 0 K, the phase diagram remains similar but at the vicinity of the multiphase lines sequences of $\langle 12^{k}\rangle$ and $\langle 2^{k}3\rangle$ appear, and even more new phases at higher temperatures (though the proportion of disordered layers increase too) Fisher and Selke (1980).