3. Atomic structures

In the same way that the finite-element solvers of QTCAD store and solve for quantities of interest on a mesh, the atomistic solvers of QTCAD store and solve for quantities on an atomic structure, which are stored in Atoms objects. An atomic structure is a finite collection of atoms, each of which is defined by its chemical species and its position in space.

At the most basic level, an atomic structure may be constructed by tiling a crystallographic unit cell (see Unit cells) over a given volume of space. This procedure results in the atomic structure of a perfect crystal, which may be experimentally relevant for unstrained structures composed of a single material. More generally, it results in an atomic structure which serves as a starting point to generate more realistic structures. This procedure in described in Unit cell tiling.

In the context of semiconductor heterostructure simulations, the arrangement of atoms in an atomic structure notably defines, at the atomistic level:

• the chemical nature and geometry of heterostructure layers;

• disorder in the heterostructure layers, such as random alloying;

• rough interfaces between heterostructure layers;

• strain.

Rough interfaces are described in Random rough surfaces, random alloying is described in Alloys, and strain is described in Atomistic strain tensor.

For strained structures and/or heterostructures of materials with dissimilar lattice constants, the atomic coordinates obtained by the unit cell tiling procedure may be displaced according to a given strain tensor field, as described in Displacing atoms according to a strain tensor field. Alternatively, the atomic coordinates may be relaxed within the Keating valence force-field model (see Keating valence force-field model solver).

Finally, atomic structure are commonly stored in the .xyz file format, which is described in The .xyz file format.

Unit cell tiling

The most basic way to build an atomic structure in QTCAD is by tiling a unit cell so as to fill a given volume of space. This is done within the constructor method of the Atoms class and involves the following steps:

• the minimal and maximal \(x\), \(y\), and \(z\) coordinates of the box to be filled with atoms are specified; this is done through the x, y, and z arguments of the constructor method;

• the unit cell to be tiled is specified; this is done through the unit_cell argument of the constructor method;

• the crystal directions to be made parallel to the \(x\), \(y\), and \(z\) axes of the box are specified; this is done through the crystal_dir argument of the constructor method.

Algorithmically, the unit cell tiling procedure consists of the following steps:

• rotate the unit cell unit_cell so that the crystal directions specified in crystal_dir are parallel to the \(x\), \(y\), and \(z\) axes;

• tile the unit cell unit_cell over the volume of the box whose boundaries are defined in x, y, and z;

• eliminate any surface atom that has fewer than two nearest neighbor, where the notion of nearest neighbor is defined by the interatomic distance set through the argument vff_params_name of the constructor method.

This procedure produces the atomic structure of a perfect crystal lattice of a single material (defined through unit_cell). The next step in modeling heterostructures is to define the chemical composition and geometry of the heterostructure layers.

Note

Silicon oxide is generally an amorphous material. Amorphous materials cannot be directly modeled in the atomistic simulations performed in the qtcad.atoms module. Thus, the silicon oxide material is instead modeled within the virtual crystal approximation, wherein fictitious \(\mathrm{SiO}_{2}\) “atoms” are arranged on a diamond lattice with the same lattice constant as that of bulk silicon, following the procedure described in Ref. [KLP+11].

Alloys

An alloy is a mixture of two or more chemical species. It is characterized by the ratio of each chemical species in the alloy. For example, in the \(\mathrm{Si}_{x}\mathrm{Ge}_{1-x}\) alloy, a ratio \(0 \leq x \leq 1\) of all atoms are silicon atoms, while a ratio \(1-x\) of all atoms are germanium atoms. QTCAD can generate three types of alloys:

• ordered alloys;

• random alloys with fixed ratios;

• random alloys with stochastic ratios.

Below, we describe how alloy generation is performed. We consider alloys denoted by \(\mathrm{A}^{(1)}_{x_1}\mathrm{A}^{(2)}_{x_2}\ldots\), where \(\mathrm{A}^{(i)}\) is the \(i^{\mathrm{th}}\) chemical species in the alloy and \(x_i\) is its ratio.

Ordered alloys

Note that the conventional unit cell of a zincblende crystal structure contains eight atoms. Since this unit cell is repeated periodically in space to build atomic structures (see Unit cell tiling), each atom in an atomic structure is a periodic image of exactly one of the eight atoms in the conventional unit cell. In other words, an atomic structure can be partitioned into exactly eight sublattices, each of which corresponds to one of the eight atoms in the conventional unit cell. The sublattice indices of atoms are stored in the attribute atom_sublat of the Atoms objects. Ordered alloys are thus generated by assigning each of the eight sublattices to a different chemical species.

