2. Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting

2.1. Requirements

2.1.1. Software components

  • QTCAD

  • Paraview

  • VESTA

2.1.2. Python script

  • qtcad/examples/tutorials/multiscale_2.py

2.1.3. References

2.2. Briefing

In the previous tutorial of this tutorial series, we resolved the device electrostatics of a SiGe–Si–SiGe quantum dot (QD) using the linear Poisson solver. We then solved the Schrödinger equation in the QD using the effective mass approximation. Within this approximation, the six “valleys” at which the conduction bands of silicon are minima are treated as equivalent. Thus, each eigenstate computed within the effective mass approximation has a six-fold degeneracy. However, in reality, the six valleys are not equivalent. Indeed, the strong quantum confinement that arises in QDs, electric fields, and atomistic surface roughness typically lift this degeneracy. The difference between the two lowest energy levels of the QD, corresponding to two different valley states, is known as the valley splitting. Ensuring that the valley splitting is sufficiently high is a key requirement for the operation of certain silicon-based quantum computers.

There are three main computational approaches to model valley splitting in Si QDs: the multivalley effective-mass theory (MVEMT), the tight-binding (TB) model, and density functional theory (DFT). TB calculations strike the right balance of physical realism, computational tractability, and direct atomistic detail needed to accurately capture valley splitting in semiconductor nanostructures. Unlike MVEMT, which simplifies the band structure and cannot easily accommodate interface complexity or subtle atomistic details, and DFT, which is computationally expensive for large systems, the TB model incorporate the necessary local bonding environments and symmetries while remaining manageable for large device geometries. As a result, the TB approach better reflects how local atomistic variations—such as surface roughness, strain, and random alloying—modify the degeneracy and separation of valley states, making it the preferred method for realistic valley splitting simulations.

The goals of this tutorial are as follows:

  • construct a realistic atomic structure describing the SiGe–Si–SiGe QD, including atomistic nonidealities and randomness arising from strain, random alloying, and surface roughness;

  • construct a TB Hamiltonian for the QD, including the external potential created by the gate electrodes, and solve the Schrödinger equation within the TB model

  • compute the valley splitting of the QD and compute the corresponding valley states;

  • extract from the TB calculations atomistic parameters needed for the quantum control simulations of the next tutorial.

Note

This tutorial involves a finite-element mesh containing around one million nodes. We expect the simulation presented here to run in a few minutes on a relatively advanced desktop workstation (as described in Software requirements). Less powerful machine may have memory constraints when running this simulation.

2.4. Constructing the atomic structure

2.4.1. Unit cell and bandstructure

We start by harnessing the results of the previous tutorial to create a finite-element device describing the SiGe–Si–SiGe QD, including the electrostatic potential created by the gate electrodes.

# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_refined_mesh = str(script_dir / "meshes" / "multiscale.msh")
scaling = 1e-9
mesh = Mesh(scaling, path_refined_mesh)

# Create device and set potential in device to result of FEM simulation
d = Device(mesh, conf_carriers="e")
path_phi_hdf5 = str(script_dir / "output" / "multiscale_phi.hdf5")
d.set_potential(io.load(path_phi_hdf5, "phi"))

Next, we construct the unit cells of Si and Ge. To do so, we create a UnitCellZincblende object, which describes the primitive unit cell of a material with diamond or zincblende crystal structure. To instantiate such an object, two variables are required: the lattice constant of the material (namely, the side length of its cubic unit cell) and an array containing the chemical symbols of the two atoms in the Wigner–Seitz cell. Here, we do so for Si and Ge.

# Create unit cells of Si and Ge
unit_cell_Si = UnitCellZincblende(lattice_constant=5.431*1e-10,
    atom_species_ws=np.array(["Si", "Si"]))
unit_cell_Ge = UnitCellZincblende(lattice_constant=5.658*1e-10,
    atom_species_ws=np.array(["Ge", "Ge"]))

