5. Atomistic disorder in a Si spherical quantum dot in a Ge matrix
5.1. Requirements
5.1.1. Software components
QTCAD
Paraview
5.1.2. Python script
qtcad/examples/tutorials/spherical_dot.py
5.1.3. References
5.2. Briefing
In previous tutorials of the atoms package, we saw how the
granularity of matter at the atomic scale impacts the properties of quantum
dots (QDs). Specifically, we investigated the impact of atomistic disorder in
the forms of (1) random alloy disorder in Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting and (2) interface
roughness in Computing g-tensors in an FD-SOI device.
In this tutorial, we investigate the impact of two other forms of atomistic
disorder on QD properties: point charge defects and alloy concentration
gradients, both of which are known to occur in real devices. To illustrate
these effects, we consider a spherical Si QD embedded in a Ge matrix. The
geometry of this device, unlike those of the devices considered in the previous
tutorials, is not defined by layers grown along a specific crystallographic
direction. Instead, the QD has a more complex 3D geometry, which requires
features of the atomic structure builder of the
Atoms class not seen in previous tutorials.
5.3. Header
We start by importing the required classes and functions.
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh
from qtcad.device.materials import Si
from qtcad.atoms import Atoms, SubAtoms
from qtcad.atoms.unit_cell import UnitCellZincblende
from qtcad.atoms import Atoms
from qtcad.atoms.keating import Solver as KeatingSolver
from qtcad.atoms.keating import SolverParams as KeatingSolverParams
from qtcad.atoms.schrodinger import Solver as SchrodingerSolver
from qtcad.atoms.schrodinger import SolverParams as SchrodingerSolverParams
from qtcad.atoms.analysis import save_vtu, analyze_dot
All of these imports were seen in previous tutorials. For reference, see Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting.
5.4. A spherical dot without disorder
5.4.1. Defining a complex 3D geometry
We start by defining the geometry parameters of our device. It consists of a Si QD with a diameter of \(4\,\mathrm{nm}\) embedded in a cubic Ge matrix. We assume that this cubic Ge matrix has a side length of \(8\,\mathrm{nm}\).
# Geometry of the spherical Si quantum dot in a Ge matrix
diameter = 4e-9
box_size = diameter * 2
x = np.array([-box_size/2, box_size/2])
y = x.copy()
z = x.copy()
Here, the arrays x, y, and z define the bounding box of our atomic
structure.
Next, we define the unit cell that will be tiled to construct the atomic structure. Noting that the Ge matrix completely surrounds and is physically larger than the Si QD, we assume that crystallographic order in the atomic structure is primarily determined by the geometry of bulk Ge. Thus, we take the unit cell to be that of bulk Ge.
# Unit cell of bulk Ge
lattice_constant = 5.647e-10
unit_cell = UnitCellZincblende(lattice_constant=lattice_constant,
atom_species_ws=np.array(["Ge", "Ge"]))
Here, \(5.647\,\mathrm{\mathring{A}}\) is the lattice constant of bulk Ge.
Since we wish to model a Si QD that is spherical in shape, the
new_layer method seen in
Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting is not sufficient to define the atomic structure. Indeed,
that method only allows defining layers that are planar and normal to the
\(z\) axis. Instead, we may instantiate the atomic structure
from a mesh and harness the mesh’s physical groups to define regions of the
atomic structure.
# Construct atomic structure from mesh
script_dir = Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes" / "spherical_dot.msh"
mesh = Mesh(scaling_factor=diameter, filename=str(path_mesh))
atoms_from_mesh = Atoms(mesh=mesh, unit_cell=unit_cell)
atoms_from_mesh.new_region(atom_species=np.array(["Si"]*8), tag="ball")
path_vtu = script_dir / "output" / "spherical_dot_from_mesh.vtu"
save_vtu(atoms=atoms_from_mesh, out_dict={}, path=path_vtu)
Here, the mesh stored in the file spherical_dot.msh consists of a cube of
side length \(2\) with a spherical region of radius \(1\) at its
center. The cube and spherical regions are assigned the physical tags
"cube_shell" and "ball", respectively. The mesh is scaled by the input
scaling_factor=diameter to obtain the desired physical dimensions.
When instantiating the Atoms`<qtcad.atoms.atoms.Atoms> object from the
``mesh` input, the atomic structure is first constructed by tiling the
unit cell unit_cell over the volume of the mesh mesh; see the
constructor method of the
Atoms class for details. We then define the
Si QD. Similar to the
new_layer method, the
new_region method takes as input
an array of atomic species to fill the new region. Here, we specify that each
of the eight atoms of the conventional unit cell that was tiled to fill the
atomic structure should be replaced with Si atoms, i.e.
atom_species=np.array(["Si"]*8). The geometry of the new region is defined
by the input tag, which is the physical tag of the volume of the mesh that
defines the region. Here, we set tag="ball", which corresponds to the
spherical region at the center of the mesh that defines the Si QD. The resulting
atomic structure is then saved to a VTU file for visualization in Paraview
using the save_vtu function.
Fig. 5.4.1 Atomic structure of the spherical Si QD in a Ge matrix constructed from a mesh. A slice through the center of the QD and normal of the \(x\) axis is shown. Si (Ge) atoms are represented as blue (red) spheres.
We may also define complex 3D geometries without using a mesh, as we will see below.
To model the Ge matrix as extending beyond the boundaries of the simulated atomic structure, we apply periodic boundary conditions (PBCs) along all three spatial directions.
# Axes along which periodic boundary conditions are applied for the atomic
# structure relaxation
pbc_directions = ["x", "y", "z"]
Far from the QD–matrix interface, these PBCs effectively impose the crystal geometry to be that of bulk Ge.
Next, we create the atomic structure through the
Atoms class. First, the unit cell of bulk Ge
is tiled to fill the bounding box defined by the x, y, and z arrays.
Then, a spherical region at the center of the atomic structure is defined to be
filled with Si atoms instead of Ge atoms, thereby defining the Si QD. This is
done through the new_region method,
again, but this time without using a mesh.
# Create atomic structure, defining the Si spherical dot via the `is_inside`
# argument of the `new_region` method
atoms = Atoms(x=x, y=y, z=z, unit_cell=unit_cell,
pbc_directions=pbc_directions)
def is_inside_dot(x, y, z):
r = np.sqrt(x**2 + y**2 + z**2)
return r <= diameter / 2
atoms.new_region(atom_species=np.array(["Si"]*8), is_inside=is_inside_dot)
The geometry of the new region is instead defined
by the is_inside argument, which takes as input a function that returns
True for points inside the new region and False for points outside the new
region. Here, we define the function is_inside_dot to return True only for
points inside a sphere of \(2\,\mathrm{nm}\) radius centered at the origin.
This results in the same atomic structure as before, except that PBCs are now
applied along all three spatial directions.
We note that when considering PBCs, the atomic structure stored in the
Atoms object represents only the
fundamental domain of the periodic atomic structure. Only the
atoms contained within the bounding box defined by the x, y, and z
arrays are stored. In more technical terms, the atomic structure thus consists
of a single supercell of the periodic atomic structure, namely a single unit
cell of the superlattice defined by the PBCs.
5.4.2. Atomic structure relaxation
Next, we relax the atomic structure using the Keating valence force-field model. This procedure is identical to that seen in Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting.
# Relax atomic structure with periodic boundary conditions
keating_solver_params = KeatingSolverParams()
keating_solver_params.verbose = True
keating_solver = KeatingSolver(atoms=atoms,
solver_params=keating_solver_params)
keating_solver.solve()
Importantly, we note that the PBCs defined earlier are considered during the relaxation procedure. Thus, atoms at the boundaries of the simulated atomic structure interact with their periodic images, effectively simulating an atomic structure that extends beyond the simulated bounding box.
Through this procedure, we obtain an atomic structure whose neighborhood around the Si QD looks as shown below.
Fig. 5.4.2 Atomic structure of the spherical Si QD in a Ge matrix. A slice through the center of the QD and normal of the \(x\) axis is shown. Si (Ge) atoms are represented as blue (red) spheres.
5.4.3. Atomistic electronic structure
Having obtained the relaxed atomic structure, we may now resolve the electronic
structure of the QD using the atomistic tight-binding model. We first impose
an external electric field of magnitude \(10^7\,\mathrm{V/m}\)
pointing in the \(-z\) direction. To do so, we impose a corresponding
electronic potential on the atomic structure using the
set_potential method.
# Impose electric potential corresponding to an electric field of 10 V/μm
# pointing in the -z direction
def phi(x, y, z):
return 1e7 * z
atoms.set_potential(phi=phi)
Such an electric field may arise, for example, through a voltage applied on a back gate.
Since the wavefunctions of electrons confined in the QD are expected to penetrate only moderately into the Ge matrix, we define a simulation domain for the atomistic Schrödinger equation that is smaller than the entire atomic structure. Specifically, we define the simulation domain to be a cube of side length \(6\,\mathrm{nm}\) centered at the origin.
# Define geometry of box where the Schrödinger equation is solved
box_size_tb = diameter * 1.5
x_tb = np.array([-box_size_tb/2, box_size_tb/2])
y_tb = x_tb.copy()
z_tb = x_tb.copy()
We may then proceed to solve the atomistic Schrödinger equation within this simulation domain using the procedure seen in Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting.
# Solve the Schrödinger equation without periodic boundary conditions
subatoms = SubAtoms(parent=atoms, x=x_tb, y=y_tb, z=z_tb, pbc_directions=[])
schrodinger_solver_params = SchrodingerSolverParams()
schrodinger_solver_params.verbose = True
schrodinger_solver_params.num_states = 5
schrodinger_solver = SchrodingerSolver(atoms=subatoms,
solver_params=schrodinger_solver_params)
energy_target = 0.7 * ct.e
schrodinger_solver.solve(energy_target=energy_target)
subatoms.print_energies()
analyze_dot(atoms=subatoms, psi=subatoms.eigenfunctions[0], verbose=True)
vtu_name = "spherical_dot.vtu"
vtu_path = script_dir / "output" / vtu_name
out_dict = {
"Potential": subatoms.get_potential(),
"Ground state": subatoms.eigenfunctions[0].get_prob_on_atoms(),
}
save_vtu(atoms=subatoms, out_dict=out_dict, path=vtu_path)
Importantly, we did not apply PBCs when solving the Schrödinger equation. This
was done by setting the pbc_directions argument of the
SubAtoms constructor method to an
empty list ([]). Removing PBCs ensures that the wavefunctions of the
electrons confined in the QD do not interact with spurious periodic images of
themselves. We also note that the value of the energy_target argument passed
to the solve method was set to
\(0.7\,\mathrm{eV}\), which is just below the conduction band edge of bulk
Si. This ensures that the solver focuses on resolving states near the
conduction-band edge of the Si QD.
The print_energies method
prints the energies of the resolved states to the console.
Atomistic tight-binding eigenenergies (eV):
[0.91119092 0.91157244 0.91207730 0.91468245 0.91483111]
The analyze_dot function prints
various geometric properties of the ground state wavefunction to the console.
--------------------------------------------------------------------------------
Geometric properties of the tight-binding wavefunction:
--------------------------------------------------------------------------------
position: [-0.03371784 -0.03340513 1.4928158 ] A
std: [ 8.7068183 8.70857933 10.01192838] A
size: [34.82727322 34.83431732 40.04771352] A
--------------------------------------------------------------------------------
Finally, the save_vtu function saves the
atomic structure, potential, and ground state probability density to a VTU
file for visualization in Paraview. Below, we show the ground state probability
density visualized in Paraview.
Fig. 5.4.3 Ground state probability density of the spherical Si QD in a Ge matrix. A slice through the center of the QD and normal of the \(x\) axis is shown. Red (blue) indicates high (low) probability density.
5.5. Modeling a point charge defect
Having obtained the electronic structure of an ideal spherical Si QD in a Ge matrix, we may now investigate the impact of atomistic disorder on the QD’s electronic properties. We start by considering the impact of a point charge defect located at the center of the Si QD. Such point charge defects may arise from impurities or vacancies in the crystal lattice.
We model a positive point charge defect through the following Coulomb-like potential.
where \(e\) is the charge of the defect, \(\varepsilon\) is the permittivity of Si, \(R=1\,\mathrm{\mathring{A}}\) is a small cutoff radius to avoid the singularity at \(r=0\), and \(r\) is the distance from the defect located at the origin.
To include this defect potential in our simulation, we first define a function
phi_defect that implements the potential given by (5.5.1). We
then plot it, and add it to the existing potential of the atomic structure
using the add_to_potential
method.
# Define and plot potential describing a positive charge defect at the origin
eps = Si.eps
def phi_defect(x, y, z):
R = 1e-10
r = np.sqrt(x**2 + y**2 + z**2)
if r < R:
return ct.e / (4*np.pi*eps*R)
else:
return ct.e / (4*np.pi*eps*r)
r_num = 100
r_vec = np.linspace(0, box_size_tb/2, r_num)
phi_vec = np.array([phi_defect(r, 0, 0) for r in r_vec])
plt.figure()
plt.plot(r_vec*1e9, phi_vec)
plt.xlabel("Distance from center of QD (nm)")
plt.ylabel("Electric potential (V)")
plt.title("Electric potential due to positive charge defect in spherical QD")
plt.grid()
plt.show()
atoms.add_to_potential(phi=phi_defect)
Below, we show the resulting plot of the electric potential associated with the defect.
Fig. 5.5.2 Electric potential due to positive charge defect in spherical Si QD.
The electronic structure of the QD in the presence of the defect may then be resolved using the same procedure seen earlier.
# Solve the Schrödinger equation with a positive charge defect at the center of
# the dot
subatoms_defect = SubAtoms(parent=atoms, x=x_tb, y=y_tb, z=z_tb,
pbc_directions=[])
schrodinger_solver_defect = SchrodingerSolver(atoms=subatoms_defect,
solver_params=schrodinger_solver_params)
schrodinger_solver_defect.solve(energy_target=energy_target)
subatoms_defect.print_energies()
analyze_dot(atoms=subatoms_defect, psi=subatoms_defect.eigenfunctions[0],
verbose=True)
vtu_name = "spherical_dot_defect.vtu"
vtu_path = script_dir / "output" / vtu_name
out_dict = {
"Potential": subatoms_defect.get_potential(),
"Ground state": subatoms_defect.eigenfunctions[0].get_prob_on_atoms(),
}
save_vtu(atoms=subatoms_defect, out_dict=out_dict, path=vtu_path)
We note that the eigenenergies are around \(0.1\,\mathrm{eV}\) lower than those seen earlier for the ideal QD without the defect. This is expected, since the positive charge defect attracts electrons, thereby lowering their energies.
Atomistic tight-binding eigenenergies (eV):
[0.79645030 0.79695321 0.79717869 0.79999066 0.80042950]
Similarly, the size of the ground state wavefunction is seen to decrease relative to that of the ideal QD, indicating that the wavefunction is more tightly confined due to the presence of the defect.
--------------------------------------------------------------------------------
Geometric properties of the tight-binding wavefunction:
--------------------------------------------------------------------------------
position: [0.0546929 0.05476448 0.94311042] A
std: [7.95284563 7.95338686 8.81115229] A
size: [31.81138251 31.81354744 35.24460914] A
--------------------------------------------------------------------------------
Below, we show the ground state probability density visualized in Paraview.
Fig. 5.5.3 Ground state probability density of the spherical Si QD in a Ge matrix in the presence of a positive charge defect. A slice through the center of the QD and normal of the \(x\) axis is shown. Red (blue) indicates high (low) probability density.
5.6. Modeling Si and Ge concentration profiles
Next, we investigate the impact of alloy concentration gradients on the electronic properties of the spherical Si QD in a Ge matrix. Such concentration gradients are known to occur in real devices due to interdiffusion of Si and Ge atoms during growth or annealing.
To do so, we first need to parametrize the Si concentration profile. A simple model consists in (1) defining a diffusion radius, here set to \(4\,\mathrm{\mathring{A}}\), (2) setting the Si concentration to be \(100\%\) inside a sphere of radius equal to the QD radius minus the diffusion radius, i.e. \(2\,\mathrm{nm}-0.4\,\mathrm{nm}=1.6\,\mathrm{nm}\), (3) setting the Si concentration profile to be \(0\%\) outside a sphere of radius equal to the QD radius plus the diffusion radius, i.e. \(2\,\mathrm{nm}+0.4\,\mathrm{nm}=2.4\,\mathrm{nm}\), and (4) linearly varying the Si concentration between these two radii. Below, we define and plot this concentration profile.
# Define and plot functions specifying Si and Ge concentration profiles
diffusion_radius = diameter / 10
def Si_concentration(x, y, z):
r = np.sqrt(x**2 + y**2 + z**2)
if r <= diameter/2 - diffusion_radius:
return 1.0
elif r >= diameter/2 + diffusion_radius:
return 0.0
else:
a = -1 / (2 * diffusion_radius)
b = (diameter/2 + diffusion_radius) / (2 * diffusion_radius)
return a * r + b
def Ge_concentration(x, y, z):
return 1.0 - Si_concentration(x, y, z)
Si_vec = np.array([Si_concentration(r, 0, 0) for r in r_vec])
Ge_vec = np.array([Ge_concentration(r, 0, 0) for r in r_vec])
plt.figure()
plt.plot(r_vec*1e9, Si_vec, label="Si")
plt.plot(r_vec*1e9, Ge_vec, label="Ge")
plt.xlabel("Distance from center of QD (nm)")
plt.ylabel("Concentration of chemical species")
plt.title("Si and Ge concentration profiles in spherical QD")
plt.legend()
plt.grid()
plt.show()
We note that the Ge concentration profile is simply given by one minus the Si concentration profile. The Si and Ge concentration profiles are shown below.
Fig. 5.6.7 Si and Ge concentration profiles in spherical Si QD.
The Si and Ge atoms will then be randomly distributed in the atomic structure within the constraints set by the given concentration profiles. The randomness of this distribution is set by a random number generation seed, which by default is determined by operating system entropy and does not need to be specified by the user. Here, we fix the seed to ensure reproducibility of this tutorial’s results. However, we stress that users do not need to set this seed for their simulations.
# Seed random number generator for reproducibility of the outputs of this
# tutorial; not required for simulations, in practice
rng_seed = 0
The atomic structure may then be built using the same procedure as before.
# Include gradients of Si and Ge concentration around in the dot-matrix
# interface
atoms_gradient = Atoms(x=x, y=y, z=z, unit_cell=unit_cell,
pbc_directions=pbc_directions, rng_seed=rng_seed)
atoms_gradient.new_region(atom_species=np.array([["Si", Si_concentration], \
["Ge", Ge_concentration]]))
This time, the input atom_species of the
new_region method is passed as a
2D array, where each row specifies an atomic species and its corresponding
concentration profile function.
Having done that, we may relax the atomic structure and solve the Schrödinger equation using the same procedures as before.
# Relax atomic structure and solve the Schrödinger equation
keating_solver_gradient = KeatingSolver(atoms=atoms_gradient,
solver_params=keating_solver_params)
keating_solver_gradient.solve()
atoms_gradient.set_potential(phi=phi)
subatoms_gradient = SubAtoms(parent=atoms_gradient, x=x_tb, y=y_tb, z=z_tb,
pbc_directions=[])
schrodinger_solver_gradient = SchrodingerSolver(atoms=subatoms_gradient,
solver_params=schrodinger_solver_params)
schrodinger_solver_gradient.solve(energy_target=energy_target)
subatoms_gradient.print_energies()
analyze_dot(atoms=subatoms_gradient, psi=subatoms_gradient.eigenfunctions[0],
verbose=True)
vtu_name = "spherical_dot_gradient.vtu"
vtu_path = script_dir / "output" / vtu_name
out_dict = {
"Potential": subatoms_gradient.get_potential(),
"Ground state": subatoms_gradient.eigenfunctions[0].get_prob_on_atoms(),
}
save_vtu(atoms=subatoms_gradient, out_dict=out_dict, path=vtu_path)
Below, we show the atomic structure resulting from this procedure.
Fig. 5.6.8 Atomic structure of the spherical Si QD in a Ge matrix with Si and Ge concentration gradients. A slice through the center of the QD and normal of the \(x\) axis is shown. Si (Ge) atoms are represented as blue (red) spheres.
This atomic structure exhibits Si and Ge concentration gradients at the QD–matrix interface, as intended. The eigenenergies and geometric properties of the ground state wavefunction are shown below.
Atomistic tight-binding eigenenergies (eV):
[0.90907681 0.91038083 0.91109460 0.91387806 0.91520954]
--------------------------------------------------------------------------------
Geometric properties of the tight-binding wavefunction:
--------------------------------------------------------------------------------
position: [ 0.06652204 -0.36254672 1.24659106] A
std: [10.22386865 9.18120988 8.2390148 ] A
size: [40.89547459 36.7248395 32.95605921] A
--------------------------------------------------------------------------------
While the eigenenergies are similar to those seen earlier for the ideal QD without concentration gradients, the size of the ground state wavefunction is seen to increase, indicating that the wavefunction is less tightly confined due to the presence of the concentration gradients. This is also seen visually below.
Fig. 5.6.9 Ground state probability density of the spherical Si QD in a Ge matrix with Si and Ge concentration gradients. A slice through the center of the QD and normal of the \(x\) axis is shown. Red (blue) indicates high (low) probability density.
5.7. Full code
__copyright__ = "Copyright 2022-2025, Nanoacademic Technologies Inc."
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh
from qtcad.device.materials import Si
from qtcad.atoms import Atoms, SubAtoms
from qtcad.atoms.unit_cell import UnitCellZincblende
from qtcad.atoms import Atoms
from qtcad.atoms.keating import Solver as KeatingSolver
from qtcad.atoms.keating import SolverParams as KeatingSolverParams
from qtcad.atoms.schrodinger import Solver as SchrodingerSolver
from qtcad.atoms.schrodinger import SolverParams as SchrodingerSolverParams
from qtcad.atoms.analysis import save_vtu, analyze_dot
# Geometry of the spherical Si quantum dot in a Ge matrix
diameter = 4e-9
box_size = diameter * 2
x = np.array([-box_size/2, box_size/2])
y = x.copy()
z = x.copy()
# Unit cell of bulk Ge
lattice_constant = 5.647e-10
unit_cell = UnitCellZincblende(lattice_constant=lattice_constant,
atom_species_ws=np.array(["Ge", "Ge"]))
# Construct atomic structure from mesh
script_dir = Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes" / "spherical_dot.msh"
mesh = Mesh(scaling_factor=diameter, filename=str(path_mesh))
atoms_from_mesh = Atoms(mesh=mesh, unit_cell=unit_cell)
atoms_from_mesh.new_region(atom_species=np.array(["Si"]*8), tag="ball")
path_vtu = script_dir / "output" / "spherical_dot_from_mesh.vtu"
save_vtu(atoms=atoms_from_mesh, out_dict={}, path=path_vtu)
# Axes along which periodic boundary conditions are applied for the atomic
# structure relaxation
pbc_directions = ["x", "y", "z"]
# Create atomic structure, defining the Si spherical dot via the `is_inside`
# argument of the `new_region` method
atoms = Atoms(x=x, y=y, z=z, unit_cell=unit_cell,
pbc_directions=pbc_directions)
def is_inside_dot(x, y, z):
r = np.sqrt(x**2 + y**2 + z**2)
return r <= diameter / 2
atoms.new_region(atom_species=np.array(["Si"]*8), is_inside=is_inside_dot)
# Relax atomic structure with periodic boundary conditions
keating_solver_params = KeatingSolverParams()
keating_solver_params.verbose = True
keating_solver = KeatingSolver(atoms=atoms,
solver_params=keating_solver_params)
keating_solver.solve()
# Impose electric potential corresponding to an electric field of 10 V/μm
# pointing in the -z direction
def phi(x, y, z):
return 1e7 * z
atoms.set_potential(phi=phi)
# Define geometry of box where the Schrödinger equation is solved
box_size_tb = diameter * 1.5
x_tb = np.array([-box_size_tb/2, box_size_tb/2])
y_tb = x_tb.copy()
z_tb = x_tb.copy()
# Solve the Schrödinger equation without periodic boundary conditions
subatoms = SubAtoms(parent=atoms, x=x_tb, y=y_tb, z=z_tb, pbc_directions=[])
schrodinger_solver_params = SchrodingerSolverParams()
schrodinger_solver_params.verbose = True
schrodinger_solver_params.num_states = 5
schrodinger_solver = SchrodingerSolver(atoms=subatoms,
solver_params=schrodinger_solver_params)
energy_target = 0.7 * ct.e
schrodinger_solver.solve(energy_target=energy_target)
subatoms.print_energies()
analyze_dot(atoms=subatoms, psi=subatoms.eigenfunctions[0], verbose=True)
vtu_name = "spherical_dot.vtu"
vtu_path = script_dir / "output" / vtu_name
out_dict = {
"Potential": subatoms.get_potential(),
"Ground state": subatoms.eigenfunctions[0].get_prob_on_atoms(),
}
save_vtu(atoms=subatoms, out_dict=out_dict, path=vtu_path)
# Define and plot potential describing a positive charge defect at the origin
eps = Si.eps
def phi_defect(x, y, z):
R = 1e-10
r = np.sqrt(x**2 + y**2 + z**2)
if r < R:
return ct.e / (4*np.pi*eps*R)
else:
return ct.e / (4*np.pi*eps*r)
r_num = 100
r_vec = np.linspace(0, box_size_tb/2, r_num)
phi_vec = np.array([phi_defect(r, 0, 0) for r in r_vec])
plt.figure()
plt.plot(r_vec*1e9, phi_vec)
plt.xlabel("Distance from center of QD (nm)")
plt.ylabel("Electric potential (V)")
plt.title("Electric potential due to positive charge defect in spherical QD")
plt.grid()
plt.show()
atoms.add_to_potential(phi=phi_defect)
# Solve the Schrödinger equation with a positive charge defect at the center of
# the dot
subatoms_defect = SubAtoms(parent=atoms, x=x_tb, y=y_tb, z=z_tb,
pbc_directions=[])
schrodinger_solver_defect = SchrodingerSolver(atoms=subatoms_defect,
solver_params=schrodinger_solver_params)
schrodinger_solver_defect.solve(energy_target=energy_target)
subatoms_defect.print_energies()
analyze_dot(atoms=subatoms_defect, psi=subatoms_defect.eigenfunctions[0],
verbose=True)
vtu_name = "spherical_dot_defect.vtu"
vtu_path = script_dir / "output" / vtu_name
out_dict = {
"Potential": subatoms_defect.get_potential(),
"Ground state": subatoms_defect.eigenfunctions[0].get_prob_on_atoms(),
}
save_vtu(atoms=subatoms_defect, out_dict=out_dict, path=vtu_path)
# Define and plot functions specifying Si and Ge concentration profiles
diffusion_radius = diameter / 10
def Si_concentration(x, y, z):
r = np.sqrt(x**2 + y**2 + z**2)
if r <= diameter/2 - diffusion_radius:
return 1.0
elif r >= diameter/2 + diffusion_radius:
return 0.0
else:
a = -1 / (2 * diffusion_radius)
b = (diameter/2 + diffusion_radius) / (2 * diffusion_radius)
return a * r + b
def Ge_concentration(x, y, z):
return 1.0 - Si_concentration(x, y, z)
Si_vec = np.array([Si_concentration(r, 0, 0) for r in r_vec])
Ge_vec = np.array([Ge_concentration(r, 0, 0) for r in r_vec])
plt.figure()
plt.plot(r_vec*1e9, Si_vec, label="Si")
plt.plot(r_vec*1e9, Ge_vec, label="Ge")
plt.xlabel("Distance from center of QD (nm)")
plt.ylabel("Concentration of chemical species")
plt.title("Si and Ge concentration profiles in spherical QD")
plt.legend()
plt.grid()
plt.show()
# Seed random number generator for reproducibility of the outputs of this
# tutorial; not required for simulations, in practice
rng_seed = 0
# Include gradients of Si and Ge concentration around in the dot-matrix
# interface
atoms_gradient = Atoms(x=x, y=y, z=z, unit_cell=unit_cell,
pbc_directions=pbc_directions, rng_seed=rng_seed)
atoms_gradient.new_region(atom_species=np.array([["Si", Si_concentration], \
["Ge", Ge_concentration]]))
# Relax atomic structure and solve the Schrödinger equation
keating_solver_gradient = KeatingSolver(atoms=atoms_gradient,
solver_params=keating_solver_params)
keating_solver_gradient.solve()
atoms_gradient.set_potential(phi=phi)
subatoms_gradient = SubAtoms(parent=atoms_gradient, x=x_tb, y=y_tb, z=z_tb,
pbc_directions=[])
schrodinger_solver_gradient = SchrodingerSolver(atoms=subatoms_gradient,
solver_params=schrodinger_solver_params)
schrodinger_solver_gradient.solve(energy_target=energy_target)
subatoms_gradient.print_energies()
analyze_dot(atoms=subatoms_gradient, psi=subatoms_gradient.eigenfunctions[0],
verbose=True)
vtu_name = "spherical_dot_gradient.vtu"
vtu_path = script_dir / "output" / vtu_name
out_dict = {
"Potential": subatoms_gradient.get_potential(),
"Ground state": subatoms_gradient.eigenfunctions[0].get_prob_on_atoms(),
}
save_vtu(atoms=subatoms_gradient, out_dict=out_dict, path=vtu_path)