This alloying method is called “ordered” because the chemical species are arranged in a predictable way. See, for example, the \(\mathrm{Si}_{0.875}\mathrm{Ge}_{0.125}\) alloy shown below. Ordered alloys are typically difficult, if not impossible, to realize experimentally. However, in the context of simulations, they may be useful for testing purposes, such as for assessing the impact of random alloying on quantum dot characteristics. Note that this alloy generation method can only model alloys where, for each chemical species \(\mathrm{A}^{(i)}\) present in the alloy, the ratio \(x_i\) is a multiple of \(\frac{1}{8}\).

Fig. 3.4 Ordered \(\mathrm{Si}_{0.875}\mathrm{Ge}_{0.125}\) alloy produced by QTCAD. Blue spheres are silicon atoms and purple spheres are germanium atoms.

Random alloys with fixed ratios

The atomic structures modeled in QTCAD are finite in size and, as such, contain a finite number of atoms. To generate a random alloy with fixed ratios in a given region of an atomic structure containing \(N\) atoms, for each chemical species \(\mathrm{A}^{(i)}\) present in the alloy:

• \(x_i N\) atoms (up to rounding error, assuming \(x_i N\) is not an integer) are randomly selected from the \(N\) atoms in the region;

• the selected atoms are assigned to be of the chemical species \(\mathrm{A}^{(i)}\).

This random alloying method leads to a random distribution of the chemical species present in the alloy and guarantees that the ratio of each chemical species is equal to the desired ratio, up to rounding error, as for the \(\mathrm{Si}_{0.875}\mathrm{Ge}_{0.125}\) alloy shown below. In structures too small for averaging effects to play an important role, this type of random alloy may be difficult to realize experimentally. However, for simulation purposes, this alloy generation method may be useful to assess the impact of the deviations from the desired alloy ratios on quantum dot characteristics.

Fig. 3.5 Random \(\mathrm{Si}_{0.875}\mathrm{Ge}_{0.125}\) alloy with fixed ratios produced by QTCAD. Blue spheres are silicon atoms and purple spheres are germanium atoms.

Random alloys with stochastic ratios

To generate a random alloy with stochastic ratios in a given region of an atomic structure, each of its \(N\) atoms is assigned a probability \(x_1\) to be of chemical species \(\mathrm{A}^{(1)}\), a probability \(x_2\) to be of chemical species \(\mathrm{A}^{(2)}\), and so on. The alloy is then generated by randomly assigning each atom a chemical species according to the given probabilities.

In the limit of infinitely large structures, due to averaging effects, this random alloying method is equivalent to the random alloying method with fixed ratios. However, in small structures, this random alloying method may lead to an alloy composition that deviates from the desired ratios. For example, in the \(\mathrm{Si}_{0.875}\mathrm{Ge}_{0.125}\) alloy shown below, only \(4\) out of \(62\) atoms are germanium atoms, leading to a real germanium ratio that is only around half of the desired ratio. This alloy generation method is the most realistic of the three, as it captures the full extent of disorder and variability that alloys typically exhibit experimentally.

Fig. 3.6 Random \(\mathrm{Si}_{0.875}\mathrm{Ge}_{0.125}\) alloy with stochastic ratios produced by QTCAD. Blue spheres are silicon atoms and purple spheres are germanium atoms.

Defining a heterostructure layer

A heterostructure layer is characterized by its geometry and material composition. In QTCAD, a heterostructure layer may be defined using the new_layer method.

This method takes as input z_bot and z_top, which specify the bottom and top surfaces of the layer along the \(z\) axis. For planar/pristine heterointerfaces, z_top and z_bot may be given as numbers representing the \(z\) coordinate of the top and bottom interfaces of the layer. For rough heterointerfaces, z_top and z_bot may be given as RoughSurface objects, which are described in Random rough surfaces.

This method also takes as input atom_species and method, which specify the chemical species of the atoms in the layer and the alloying method, as described in Alloys.

Algorithmically, the chemical species of the atoms in a heterostructure are adjusted as follows upon calling the new_layer method:

• the coordinate along \(z\) of each atom is compared to the height along \(z\) of the top and bottom surfaces of the layer specified in z_top and z_bot;

• if the coordinate along \(z\) of an atom is between the top and bottom surfaces of the layer, then the atom is assigned to be of one of the chemical species specified in atom_species with the alloying method specified in method.

Note that the atomic positions of the atoms are not modified upon calling the new_layer method.

Atomistic strain tensor

The concept of strain tensor originally stems from continuum mechanics, wherein materials are treated as continuum media. In heterostructure quantum dots, wherein the strain “felt” by an indidual atom or a small group of atoms may greatly influence the physical properties of the dot, this continuum mechanics approach may not be appropriate.