The bandstructure of these two materials may be computed and plotted using the plot_bandstructure method. Note that the bandstructure calculation is based on the TB model.

# Plot band structures of Si and Ge
path_bs_Si = str(script_dir / "output" / "multiscale_bs_Si.png")
unit_cell_Si.plot_bandstructure(title="Si", path=path_bs_Si, show_figure=False)
path_bs_Ge = str(script_dir / "output" / "multiscale_bs_Ge.png")
unit_cell_Ge.plot_bandstructure(title="Ge", path=path_bs_Ge, show_figure=False)

The plots that result from this operation are shown below.

Bandstructure of Si computed within the TB model.

Fig. 2.4.6 Bandstructure of Si computed within the TB model.

Bandstructure of Ge computed within the TB model.

Fig. 2.4.7 Bandstructure of Ge computed within the TB model.

The bandstructure plots may be inspected to reveal that the bandgap of Si is predicted from the TB model to be around \(1.19\;\mathrm{eV}\), while that of Ge is predicted to be around \(0.83\;\mathrm{eV}\). These values are in good agreement with experimental data. In addition, while the valence band maximum of Si lies around \(-0.60\;\mathrm{eV}\), that of Ge lies around \(0.0\;\mathrm{eV}\), which reflects the band offset between these two materials.

2.4.2. Atomic structure of the quantum dot

We now construct the atomic structure of the SiGe–Si–SiGe QD. To do so, we create an Atoms object. To instantiate it, we specify the boundaries of a box in space to be filled with atoms. Here, we choose this box to fully encapsulates the SiGe spacer, Si quantum well (QW), and SiGe buffer layers of the heterostructure. Along the plane normal to the heterostructure growth axis, we set the boundaries of the box to encapsulate the ground state that was computed in the previous tutorial. These boundaries are specified in three arrays, x, y, and z, which contain the minimum and maximum values of the box along the \(x\), \(y\), and \(z\) axes, respectively. We also set the random number generator seed (rng_seed), which is used to generate the random SiGe alloys in the atomic structure. Note that a typical simulation does not require the user to set this seed, as it is, by default, automatically set to a value obtained via from fresh, unpredictable entropy from the operating system. Here, we set the seed so as to ensure reproducibility of the results of the tutorial.

# Create atomic structure for heterostructure
x = np.array([-10., 10.]) * 1e-9
y = np.array([-10., 10.]) * 1e-9
z = np.array([-50., -5.]) * 1e-9
rng_seed = 100
atoms = Atoms(x=x, y=y, z=z, rng_seed=100)

To model some of the nonidealities and randomness that arise in real atomistic structures, we create a random rough surface to model the interface between the Si cap and SiGe spacer layers. To do so, we create a RoughSurface object. We specify the Hurst exponent of the surface (which may be estimated from experimental transmission electron microscopy images [CTG+24]), its mean height (mean), and its root-mean-square displacement from the mean height (rms).

# Create rough surface describing top and bottom surfaces of the well
well_top = RoughSurface(hurst=0.3, mean=-35.*1e-9, rms=5.*1e-10, rng_seed=0)
well_bot = RoughSurface(hurst=0.3, mean=-38.*1e-9, rms=5.*1e-10, rng_seed=1)

Note that the rough surface is generated using a random number generator. We set the seed of this generator to the same value as the seed of the atomic structure builder through the keyword argument rng_seed. Again, this does not need to be done in a typical simulation. We do it here to ensure reproducibility.

To visualize this rough surface, we call the plot method. Either 2D or 3D plots may be generated, as is done below.

# Plot rough surface
x_plot = np.linspace(x[0], x[1], 100)
y_plot = np.linspace(y[0], y[1], 100)
path_rough_2d = str(script_dir / "output" / "multiscale_rough_2d.png")
well_top.plot(x=x_plot, y=y_plot, type="2D", path=path_rough_2d,
    show_figure=False)
path_rough_3d = str(script_dir / "output" / "multiscale_rough_3d.png")
well_top.plot(x=x_plot, y=y_plot, type="3D", path=path_rough_3d,
    show_figure=False)

