1. Introduction 
2. Method of Calculation 
   2.1 Simulation Domain 
   2.2 Geometry Construction 
   2.3 Strain Relaxation 
   2.4 Electronic Structure Calculation 
   2.5 Piezoelectricity 
   2.6 Optical Transition Strengths 
3. Multi-Particle States 
   3.1 Two Particle States from CI 
4. NEMO 3-D Flow Chart Diagram 
5. NEMO 3-D Scaling 


The theoretical knowledge of the electronic structure of nanoscale semiconductor devices is the first and most essential step towards the interpretation and the understanding of the experimental data and reliable device design at the nanometer scale. The following is a list of the modeling and simulation challenges in the design and analysis of realistically sized engineered nano-devices. 

(A) Full Three-Dimensional Atomistic Representation: The lack of spatial symmetry in the overall geometry of the nano-devices usually requires explicit three-dimensional representation. For example, Stranski-Krastanov growth techniques tend to produce self-assembled InGaAs/GaAs quantum dots [1-4] with some rotational symmetry, e.g. disks, truncated cones, domes, or pyramids [5]. These structures are generally not perfect geometric objects, since they are subject to interface inter-diffusion, and discretization on an atomic lattice. There is no such thing as a round disk on a crystal lattice! The underlying crystal symmetry imposes immediate restrictions on the realistic geometry and influences the quantum mechanics. Continuum methods such as effective mass [6] and k.p [7] typically ignore such crystal symmetry and atomistic resolution. 

The required simulation domain sizes containing millions of atoms prevent the usage of ab-initio methods. Empirical methods which eliminate enough unnecessary details of core electrons, but are finely tuned to describe the atomistically dependent behaviour of valence and conduction electrons, are needed. The current state-of-the-art leaves two choices: 1) pseudo-potentials [8] and 2) Tight Binding [9]. Both methods have their advantages and disadvantages. Pseudo-potentials use plane waves as a fundamental basis choice. Realistic nanostructures contain high frequency features such as alloy-disorder or hetero-interfaces. This means that the basis needs to be adjusted (by an expert) for every different device, which limit the potential impact for non-expert users. Numerical implementations of pseudo-potential calculations typically require a Fourier transform between real and momentum space which demand full matrix manipulations and full transposes. This typically requires high bandwidth communication capability (i.e. extremely expensive) parallel machines, which limit the practical dissemination of the software to end users with limited compute resources. Tight-binding is a local basis representation, which naturally deals with finite device sizes, alloy-disorder and hetero-interfaces and it results in very sparse matrices. The requirements of storage and processor communication are therefore minimal compared to pseudo-potentials and actual implementations perform extremely well on inexpensive clusters [9]. 

Tight-binding has the disadvantage that it is based on empirical fitting and some in the community continue to question the fundamental applicability of tight-binding. The NEMO 3-D team has spent a significant effort to expand and document the tight-binding capabilities with respect to handling of strain [10], electromagnetic fields [11], and Coulomb matrix elements [12] and fit them to well known and accepted bulk parameters [9]. With tight-binding, the NEMO team was able early on to match experimentally verified, high-bias current-voltage curves of resonant tunneling [13] that could not get modeled by ether effective mass (due to the lack of physics) or pseudo-potential methods (due to the lack of open boundary conditions). We continue to learn about the tight-binding method capabilities, and we are in the process of benchmarking it against more fundamental ab initio approaches and pseudo-potential approaches. Our current Si/Ge parametrisation is described in references [14]. Figure 1 depicts a range of phenomena that represent new challenges presented by new trends in Nanoelectronic and lays out the NEMO 3-D modeling agenda. 

Figure 1: NEMO 3-D modeling agenda: map electronic properties of individual atoms into realistic structures 
containing millions of atoms, computation of nanoscale quantum dots that maps into real applications. 

(B) Atomistic Strain: Strain that originates from the assembly of lattice-mismatched semiconductors strongly modifies the energy spectrum of the system. In case of the InGaAs/GaAs quantum dots, this mismatch is around 7% and leads to a strong long-range strain field within the extended neighborhood (typically ~ 25 nm) of each quantum dot [15]. Si/Ge core/shell structured Nanowire are another example of strain dominated atom arrangements [16]. Si quantum wells and SiGe quantum computing architectures rely on strain for state separation [17]. The strain can be atomistically inhomogeneous, involving not only biaxial components but also non-negligible shear components. Strain strongly influences the core and barrier material band structures, modifies the energy band gaps, and lifts the heavy hole-light hole degeneracy at the zone center. In the nanoscale regime, the classical harmonic linear/continuum elasticity model for strain is inadequate, and device simulations must include the fundamental quantum character of charge carriers and the long-distance atomistic strain effects with proper boundary conditions on equal footing [18]. 