In this section, we introduce the atomistic strain tensor, which is a generalization of the continuum mechanics notion of strain tensor. It is this atomistic strain tensor that is (implicitly) used in the Keating valence force field model (see Keating valence force-field model solver) as well as in the tight-binding model (see Tight-binding Schrödinger solver) to account for the effect of strain on the physics of an atomic structure. It may be computed using the get_atomistic_strain_tensor method.

The atomistic strain tensor, \(\boldsymbol{\epsilon}\), of a tetrahedrally-bonded atom is defined by the following equation [PKW+98],

(3.39)\[\begin{split}\begin{pmatrix} \epsilon_{xx} & \epsilon_{yx} & \epsilon_{zx}\\ \epsilon_{xy} & \epsilon_{yy} & \epsilon_{zy}\\ \epsilon_{xz} & \epsilon_{yz} & \epsilon_{zz} \end{pmatrix} = \begin{pmatrix} R_{12,x} & R_{23,x} & R_{34,x}\\ R_{12,y} & R_{23,y} & R_{34,y}\\ R_{12,z} & R_{23,z} & R_{34,z} \end{pmatrix} \begin{pmatrix} R^{0}_{12,x} & R^{0}_{23,x} & R^{0}_{34,x}\\ R^{0}_{12,y} & R^{0}_{23,y} & R^{0}_{34,y}\\ R^{0}_{12,z} & R^{0}_{23,z} & R^{0}_{34,z} \end{pmatrix}^{-1} - I\end{split}\]

To understand the parameters of this equation, consider that the atom sits at the center of a tetrahedron whose vertices are its four nearest-neighbouring atoms. The indices \(1\), \(2\), \(3\), and \(4\) refer to the four nearest-neighbouring atoms of the considered atom, which we index by \(0\), for convenience. The vector \(\mathbf{R}_{ij}\) points from atom \(i=1,2,3\) to atom \(j=2,3,4\). The vectors \(\mathbf{R}_{12}\), \(\mathbf{R}_{23}\), and \(\mathbf{R}_{34}\) are thus edges of the tetrahedron formed by the four nearest-neighbouring atoms of atom \(0\). In strained structure, this tetrahedron is not regular, but rather distorted.

The vector \(\mathbf{R}^{0}_{ij}\) are the corresponding edges of a “reference” tetrahedron. The strain tensor of Eq. (3.39) thus simply describes the distortion of the tetrahedron from its “reference” configuration.

Assuming the atoms \(1\), \(2\), \(3\), and \(4\) are all of the same chemical species, the “reference” tetrahedron defined by \(\mathbf{R}^{0}_{ij}\) is simply a regular tetrahedron in which the distance between atom \(0\) and atoms \(1, 2, 3, 4\) is set to the equilibrium bond length specified through the attribute vff_params_name of the Atoms object.

The scientific literature does not provide a consensus on the definition of the atomistic strain tensor in cases where the atoms \(1\), \(2\), \(3\), and \(4\) are not all of the same chemical species. In this case, we make the following assumption in QTCAD. The “reference” tetrahedron is instead set so that the atoms \(1\), \(2\), \(3\), and \(4\) lie along the same directions as the corners of the regular tetrahedron, but with bond lengths set to the equilibrium bond lengths of the corresponding chemical species, as specified through the attribute vff_params_name of the Atoms object.

The matrix \(I\) is simply the \(3 \times 3\) identity matrix.

We note that for atoms that are not tetrahedrally bonded, namely, surface atoms, the atomistic strain tensor is set to zero. In addition, it can be shown that the strain tensor resulting from Eq. (3.39) is independent of the ordering of the atoms \(1\), \(2\), \(3\), \(4\).

Finally, we note that, for small strain, the continuum concept of strain tensor is equivalent to the atomistic strain tensor [PKW+98].

Displacing atoms according to a strain tensor field

Strain in heterostructures may be induced by thermal misfit during the growth of heterostructure layers or by external forces. To capture such effects, it is possible to displace atoms in an atomic structure according to a given strain tensor field. This may be done through the strain_atomic_structure method.

In infinitesimal (continuum) strain theory, the strain tensor field is defined as [BS02]:

(3.40)\[\varepsilon_{\alpha\beta} \left(\mathbf{s}\right) = \frac{1}{2} \left[ \frac{\partial u_{\alpha}}{\partial x_{\beta}} \left(\mathbf{s}\right) + \frac{\partial u_{\beta}}{\partial x_{\alpha}} \left(\mathbf{s}\right) \right]\,,\]

where \(\mathbf{s}\) is a position of a point in space, the subscripts \(\alpha\) and \(\beta\) run over the Cartesian coordinates (\(x\), \(y\), and \(z\)), and \(\mathbf{u}\left(\mathbf{s}\right)\) is the displacement vector at point \(\mathbf{s}\).