This results in the following images.

Contour plot of the rough interface between the SiGe spacer and the Si well.

Fig. 2.4.8 Contour plot of the rough interface between the SiGe spacer and the Si well.

Height profile of the rough interface between the SiGe spacer and the Si well.

Fig. 2.4.9 Height profile of the rough interface between the SiGe spacer and the Si well.

The rough surface exhibits fractal character. The Hurst exponent of the rough surface, \(H\), is related to the surface’s fractal dimension \(D_f\) through \(D_f = 3 - H\).

Another aspect of atomistic randomness that we model is the random alloying in the SiGe spacer and SiGe buffer layers. To do so, we first define a random alloy with \(70\%\) Si and \(30\%\) Ge concentrations, and set two layers of the heterostructure to be filled with this random alloy. This is done via the new_layer method.

# Define the material stack of the heterostructure
SiGe_alloy = np.array([["Si", 0.7], ["Ge", 0.3]])
atoms.new_layer(z_top=-5.*1e-9, z_bot=well_top, atom_species=SiGe_alloy) # cap
atoms.new_layer(z_top=well_bot, z_bot=z[0], atom_species=SiGe_alloy) # buffer

Next, the atomic structure may be relaxed using the Keating valence force-field model. This is done via the Solver and SolverParams classes of the qtcad.atoms.keating module.

# Relax atomic coordinates
keating_solver_params = KeatingSolverParams()
keating_solver_params.verbose = True
keating_solver = KeatingSolver(atoms=atoms,
    solver_params=keating_solver_params)
keating_solver.solve()
path_xyz = str(script_dir / "output" / "multiscale.xyz")
atoms.save_xyz(path=path_xyz)

The .xyz file of the atomic structure was saved via the save_xyz method and may be opened by an atomic structure visualization software such as VESTA.

Atomic structure visualized in VESTA.

Fig. 2.4.10 Atomic structure visualized in VESTA.

Here, the Si atoms are shown in blue and the Ge atoms are shown in green. The Si QW layer, which is uniformly blue-colored, is clearly visible. Its rough interfaces with the SiGe layers is also apparent. In addition, the Ge atoms are randomly scattered throughout the SiGe spacer and SiGe buffer layers, as expected from the random alloying that we set up in the structure. We also notice a slight bulging of the SiGe layers compared to the Si layers. This is due to the fact that the lattice constant of Ge is greater than that of Si, which causes the SiGe layers to expand. In turn, this leads to strain in the atomic structure.

2.5. Solving the Schrödinger equation within the TB model

While strain may still be significant several tens or even hundreds of nanometers away from a heterointerface, the QD wavefunctions only spread a few nanometers, as we saw in the previous tutorial. Thus, to reduce the computational cost of the TB simulation of the QD with minimal impact on accuracy, we construct a restricted atomic structure that encapsulates atoms in and close to the Si QD. This is done via the SubAtoms class.

# Contruct restricted atomic structure for QD
x_QD = np.array([-8., 8.]) * 1e-9
y_QD = np.array([-8., 8.]) * 1e-9
z_QD = np.array([-40., -33.]) * 1e-9
atoms_QD = SubAtoms(parent=atoms, x=x_QD, y=y_QD, z=z_QD)

Here, we restricted the atomic structure along the \(x\), \(y\) and \(z\) axes to the spread of the ground state wavefunction of the QD computed in the previous tutorial. This is done through the inputs x, y, and z. We note that the original atomic structure contains \(908,701\) atoms, while the restricted atomic structure contains \(86,700\) atoms, as is automatically printed to screen when the structures are built.

Total number of atoms in atomic structure: 908701.
...
Total number of atoms in restricted atomic structure: 86700.

The atoms and atoms_QD objects both have attributes atom_species and atom_pos, which are arrays that respectively contains the chemical species of the atoms in the atomic structure and the positions of these atoms.