(C) Piezoelectric Field: A variety of III-V materials such as GaAs, InAs, GaN, etc are piezoelectric. Any spatial non-symmetric distortion in nanostructures made of these materials will create piezoelectric fields, which will modify the electrostatic potential landscape. Recent spectroscopic analyses of self-assembled QDs demonstrate polarised transitions between confined hole and electron levels [5]. While the continuum models (effective mass or k.p) can reliably predict aspects of the single-particle energy states, they fail to capture the observed non-degeneracy and optical polarisation anisotropy of the excited energy states in the (001) plane. These methods fail because they use a confinement potential which is assumed to have only the shape symmetry of the nanostructure, and they ignore the underlying crystal symmetry. The experimentally measured symmetry is significantly lower than the assumed continuum symmetry because of (a) underlying crystalline symmetry, (b) atomistic strain relaxation and (c) piezoelectric field. For example, in the case of pyramid shaped quantum dots with square bases, continuum models treat the underlying material in C4 symmetry while the atomistic representation lowers the crystal symmetry to C2. The piezoelectric potential originating from the non-zero shear component of the strain field must be taken into account to properly model the associated symmetry breaking and the introduction of a global shift in the energy spectra of the system. 
 Back to Top



NEMO 3-D [9, 19] bridges the gap between the large size, classical semiconductor device models and the molecular level modeling. This package currently allows calculating single-particle electronic states and optical response of various semiconductor structures including bulk materials, quantum dots, quantum wires, quantum wells and nanocrystals. NEMO 3-D includes spin in its fundamental atomistic tight binding representation. Spin is therefore not added in as an afterthought into the theory, but spin-spin interactions are naturally included in the Hamiltonian. Effects of interaction with external electromagnetic fields are also included [9, 11]. A schematic view of InAs quantum dot embedded in a GaAs barrier material the sample is presented in Figure 2. The quantum dot is positioned on a 0.6 nm thick wetting layer (dark region). The simulation of strain is carried out in the large computational box Dstrain, while the electronic structure computation is restricted to the smaller domain Delec. Strain is long-ranged and penetrates around 25 nm into the dot substrate thus stressing the need for using large substrate thickness in the simulations. NEMO 3-D enables the computation of strain and electronic structure in an atomistic basis for over 64 and 52 million atoms, corresponding to volumes of (110nm)3 and (101nm)3, respectively. These volumes can be spread out arbitrarily over any closed geometry. For example, if a thin layer of 15 nm height is considered, the corresponding widths in the x-y plane correspond to 298 nm for strain calculations and 262 nm for electronic structure calculations. No other atomistic tool can currently handle such volumes needed for realistic device simulations. NEMO 3-D runs on serial and parallel platforms, local cluster computers as well as the NSF TeraGrid. 

Figure 2. Simulated dome shaped InAs quantum dot buried in GaAs. Two simulation domains are shown, 
Delec: central smaller domain for electronic structure calculation, and Dstrain: outer larger domain for strain calculation.
 In the figure: "s" is the substrate height, "c" is the cap layer thickness, "h" is the dot height, "d" is the dot diameter.
Back to Top

The first part is the geometry constructor, whose purpose is to represent the treated nanostructure in atomistic detail in the memory of the computer. Each atom is assigned three single-precision numbers representing its coordinates, stored is also its type (atomic number in short integer), information whether the atom is on the surface or in the interior of the sample (important later on in electronic calculations), what kind of computation it will take part of (strain only or strain and electronic), and what its nearest neighbor relation in a unit cell is. The arrays holding this structural information are initialized for all atoms on all CPUs, i.e., the complete information on the structure are available on each CPU. By default most of this information can be stored in short integer arrays or as single bit arrays, which does not require significant memory. This serial memory allocation of the atom positions, however, becomes significant for very large systems which must be treated in parallel.
Back to Top