Noting that \(\frac{\partial u_{\alpha}}{\partial x_{\beta}} = \frac{1}{2} \left( \frac{\partial u_{\alpha}}{\partial x_{\beta}} + \frac{\partial u_{\beta}}{\partial x_{\alpha}}\right) + \frac{1}{2} \left( \frac{\partial u_{\alpha}}{\partial x_{\beta}} - \frac{\partial u_{\beta}}{\partial x_{\alpha}}\right)\), it can be seen that \(\frac{\partial u_{\alpha}}{\partial x_{\beta}} = \frac{\partial u_{\beta}}{\partial x_{\alpha}}\) if:

(3.41)\[\frac{\partial u_{\alpha}}{\partial x_{\beta}} - \frac{\partial u_{\beta}}{\partial x_{\alpha}}=0\,,\]

namely, if the displacement vector field is irrotational, or in other words, curl-free. In this case, the strain tensor field may be expressed as:

(3.42)\[\varepsilon_{\alpha\beta} \left(\mathbf{s}\right) = \frac{\partial u_{\alpha}}{\partial x_{\beta}} \left(\mathbf{s}\right)\,.\]

The inverse problem to Eq. (3.40), namely the problem of finding a displacement vector field that corresponds to a given strain tensor field, is non-trivial. In QTCAD, this inverse problem is solved under the assumption of Eq. (3.41). In this case, the displacement vector field may then be obtained by integrating Eq. (3.42):

(3.43)\[\mathbf{u} \left(\mathbf{s}\right) = \int\limits_{\gamma\left(\mathbf{s}_0,\mathbf{s}\right)} \boldsymbol{\varepsilon} \left(\mathbf{t}\right) \cdot d\mathbf{t}\,,\]

where \(\mathbf{s}_0\) is a reference point, and \(\gamma\left(\mathbf{s}_0,\mathbf{s}\right)\) is any smooth path starting at \(\mathbf{s}_0\) and ending at \(\mathbf{s}\). From there, the atoms in an atomic structure may be displaced according to the dispacement vector field. Specifically, the position \(\mathbf{r}_i\) of an atom labeled by \(i\) in the atomic structure is updated as:

(3.44)\[\mathbf{r}_i \rightarrow \mathbf{r}_i + \mathbf{u} \left(\mathbf{r}_i\right)\,.\]

We note that the strained atomic structure is translated so that its center of mass matches that of the unstrained atomic structure. In addition, it is rotated so as to minimize its root-mean-square deviation from the unstrained atomic structure using the Kabsch algorithm [Kab76]. This is done to ensure that the relaxed atomic structure maintains the same position in space and crystal orientation as the unstrained atomic structure. In addition, this ensures that the resulting strained atomic structure is independent of the choice of the reference point \(\mathbf{s}_0\).

The .xyz file format

The .xyz file format is commonly used in computational chemistry, material physics, and device engineering to store atomic structures. It is a simple text file format that contains information about the chemical species and positions of atoms in an atomic structure. While no standard exists for the .xyz file format, QTCAD uses one of the most common conventions, which is as follows:

• the first line of an .xyz file contains the number of atoms in the atomic structure; call this number \(N\);

• the second line is a comment line, which may be empty; it typically contains a name and/or description of the atomic structure;

• the following \(N\) lines each contain the following four columns of data, which are tab-separated:

  • the first column contains the chemical species (e.g., Si, Ge, etc.) of an atom;

  • the second, third, and fourth columns contain the \(x\), \(y\), and \(z\) coordinates of the atom in units of angstroms.

For example, the following .xyz file describes the conventional unit cell of bulk silicon.

8
Conventional unit cell of bulk silicon
Si    0.0000000000000 0.0000000000000 0.0000000000000
Si    2.7154977091011 2.7154977091011 0.0000000000000
Si    2.7154977091011 0.0000000000000 2.7154977091011
Si    0.0000000000000 2.7154977091011 2.7154977091011
Si    1.3577488545506 1.3577488545506 1.3577488545506
Si    4.0732465636516 4.0732465636516 1.3577488545506
Si    4.0732465636516 1.3577488545506 4.0732465636516
Si    1.3577488545506 4.0732465636516 4.0732465636516

Atomic structures may be visualized using various softwares, for example VESTA. Below, we show the silicon unit cell as visualized in VESTA.

Fig. 3.7 Conventional unit cell of bulk silicon visualized in VESTA.

Note that, for practical simulations, the atomic structures of semiconductor heterostructure quantum dots may contain up to millions of atoms; such large structures may be slow to process by atomic structure visualization softwares.

The .xyz file associated with an atomic structure may be generated by calling the write_xyz method.