Having obtained an atomic structure for the QD, we can now compute the TB Hamiltonian of the QD and obtain a few eigenpairs close to a specified target energy.

To do so, we first set the electric potential in the atomic structure to be determined through the results of the finite-element simulations of the previous tutorial. This is done via the set_potential_from_device method, which takes as input the device in question (d). Next, we set up the parameters of the TB Schrödinger equation solver through the SolverParams. This object contains several attributes, such as the number of eigenpairs to be computed (num_states), and whether messages pertaining to the solver are printed to the screen (verbose). We then create a Solver object, which takes as input the atomic structure (atoms) and the solver parameters (solver_params).

# Create Schrödinger equation solver within TB formalism
atoms_QD.set_potential_from_device(d=d)
schrodinger_solver_params = SchrodingerSolverParams()
schrodinger_solver_params.verbose = True
schrodinger_solver_params.num_states = 10
schrodinger_solver = SchrodingerSolver(atoms=atoms_QD,
    solver_params=schrodinger_solver_params)

We can then call the solve method of schrodinger_solver, which internally computes the TB Hamiltonian of the atomic structure. It is expressed in a basis of atomic orbitals. By default, there are ten orbitals per atom: \(s\), \(p_x\), \(p_y\), \(p_z\), \(d_{yz}\), \(d_{xz}\), \(d_{xy}\), \(d_{x^2-y^2}\), \(d_{3z^2-r^2}\), and \(s^{\star}\). This matrix, while sparse, is very large. In the case of this tutorial, it is a \(867,000 \times 867,000\) matrix. We thus use a specialized sparse eigensolver to extract a few eigenpairs close to a specified target energy, the energy_target input of the solve method. Ideally, this target energy should be set to lie within the bandgap of the QD, close to the conduction band edge (valence band edge) when studying the properties of electrons (holes). From the previous tutorial, the conduction band edge in the center of the QD is around \(0\;\mathrm{eV}\). We thus set the target energy to \(-0.1\;\mathrm{eV}\).

# Solve the Schrödinger equation
energy_target = -0.1 * ct.e
schrodinger_solver.solve(energy_target=energy_target)
atoms_QD.print_energies()

We obtain the following eigenenergies.

Atomistic tight-binding eigenenergies (eV):
[0.01009784 0.01043149 0.03336937 0.03356846 0.03705035 0.03719805 0.05708449 0.05717669 0.06442718 0.06481970]

The eigenenergies and eigenfunctions obtained via this procedure are respectively stored in the energies and eigenfunctions attributes of atoms_QD.

2.6. Analyzing the TB simulation results

2.6.1. Orbital, chemical, and geometrical properties

Having obtained eigenpairs of the TB Hamiltonian, we can now compute parameters that are used to model quantum control in the QD. The first of these is the valley splitting.

# Compute valley splitting
valley_splitting = atoms_QD.energies[1] - atoms_QD.energies[0]
print("Valley splitting: " + str(valley_splitting/ct.e) + " eV")

We selected the target energy for the eigensolver, energy_target, to lie within the band gap of the Si QD. As such the two first eigenpairs correspond to the ground and first-excited conduction-band states. Their difference is therefore the valley splitting. We obtain the following valley splitting.

Valley splitting: 0.0003336573366181233 eV

Next, we analyze the orbital, chemical, and geometrical properties of the two valley states. We do so by calling the get_prob_on_orbitals and get_prob_on_chemical_species methods of the <Wavefunction>qtcad.atoms.wavefunction.Wavefunction class, as well as the analyze_dot function. They all take as input verbose, a boolean that determines whether the results are printed to screen. In order to project the wavefunction onto the chemical species present in the system, the get_prob_on_chemical_species method also takes as input atoms, which contains the chemical species of the atoms in the QD. The analyze_dot function also takes as input atoms, which contains the positions of the atoms in the QD, and psi, which contains the wavefunction to be analyzed. They respectively compute the modulus-square of the projections of the states on the ten TB model atomic orbitals, the modulus-square of the projections of the state on the chemical species present in the system, and the geometrical properties (mean position, standard deviation, and size) of the wavefunction.