The materials making up the QD nanostructure may differ in their lattice constants; for the InAs/GaAs system this difference is of the order of 7%. This lattice mismatch leads to the appearance of strain: atoms throughout the sample are displaced from their bulk positions. Knowledge of equilibrium atomic positions is crucial for the subsequent calculation of QD's electronic properties, which makes the computation of strain a necessary step in realistic simulations of these nanostructures. NEMO 3-D computes strain field using an atomistic valence force field (VFF) method [20] with the Keating Potential. In this approach, the total elastic energy of the sample is computed as a sum of bond-stretching and bond-bending contributions from each atom. The local strain energy at atom i is given by a phenomenological formula: 

where the sum is carried out over the "n" nearest neighbours "j" of atom "i", and dij, Rij are the bulk and actual (distorted) distances between neighbour atoms, respectively, and α and β are empirical material-dependent elastic parameters. The equilibrium atomic positions are found by minimising the total elastic energy of the system. Several other strain potentials [18] are also implemented in NEMO 3-D. While they modify some of the strain details they roughly have the same computational efficiency.
Back to Top

The single-particle energies and wave functions are calculated using an empirical nearest-neighbour tight-binding model. The underlying idea of this approach is the selection of a basis consisting of atomic orbitals (such as s, p, d, and s*) centred on each atom. These orbitals are further treated as a basis set for the Hamiltonian, which assumes the following form: 

where c+i,v (ci,v) is the creation (annihilation) operator of an electron on the orbital localised on atom "i". In the above equation, the first term describes the onsite orbital terms, found on the diagonal of the Hamiltonian matrix. The second term describes coupling between different orbitals localised on the same atom (only the spin-orbit coupling between p-orbitals), and the third term describes coupling between different orbitals on different atoms. The restriction in the summation of the last term is that the atoms "i" and "j" be nearest neighbours. 

The characteristic parameters ε and t are treated as empirical fitting parameters for each constituent material and bond type. They are usually expressed in terms of energy constants of σ and π bonds between the atomic orbitals. For example, for a simple cubic lattice, the interaction between the s-orbital localised on the atom "i" at origin and the orbital px localised on the atom "j" with coordinate d=ax with respect to the atom "i" would simply be expressed as tij=Vspσ. Most of the systems under consideration, however, crystallise in the zinc-blende lattice, which means that the distance between the nearest neighbours is described by a 3-D vector dij=lx + my + nz , with l, m, n being the directional cosines. These cosines rescale the interaction constants, so that the element describing the interaction of the orbitals s and px is tij=lVspσ. The parameterization of all bonds using analytical forms of directional cosines for various tight-binding models is given in Ref 21. NEMO 3-D provides the user with choices of the sp3d5s*, sp3s*, and single s-orbital models with and without spin, in Zincblende, Wurtzite, and simple cubic lattices. 

Additional complications arise in strained structures, where the atomic positions deviate from the ideal (bulk) crystal lattice [22]. The presence of strain leads to distortions not only of bond directions, but also bond lengths. In this case, the discussed interaction constant tij=l'Vspσ(d/do)η, where the new directional cosine l' can be obtained analytically from the relaxed atom positions, but the bond-stretch exponent η needs to be fitted to available data. The energy constants parameterising the on-site interaction change as well due to bond renormalisation [9]. 

The 20-band nearest-neighbour tight-binding model is thus parameterised by 34 energy constants and 33 strain parameters, which need to be established by fitting the computed electronic properties of materials to those measured experimentally. This is done by considering bulk semiconductor crystals (such as GaAs or InAs) under strain. The summation in the Hamiltonian for these systems is done over the primitive crystallographic unit cell only. The model makes it possible to compute the band structure of the semiconductor throughout the entire Brillion zone. For the purpose of the fitting procedure, however, only the band energies and effective masses at high symmetry points and along the Δ line from Γ to X are targeted, and the tight-binding parameters are adjusted until a set of values closely reproducing these target values is found. Search for optimal parameterization is done using a genetic algorithm, described in detail in Refs 9. Once it is known for each material constituting the QD, a full atomistic calculation of the single-particle energy spectrum is carried out on samples composed of millions of atoms. No further material properties are adjusted for the nanostructure, once they are defined as basic bulk material properties.
Back to Top

Piezoelectricity is defined as the generation of the electrical polarisation by the application of stress to a crystal lacking a centre of symmetry [23]. The Zinc-blende structure is one of the simplest examples of such a lattice, and the strength of the resulting polarisation is described by one parameter alone, e14 for the linear case (resulting in a polarisation P1), and three parameters, B114, B124, and B156, for the quadratic case (resulting in a polarization P2) (see the table 1) [24]. Their relation to the strain tensor field is given by: 