# Store valley states in variables nu_1 and nu_2
nu_1 = atoms_QD.eigenfunctions[0]
nu_2 = atoms_QD.eigenfunctions[1]

# Compute modulus-square of projections of ground state on the
# spds* orbitals and on the chemical species in the system
# Also compute the geometric properties of the ground state wavefunction
print("Ground state (nu_1)")
nu_1.get_prob_on_orbitals(verbose=True)
nu_1.get_prob_on_chemical_species(atoms=atoms_QD, verbose=True)
analyze_dot(atoms=atoms_QD, psi=nu_1, verbose=True)

# Repeat these steps for the first excited state
print("First excited state (nu_2)")
nu_2.get_prob_on_orbitals(verbose=True)
nu_2.get_prob_on_chemical_species(atoms=atoms_QD, verbose=True)
analyze_dot(atoms=atoms_QD, psi=nu_2, verbose=True)

For example, for the first valley state, we obtain the following results.

Ground state (nu_1)
--------------------------------------------------------------------------------
Probabilities of finding an electron on each type of orbital
s: 0.14526646909336854
px: 0.0015057940807302738
py: 0.0015468140478460816
pz: 0.4040505962758415
dyz: 0.00010996718753399113
dxz: 0.0001312717656277609
dxy: 0.1876538808865558
d(x^2-y^2): 0.00013019769842256083
d(3z^2-r^2): 0.22947991412781676
s*: 0.030125094836253388
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Probabilities of finding an electron on each chemical species
Ge:  0.035485065249205805
Si:  0.9645149347507933
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Geometric properties of the tight-binding wavefunction:
--------------------------------------------------------------------------------
position:  [  -5.41243594   -0.692056   -366.34708394] A
std:  [27.6468206  26.07562881  7.28261706] A
size:  [110.58728239 104.30251526  29.13046825] A
--------------------------------------------------------------------------------

These results give several insights into the makeup of the state. For example, the ground state is of \(40\%\) \(p_z\) character, and only of \(0.15\%\) \(p_x\) and \(p_y\) character. This is expected, since the gate electrodes create an electric field parallel to the \(z\) axis in the QD, thereby energetically favoring the \(p_z\) orbital. In addition, the ground state is mostly localized on Si atoms, with \(96.5\%\) of its probability density on Si atoms and only \(3.5\%\) on Ge atoms. This is also expected, since due to the band offset between Si and Ge, Ge atoms have higher conduction band minimum, and are therefore less energetically favorable for the electron to occupy. Finally, the ground state wavefunction has a size of \(111\;\mathrm{\mathring{A}}\) along the \(x\) axis, \(104\;\mathrm{\mathring{A}}\) along the \(y\) axis, and \(29\;\mathrm{\mathring{A}}\) along the \(z\) axis. These sizes are smaller than those obtained in the finite-element simulations of the previous tutorial, which could be explained by a greater effective confinement arising as a result of the rough surfaces of the QW.

The wavefunctions may also be visualized in 3D. This is done by first computing the modulus-square of the projections of the states onto the atoms of the system by using the get_prob_on_atoms method. The results are then saved to a .vtu file by calling the save_vtu function.

# Save modulus-square of projection of states on atoms to a .vtu file
prob_on_atoms_0 = nu_1.get_prob_on_atoms()
prob_on_atoms_1 = nu_2.get_prob_on_atoms()
prob_on_atoms_2 = atoms_QD.eigenfunctions[2].get_prob_on_atoms()
path_vtu = str(script_dir / "output" / "multiscale_TB_wf.vtu")
out_dict = {"Probability density (ground state)": prob_on_atoms_0,
    "Probability density (1st excited state)": prob_on_atoms_1,
    "Probability density (2nd excited state)": prob_on_atoms_2}
save_vtu(atoms=atoms_QD, out_dict=out_dict, path=path_vtu)

This file may be opened in Paraview. To properly visualize the wavefunctions, we recommend setting the Representation parameter to Point Gaussian and the Gaussian radius to 1. We also recommend using the Clip setting to better visualize the probabilities on atoms in the middle of the atomic structure.

Probability density of the ground state wavefunction visualized in Paraview.

Fig. 2.6.3 Probability density of the ground state wavefunction visualized in Paraview.

Probability density of the first-excited state wavefunction visualized in Paraview.

Fig. 2.6.4 Probability density of the first-excited state wavefunction visualized in Paraview.

Probability density of the second-excited state wavefunction visualized in Paraview.

Fig. 2.6.5 Probability density of the second-excited state wavefunction visualized in Paraview.

The two valley states, namely the ground and first-excited states, have very similar probability densities. The second-excited state, on the other hand, exhibits a clear two-node structure. The positions of the Ge atoms in the QD may also be visualized in Paraview.

Ge atoms visualized in Paraview.

Fig. 2.6.6 Ge atoms visualized in Paraview.

2.6.2. Calculation of matrix elements

Besides the valley splitting, the quantum control model investigated in the third tutorial requires two additional parameters.

The first of these parameters is the matrix element \(\left<\nu_1\downarrow\right|H_{\text{SOC}}\left|\nu_2\uparrow\right>\), where \(H_{\text{SOC}}\) is the spin–orbit coupling term of the TB Hamiltonian. To compute this matrix element, we first “promote” the ground state nu_1 to have a spin down component using the add_spin_to_orbital_state method. A similar procedure is repeated for the nu_2 state, which is promoted to have a spin up component. We then compute the SOC Hamiltonian using the get_hamiltonian_soc method of the schrodinger_solver object. Finally, we compute the matrix element using the get_operator_element function, which takes as inputs the two states and the SOC Hamiltonian.

# Compute the matrix element of the spin-orbit coupling Hamiltonian between
# the valley states with opposite spins
H_soc = schrodinger_solver.get_hamiltonian_soc()
nu_2_up = nu_2.add_spin_to_orbital_state(spin=True, in_place=False)
nu_1_down = nu_1.add_spin_to_orbital_state(spin=False, in_place=False)
C = get_operator_element(psi_l=nu_2_up, op=H_soc, psi_r=nu_1_down)
print("Matrix element of SOC Hamiltonian: " + str(C/ct.e) + " eV")

Here, we obtain the following value.

Matrix element of SOC Hamiltonian: (7.887230152811525e-06+1.2796884642820679e-05j) eV

The second of these parameters is the matrix element \(\left<\nu_1\right|D_{\text{BG}}\left|\nu_2\right>\), where \(D_{\text{BG}}\) is the derivative of the confinement potential with respect to barrier gate voltage. This quantity was computed in the previous tutorial. Still, we need to convert it to an operator expressed in the same basis of atomic orbitals as the TB states nu_1 and nu_2. This is done by calling the get_operator_from_fem function.

# Compute the matrix element of the derivative of the potential with respect to
# barrier gate voltage between the valley states
qd = ["cap_qd", "well_qd", "buffer_qd"] # quantum dot regions
submesh = SubMesh(mesh, qd)
path_der_hdf5 = str(script_dir / "output" / "multiscale_der.hdf5")
der = io.load(path_der_hdf5, "der")
Dfg = get_operator_from_fem(atoms=atoms_QD, tb_orbs=nu_1.tb_orbs,
    tb_spin=nu_1.tb_spin, mesh=submesh, quantity=der)
D = get_operator_element(psi_l=nu_1, op=Dfg, psi_r=nu_2)
print("Matrix element of derivative of potential with respect to barrier " + \
    "gate voltage: " + str(D/ct.e) + " eV/V")

We obtain the following parameter value.

Matrix element of derivative of potential with respect to barrier gate voltage: (-3.049793596508124e-05-3.7349199660238024e-21j) eV/V

Finally, we save the parameters to a .txt file.