The resulting piezoelectric potential is obtained by solving the Poisson's equation, taking into account the material dependence of the static dielectric constant εs(r). The value of the dielectric constant for vacuum ε0 is 8.85x10-12F/m. For InGaAs and GaAs materials, we used the relative dielectric constant values of 14.00ε0 and 12.84ε0, respectively. 

Table 1: Polarisation constants for calculation of piezoelectric potential from reference [24]. 
The values for In0.5Ga0.5As are obtained by linear interpolation between GaAs and InAs. 

How can such a piezoelectric polarisation be included in an electronic structure calculation? A self-consistent calculation of the electronic structure of a deformed solid naturally includes the field generated by piezoelectric displacements. Of course, such a self-consistent evaluation requires the inclusion of all occupied states and their response to strain effects. Therefore such an approach is limited to very small systems. For larger systems it is often impractical to calculate all occupied states, and one has to concentrate oneself to only a few states in the energy range of interest. In such cases the calculation is no longer self-consistent and the piezoelectricity does not arise naturally. Instead it has to be introduced as an external potential Vp in the Schrodinger equation. The resulting electrostatic potential is included in the TB model as a site-diagonal potential energy Vp = -eφpiezo.
Back to Top


As a first step, we consider a system with a Hamiltonian {Ho} and a time-dependent perturbation {Hp} switched on at t = 0+. Then, the total Hamiltonian {HT} becomes as shown in (eq. 7). 

Now, let's suppose that an electron or hole was sitting in energy eigenstates ΦK. Then, we can express the electron/hole states in the conduction/valence band for t > 0 as shown (eq. 8) and Fermi's Golden rule says the evaluation of corresponding optical transition rate from state K to N (Absorption) or N to K (Emission) is done with the numerical expression shown in (eq. 9), where we assume (Energy at state N) > (Energy at state K). 

If we assume the equality and symmetry of perturbation matrices {H+ , H-}, the rate for absorption and emission process between a fixed set of {initial state K, final states N} becomes equal, which can be evaluated in further detail, using the numerical expression as shown in (eq. 10). Please note that we used subscripts {i, f} to represent initial/final states instead of {K, N}. 

The momentum matrix component "M" is defined as shown in (eq. 11) and it is quite clear that the momentum matrix components fully determines the transition rate as long as we assume a completely filled initial state and empty final state. Throughout our model, we focus on the evaluation of these momentum matrix components to estimate inter-band transition rates for quantum dot systems.
Back to Top

The numerical details for the momentum matrix introduced in the previous section works for the general bulk systems. The expressions, therefore, still need further modifications because quantum dots are 3-D finite (closed) systems without any periodicities. In this section, we discuss and show the details how the numerical expression in the previous section should be modified. 

For simplicity, let's consider 1-Dimensional finite structure with length "L". Then, only a finite number of "k" points satisfy the closed boundary condition in the 1st Brillouin-zone of a unit-cell. If we restrict our interest to the 1st Brillouin-zone of a super-cell with length "L", the only available "k" point will be zero (Γ point) as described in (eq. 12). 


As we talked in the previous section, the inter-band transition rate between ground states should increase with increasing density of electrons (at CB ground state) and holes (at VB ground state) at the same "spatial" point because the recombination rate is directly proportional to the spatial density of carriers. In this work, therefore, we compute a spatial overlap of wave functions at CB and VB ground state to confirm this basic physical phenomenon, as well as to provide a qualitative explanation for the effect of strain and corresponding local polarisation on the optical property of quantum dots. This spatial overlap can be evaluated with (eq. 15), which is the first moment of a normal vector of the incident light with respect to wave functions at CB and VB ground state. 

Back to Top

Equation 15 gives optical absorption strength between state "f" and "i" in the "n" direction. Equation 15 gives first moment of spatial overlap between state "f" and "i" in "n" direction. Below, a step by step procedure is presented to solve these equations between two single particle states for the "x" direction polarization of light. Same procedure can be applied to calculate these equations in any general direction "n". The initial state is in the conduction band and final state is in the valence band: i=c, f=v. 

Steps to calculate optical matrix element in "x" direction:

1. Solve a single particle Hamiltonian H for the single particle Eigen-energies and states for the conduction and valence band. Let e1 and h1 are the energies for the ground electron and hole states. Let |φc> and |φv> be the corresponding ground electron and hole states. 
2. Find the Kramer degenerate states for the ground electron and hole states: |φc'> ,|φv'>. These states will have the same energy value and the opposite spin.
3. Calculate: x|φc> and x|φc'>. The method to calculate these is: Each atom has a real space coordinate {x, y, z}. Now the vectors |φc> and |φc'> will have twenty values for each atomic sight because of the twenty band orbital basis. These twenty entries for one atomic sight will be multiplied with the corresponding "x" value of the atomic coordinate.
4. Calculate:

v|x|φc>, <φv|x|φc'> <φv'|x|φc>, <φv'|x|φc'>

5. Equation 16 =

v|x|φc> + <φv|x|φc'> + <φv'|x|φc> + <φv'|x|φc'>

6. Calculate:

M(c,v) = <φv|[x, H]|φc> = i(e1-h1)<φv|x|φc>
M(c',v) = <φv|[x, H]|φc'> = i(e1-h1)<φv|x|φc'>
M(c,v') = <φv'|[x, H]|φc> = i(e1-h1)<φv'|x|φc>
M(c',v') = <φv'|[x, H]|φc'> = i(e1-h1)<φv'|x|φc'>

Optical Transition Rate between e1 and h1 : 

Equation 15 : 

T(e1-h1) = |M(c,v)|2 + |M(c',v)|2 + |M(c,v')|2 + |M(c',v')|2

Back to Top


The fundamental picture of the quantum dot electronic structure is represented in terms of the single particle energy states i.e. one electron is lifted from the valence band and goes into an empty conduction band. The valence band with one electron missing is represented as having one hole. The holes form new quasiparticles with the charge and spin opposite to that of the electron jumping from the valence band to the conduction band. This system of one electron and one hole is represented by the Hamiltonian: 

H = He + Hh          (17)

Here He/h are the single particle Hamiltonians of the electron/hole. However, due to the electron-electron interaction or in this case electron-hole interaction, this system of one electron and one hole cannot be simply described as the sum of the electron Hamiltonian and the hole Hamiltonian. Rather, the electron and hole are attached to each other through the Coulomb interaction and hence form a new quasiparticle called exciton. The energy of this exciton can be determined by the Hamiltonian, 

H = He + Hh + Hcoul          (18)

where Hcoul is the coulomb interaction between the electron and the hole and it lowers the total energy of an exciton with respect to the total energy of a free electron plus a free hole. The difference between the exciton energy and the free particle energy is called the "binding energy of the exciton". 

The strongest or the most important part of the Coulomb interaction is called the direct or the mean field Coulomb interaction. Let Ce-h be the Coulomb interaction between the electron and the hole charge carriers given by: 

In the single particle formulation, the exciton energy will be given by: 

E(e-h) = Ee + Eh + Ce-h           (20)

Apart from the exciton, also called a neutral exciton (X), comprised of one electron and hole, a quantum dot can also be comprised of many electrons and holes giving rise to the higher order excitons called the charged excitons. The charged excitons are also called the trions (X+, X-) comprising of one electron + two hole (X+), or two electrons + one hole (X-), and the biexcitons (XX) comprising of two electrons + two holes. The energies of these charged excitons can be calculated as follows: 

Note that the Ce-e and Ch-h are the repulsive interactions i.e. Ce-e, Ch-h > 0, whereas the Ce-h is the attractive interaction i.e. Ce-h < 0. Depending upon the (positive/negative) magnitude of these repulsive and attractive interactions, the exciton complex will be of either a bonding or an anti-bonding type. 

So far the spin part of the electron/hole has been ignored. The electrons and the holes have spin of the magnitude 1/2 and they follow Pauli's exclusion principle. This gives rise to another kind of interaction energy called the exchange interaction energy. The exchange interaction energy is given by: 

The neutral exciton energy arising from an electron and a hole interacting by the Coulomb and Exchange interactions (Eq. 19, 24) is described by the first order perturbation theory as: 

E(e-h) = Ee + Eh + Ce-h - Ke-hδS,0          (25)

where δS,0 is 1 for the triplet states and 0 for the singlet states.
Back to Top

The above picture described by the equation 20 is only based on the first order perturbation theory. It comprised of only one electron and one hole particle in the system. Normally we can have multiple electrons in the conduction band leaving behind multiple holes in the valence band. In such system, it is required to solve the Schrodinger equation taking into account the Coulomb and the Exchange interactions between all of the pairs of the electrons and the holes. This is called an "N" electron system. A more detailed analysis can be achieved by employing the configuration interaction (CI) method, which takes into account the direct Coulomb, the Exchange, and the large parts for self-consistency effects and correlation. However, the calculation of the CI method for the "N" electron system known as the complete or full CI method is very cumbersome and its computation becomes extremely time consuming. So we restrict ourselves to only two electron picture where the conduction band has two electrons and the valence band has two holes. Moreover, it is assumed that the system has only two possible energy states (each is spin degenerate). These four energy states are labeled as E↑, E↓, H↑, and H↓, where ↑/↓ represent the up/down spin of the electron in a state and E/H mark it as the conduction/valence band state. In this two particle picture, these four states will make six possible configurations for the two electrons as shown in the figure 3: 

Figure 3: Six two particle states are shown for a two electron and two single particle energy levels configuration. 
The black dots mark the position of an electron and the empty squares mark the position of a hole. 

Let us suppose that the single particle states are represented by the equations 26 to 29: 

E↑ = Χ1 = Ψ(r1).α(s1)           (26)
E↓ = Χ2 = Ψ(r1).β(s1)          (27)
H↑ = Χ3 = Ψ(r2).α(s2)           (28)
H↓ = Χ4 = Ψ(r2).β(s2)          (29)

where Ψ represent the spatial wave function and α/β represent the up/down spin of the state. 

Now the two particle states from the Slater determinants are given in the equations 30 to 35: 

State 1: Ψ(r2,r2) = 1/√(2)(Χ3 . Χ4 - Χ4 . Χ3)          (30)
State 2: Ψ(r1,r2) = 1/√(2)(Χ1 . Χ3 - Χ3 . Χ1)          (31)
State 3: Ψ(r1,r2) = 1/√(2)(Χ1 . Χ4 - Χ4 . Χ1)           (32)
State 4: Ψ(r1,r2) = 1/√(2)(Χ3 . Χ2 - Χ2 . Χ3)          (33)
State 5: Ψ(r1,r2) = 1/√(2)(Χ2 . Χ4 - Χ4 . Χ2)           (34)
State 6: Ψ(r1,r1) = 1/√(2)(Χ1 . Χ2 - Χ2 . Χ1)           (35)

These two particle states can be computed by finding the Eigen-values of a two particle Hamiltonian given by the equation 31: 

H(ij,kl) = <Ψ(i,j) | H1 + H2 + 1/r12 | Ψ(k,l)>          (36)

where H1 and H2 are the single particle Hamiltonians for the electron '1' and electron '2'. The state '1' in the figure 3 is called the 'singlet ground state' as it has the lowest energy and both of the electrons are in the valence band. The state '6' in the figure 3 is called the 'singlet excited state' because it has the highest energy and both of the electrons are in the conduction band. The states 2, 3, 4, and 5 in the figure 3 are the excited states where only one electron has jumped up from the valence band to the conduction band. If the energy of E ↑ ↓ is ε1, the energy of H ↑ ↓ is ε2, and the coulomb and exchange interactions of the electrons are given by the equations 19 and 24, then the energies of the two particle states will be given by: 

State 1: 2ε2 + C 
State 2: ε1 + ε2 + C + K 
State 3: ε1 + ε2 + C 
State 4: ε1 + ε2 + C 
State 5: ε1 + ε2 + C + K 
State 6: 2ε1 + C 

 Back to Top


The flow chart diagram of the NEMO 3-D simulator is shown below which combines all the components described in the previous sections. First, the geometry information extracted from the experimental data is provided in an input deck. The single particle Hamiltonian is solved including the piezoelectric potentials and the electric/magnet fields for the single particle energy states and wave-functions. Finally, the CI calculations based on the Hartree-Fock formulation are done to calculate the few particle excitonic states. Inter-band optical transition intensities between various excitonic states are calculated using Fermi's golden rule from the squared magnitude of the optical matrix elements summed over the degenerate spin states (see section 2.6). Further details about the NEMO 3-D simulator can be found here

Figure 4: Flow chart diagram of NEMO 3-D simulator. 
Back to TOP


The complexity and generality of physical models in NEMO 3-D can place high demands on computational resources. For example, in the 0-band electronic calculation the discrete Hamiltonian matrix is of order 20 times the number of atoms. Thus, in a computation with 20 million atoms, the matrix is of order 400 million. Computations of that size can be handled because of the parallelized design of the package. NEMO 3-D is implemented in ANSI C, C++ with MPI used for message-passing, which ensures its portability to all major high-performance computing platforms, and allows for an efficient use of distributed memory and parallel execution mechanisms. 

Although the strain and electronic parts of the computation are algorithmically different, the key element in both is the sparse matrix-vector multiplication. This allows the use of the same memory distribution model in both phases. The computational domain is divided into slabs along one dimension. All atoms from the same slab are assigned to a single CPU, so if all nearest neighbors of an atom belong to its slab, no inter-CPU communication is necessary. The inter-atomic couplings are then fully contained in one of the diagonal blocks of the matrix. On the other hand, if an atom is positioned on the interface between slabs, it will couple to atoms belonging both to its own and the neighbouring slab. This coupling is described by the off-diagonal blocks of the matrix. Its proper handling requires inter-CPU communication. However, due to the first-nearest-neighbor character of the strain and electronic models, the messages need to be passed only between pairs of CPUs corresponding to adjacent domains - even if the slabs are one atomic layer thick. Full duplex communication patterns are implemented such that all inter-processor communications can be performed in 2 steps [9]. 

Figure 5 shows the memory requirements for the dominant phase of the code (electronic structure calculations). It shows how the number of atoms that can be treated grows as a function of the number of CPUs, for a fixed amount of memory per CPU. The number of atoms can be intuitively characterized by the length of one side of a cube that would contain that many atoms. This length is shown in Figure 5, on the vertical axis on the right side of each plot. This figure shows that the number of atoms that can be treated in NEMO 3-D continues to grow for larger CPU counts. The strain calculations have so far never been memory limited. NEMO 3-D is typically size limited in the electronic structure calculation. 

Figure 5: Number of atoms that can be treated, as a function of the number of CPUs for different amounts of memory per CPU for the electronic structure calculation. 
The vertical axis on the right side of each plot gives the equivalent length in nm of one side of the cube that would contain the given number of atoms. 

To investigate the performance of NEMO 3-D package, computation was performed in a single dome shaped InAs quantum dot nanostructure embedded in a GaAs barrier material as shown in Figure 1. The HPC platform used in the performance studies are shown in Table 2. These include a Linux clusters at the Rosen Centre for Advanced Computing (RCAC) at Purdue with Intel processors (dual core Woodcrest). The other five platforms are a BlueGene at the Rensselaer Polytechnic Institute (RPI), the Cray XT3 at the Pittsburgh Supercomputing Center (PSC), the Cray XT3/4 at ORNL, JS21 at Indiana University, and a Woodcrest machine at NCSA. Table 2 provides the relevant machine details. These platforms have proprietary interconnects, that are higher performance than Gigabit Ethernet (GigE) for the three Linux clusters at Purdue. In the following, the terms processors and cores are used interchangeably. 

Table 2: Specifications for the HPC platforms used in the performance comparisons. 

In addition to the performance for the benchmark cases end-to-end runs on the PU/Woodcrest cluster are carried out next (Figure 6). This involves iterating to convergence and computing the Eigen states in the desired range (4 conduction band and 4 valence band states). For each problem size, measured in millions of atoms, the end-to-end cases were run to completion, for one choice of number of cores. 

Figure 6: Wall clock time vs. number of atoms for end-to-end computations of the electronic structure of a quantum dot, 
for various numbers of cores on the PU/Woodcrest cluster. Listed next to the number of cores are the CPU hours/Million 
of atoms needed in the simulation. 

Back to Top

[1] Petroff PM (2003) Single Quantum Dots: Fundamentals, Applications, and New Concepts. Springer, Berlin 
[2] Michler P, Kiraz A, Becher C, Schoenfeld W, Petroff P, Zhang L, Hu E, Imamoglu A (2000) A Quantum Dot Single-Photon Turnstile Device. Science 290:2282-2285 
[3] Reed M (1993) Quantum Dots. Scientific American 268:118 
[4] Reed M, Randall J, Aggarwal R, Matyi R, Moore T, and Wetsel A (1988) Observation of discrete electronic states in a zero-dimensional semiconductor nanostructure. Phys Rev Lett 60:535 
[5] Bester G, Zunger A (2005) Cylindrically shaped zinc-blende semiconductor quantum dots do not have cylindrical symmetry: Atomistic symmetry, atomic relaxation, and piezoelectric effects. Physical Review B, 71:045318. Also see references therein. 
[6] Pryor C, Kim J, Wang L, Williamson A, Zunger A (1998) Comparison of two methods for describing the strain profiles in quantum dots. J Apl Phys 83:2548 
[7] Stier O, Grundmann M, and Bimberg D (1999) Electronic and optical properties of strained quantum dots modeled by 8-band k.p theory. Phys Rev B, 59:5688 
[8] Calder`on, M J, Koiler B, Hu X, Das Sarma S (2006) Phys Rev Lett 96:096802 
[9] Klimeck G, Oyafuso F, Boykin T, Bowen R, Allman P (2002) Development of a Nanoelectronic 3-D (NEMO 3-D) Simulator for Multimillion Atom Simulations and Its Application to Alloyed Quantum Dots. Computer Modeling in Engineering and Science 3:601 
[10] Boykin T, Klimeck G, Bowen R, Oyafuso F (2002) Diagonal parameter shifts due to nearest-neighbor displacements in empirical tight-binding theory. Phys Rev B 66:125207 
[11] Boykin T, Bowen R, Klimeck G (2001) Electromagnetic coupling and gauge invariance in the empirical tight-binding method. Phys Review B 63:245314 
[12] Lee S, Kim J, Jonsson L, Wilkins J, Bryant G, Klimeck G (2002) Many-body levels of multiply charged and laser-excited InAs nanocrystals modeled by empirical tight binding. Phys Rev B 66:235307 
[13] Klimeck G, Boykin T , Chris R, Lake R, Blanks D, Moise T, Kao Y, Frensley W (1997) Quantitative Simulation of Strained InP-Based Resonant Tunneling Diodes. Proceedings of the 1997 55th IEEE Device Research Conference Digest:92 
[14] Boykin T, Kharche N, Klimeck G (2007) Brillouin zone unfolding of perfect supercells composed of non-equivalent primitive cells. Phys Rev B 76:035310 
[15] Ahmed S, Usman M, Heitzinger C, Rahman R, Schliwa A, Klimeck G (2007) Atomistic Simulation of Non-Degeneracy and Optical Polarization Anisotropy in Zincblende Quantum Dots. The 2nd Annual IEEE International Conference on Nano/Micro Engineered and Molecular Systems (IEEE-NEMS), Bangkok, Thailand 
[16] Liang G, Xiang J, Kharche N, Klimeck G, Lieber C, Lundstrom M (2006) Performance Analysis of a Ge/Si Core/Shell Nanowire Field Effect Transistor. cond-mat 0611226 
[17] Eriksson M, Friesen M, Coppersmith S, Joynt R, Klein L, Slinker K, Tahan C, Mooney P, Chu J, Koester S (2004) Spin-based quantum dot quantum computing in Silicon. Quantum Information Processing 3:133 
[18] Lazarenkova O, Allmen P, Oyafuso F, Lee S, Klimeck G (2004) Effect of anharmonicity of the strain energy on band offsets in semiconductor nanostructures. Appl Phys Lett 85:4193 
[19] Korkusinski M, Klimeck G (2006) Atomistic simulations of long-range strain and spatial asymmetry molecular states of seven quantum dots. Journal of Physics Conference Series 38:75-78 
[20] Keating P (1966) Effect of Invariance Requirements on the Elastic Strain Energy of Crystals with Application to the Diamond Structure. Phys Rev :145 
[21] Srijit Goswami, K. A. Slinker, Mark Friesen, L. M. McGuire, J. L. Truitt, Charles Tahan, L. J. Klein, J. O. Chu, P. M. Mooney, D. W. van der Weide, Robert Joynt, S. N. Coppersmith, and Mark A. Eriksson, Nature Physics, 3, 41 (2007). 
[22] Jancu J, Scholz R, Beltram F, Bassani F (1998) Empirical spds* tight-binding calculation for cubic semiconductors: General method and material parameters. Phys Rev B 57:6493 
[23] W.F. Cady, Piezoelectricity (McGraw-Hill, New York, 1946) 
[24] Andrei Schliwa, M.Winkelnkemper, and D. Bimberg, Pyhs. Rev B 76, 205324 (2007)
Back to Top
This page consists of material copied from: 
1. M. Usman, "Multi-million Atom Electronic Structure Calculations for Quantum Dots", Chapter 2 of PHD Thesis 2010, Purdue University West Lafayette IN, USA.
2. S. Ahmed, N. Kharche, R. Rahman, M. Usman et al., "Multimillion Atom Simulations with NEMO 3-D", published in Encyclopedia of Complexity and Systems Science, 2009, LXXX, p. 10370 illus. 4300, 11 Volumes, ISBN: 978-0-387-75888-6, Editor-in-chief: Meyers, Robert A., (Pre-print version here) 
Back to Top
Subpages (1): NEMO3D