# Save parameters of quantum control problem to .txt file
params = np.array([valley_splitting, C, D])
path_params = str(script_dir / "output" / "multiscale_params.txt")
np.savetxt(path_params, params)

2.7. Full code

__copyright__ = "Copyright 2022-2025, Nanoacademic Technologies Inc."

import pathlib
import numpy as np
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import Device
from qtcad.device import io
from qtcad.atoms.unit_cell import UnitCellZincblende
from qtcad.atoms import Atoms, SubAtoms
from qtcad.atoms.rough_surface import RoughSurface
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 analyze_dot, save_vtu, get_operator_element, \
    get_operator_from_fem

# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_refined_mesh = str(script_dir / "meshes" / "multiscale.msh")
scaling = 1e-9
mesh = Mesh(scaling, path_refined_mesh)

# Create device and set potential in device to result of FEM simulation
d = Device(mesh, conf_carriers="e")
path_phi_hdf5 = str(script_dir / "output" / "multiscale_phi.hdf5")
d.set_potential(io.load(path_phi_hdf5, "phi"))

# Create unit cells of Si and Ge
unit_cell_Si = UnitCellZincblende(lattice_constant=5.431*1e-10,
    atom_species_ws=np.array(["Si", "Si"]))
unit_cell_Ge = UnitCellZincblende(lattice_constant=5.658*1e-10,
    atom_species_ws=np.array(["Ge", "Ge"]))

# Plot band structures of Si and Ge
path_bs_Si = str(script_dir / "output" / "multiscale_bs_Si.png")
unit_cell_Si.plot_bandstructure(title="Si", path=path_bs_Si, show_figure=False)
path_bs_Ge = str(script_dir / "output" / "multiscale_bs_Ge.png")
unit_cell_Ge.plot_bandstructure(title="Ge", path=path_bs_Ge, show_figure=False)

# Create atomic structure for heterostructure
x = np.array([-10., 10.]) * 1e-9
y = np.array([-10., 10.]) * 1e-9
z = np.array([-50., -5.]) * 1e-9
rng_seed = 100
atoms = Atoms(x=x, y=y, z=z, rng_seed=100)

# Create rough surface describing top and bottom surfaces of the well
well_top = RoughSurface(hurst=0.3, mean=-35.*1e-9, rms=5.*1e-10, rng_seed=0)
well_bot = RoughSurface(hurst=0.3, mean=-38.*1e-9, rms=5.*1e-10, rng_seed=1)

# Plot rough surface
x_plot = np.linspace(x[0], x[1], 100)
y_plot = np.linspace(y[0], y[1], 100)
path_rough_2d = str(script_dir / "output" / "multiscale_rough_2d.png")
well_top.plot(x=x_plot, y=y_plot, type="2D", path=path_rough_2d,
    show_figure=False)
path_rough_3d = str(script_dir / "output" / "multiscale_rough_3d.png")
well_top.plot(x=x_plot, y=y_plot, type="3D", path=path_rough_3d,
    show_figure=False)

# Define the material stack of the heterostructure
SiGe_alloy = np.array([["Si", 0.7], ["Ge", 0.3]])
atoms.new_layer(z_top=-5.*1e-9, z_bot=well_top, atom_species=SiGe_alloy) # cap
atoms.new_layer(z_top=well_bot, z_bot=z[0], atom_species=SiGe_alloy) # buffer

# Relax atomic coordinates
keating_solver_params = KeatingSolverParams()
keating_solver_params.verbose = True
keating_solver = KeatingSolver(atoms=atoms,
    solver_params=keating_solver_params)
keating_solver.solve()
path_xyz = str(script_dir / "output" / "multiscale.xyz")
atoms.save_xyz(path=path_xyz)

# Contruct restricted atomic structure for QD
x_QD = np.array([-8., 8.]) * 1e-9
y_QD = np.array([-8., 8.]) * 1e-9
z_QD = np.array([-40., -33.]) * 1e-9
atoms_QD = SubAtoms(parent=atoms, x=x_QD, y=y_QD, z=z_QD)

# Create Schrödinger equation solver within TB formalism
atoms_QD.set_potential_from_device(d=d)
schrodinger_solver_params = SchrodingerSolverParams()
schrodinger_solver_params.verbose = True
schrodinger_solver_params.num_states = 10
schrodinger_solver = SchrodingerSolver(atoms=atoms_QD,
    solver_params=schrodinger_solver_params)

# Solve the Schrödinger equation
energy_target = -0.1 * ct.e
schrodinger_solver.solve(energy_target=energy_target)
atoms_QD.print_energies()

# Compute valley splitting
valley_splitting = atoms_QD.energies[1] - atoms_QD.energies[0]
print("Valley splitting: " + str(valley_splitting/ct.e) + " eV")

# Store valley states in variables nu_1 and nu_2
nu_1 = atoms_QD.eigenfunctions[0]
nu_2 = atoms_QD.eigenfunctions[1]

# Compute modulus-square of projections of ground state on the
# spds* orbitals and on the chemical species in the system
# Also compute the geometric properties of the ground state wavefunction
print("Ground state (nu_1)")
nu_1.get_prob_on_orbitals(verbose=True)
nu_1.get_prob_on_chemical_species(atoms=atoms_QD, verbose=True)
analyze_dot(atoms=atoms_QD, psi=nu_1, verbose=True)

# Repeat these steps for the first excited state
print("First excited state (nu_2)")
nu_2.get_prob_on_orbitals(verbose=True)
nu_2.get_prob_on_chemical_species(atoms=atoms_QD, verbose=True)
analyze_dot(atoms=atoms_QD, psi=nu_2, verbose=True)

# Save modulus-square of projection of states on atoms to a .vtu file
prob_on_atoms_0 = nu_1.get_prob_on_atoms()
prob_on_atoms_1 = nu_2.get_prob_on_atoms()
prob_on_atoms_2 = atoms_QD.eigenfunctions[2].get_prob_on_atoms()
path_vtu = str(script_dir / "output" / "multiscale_TB_wf.vtu")
out_dict = {"Probability density (ground state)": prob_on_atoms_0,
    "Probability density (1st excited state)": prob_on_atoms_1,
    "Probability density (2nd excited state)": prob_on_atoms_2}
save_vtu(atoms=atoms_QD, out_dict=out_dict, path=path_vtu)

# Compute the matrix element of the spin-orbit coupling Hamiltonian between
# the valley states with opposite spins
H_soc = schrodinger_solver.get_hamiltonian_soc()
nu_2_up = nu_2.add_spin_to_orbital_state(spin=True, in_place=False)
nu_1_down = nu_1.add_spin_to_orbital_state(spin=False, in_place=False)
C = get_operator_element(psi_l=nu_2_up, op=H_soc, psi_r=nu_1_down)
print("Matrix element of SOC Hamiltonian: " + str(C/ct.e) + " eV")

# Compute the matrix element of the derivative of the potential with respect to
# barrier gate voltage between the valley states
qd = ["cap_qd", "well_qd", "buffer_qd"] # quantum dot regions
submesh = SubMesh(mesh, qd)
path_der_hdf5 = str(script_dir / "output" / "multiscale_der.hdf5")
der = io.load(path_der_hdf5, "der")
Dfg = get_operator_from_fem(atoms=atoms_QD, tb_orbs=nu_1.tb_orbs,
    tb_spin=nu_1.tb_spin, mesh=submesh, quantity=der)
D = get_operator_element(psi_l=nu_1, op=Dfg, psi_r=nu_2)
print("Matrix element of derivative of potential with respect to barrier " + \
    "gate voltage: " + str(D/ct.e) + " eV/V")

# Save parameters of quantum control problem to .txt file
params = np.array([valley_splitting, C, D])
path_params = str(script_dir / "output" / "multiscale_params.txt")
np.savetxt(path_params, params)