# Petascale Electronic Structure Calculations for Nanostructures

**Lin-Wang Wang, Zhengji Zhao and Juan Meza**

## OVERVIEW

Density functional theory (DFT) is the most widely used *ab initio*
method in material simulations accounting for 75% of the NERSC
allocation time in the materials science category. DFT can be used to
calculate the electronic structure, the charge density, the total energy
and the atomic forces of a material system, and with the advance of new
algorithms, DFT can now be used to study thousand-atom systems. But
there are many problems that either require much larger systems (more
than 100,000 atoms), or many total energy calculation steps (molecular
dynamics or atomic relaxations). Examples include: grain boundary,
dislocation energies and atomic structures, impurity transport and
clustering in semiconductors, nanostructure growth, electronic
structures of nanostructures and their internal electric fields. However
due to the O(N^3^) scaling of conventional DFT algorithms (as
implemented in codes like Qbox, Paratec, Petot) these problems are even
beyond the reach of planned petascale computers. In addition, as
proposed petascale computers might have millions of processors, new
computational paradigms and algorithms are needed to solve these
problems. For a large materials science problem, a natural approach to
achieve this goal is by a divide-and-conquer approach: spatially divide
the system into many small pieces and solve each piece by a small local
group of processors. This solves the O(N) scaling and the
parallelization problem at the same time. However, the challenge of this
approach is in decomposing the system into small pieces and patching
them up again. We present a linear scaling 3 dimensional fragment
(LS3DF) method that uses a novel decomposition/patching scheme that
cancels out the artificial boundary effects of the spatial division. As
a result, the LS3DF results are essentially the same as the original
full system DFT results with errors smaller than chemical accuracy, and
the resulting computational times are thousands of times smaller than
conventional DFT methods, making calculations with 100,000 atom systems
feasible. The LS3DF method is applicable to insulator and semiconductor
systems, which covers a current gap in DOE's material science code
portfolio for *ab initio* ultrascale simulations. We will use it here to
solve the internal electric field problems for composite nanostructures.

## SCIENTIFIC PROBLEM

Nanostructures such as quantum dots and wires, composite quantum rods
and core/shell structures have been proposed for electronic devices or
optical devices like solar cells. Understanding the electronic
structures of such systems and the corresponding carrier dynamics is
essential to the successful design and deployment of such devices.
Despite more than a decade of research, some critical issues of the
electronic structure of moderately complex nanostructures are still
poorly understood. One issue is the internal electric field in a
composite colloidal nanostructure and its consequences on the electron
wavefunctions. It is well known that there are strong internal electric
fields in some of the bulk semiconductor heterostructures, like the
InN/GaN superlattice. These electric fields could be caused by surface
and interface dipoles, total dipole of the nanostructure, piezoelectric
effects, surface trapped charges and charged dopants. They induce strong
spatial localizations of the wavefunctions, thus cause different
electron-hole recombination rates, charge transports, and nonlinear
optical properties, all important to the performance of the
nanostructure electro-optical devices. Unfortunately, the continuous
model used in conventional device simulations can no longer be used for
these nanostructures due to the atomic natures of the charge, dopant and
geometry, the high order effects of different phenomenon, and the change
of dielectric functions, etc. Thus, what is needed here is a direct
atomistic *ab initio* selfconsistent calculation for the charge density
and the electric field, and the corresponding atomic relaxation for the
nanosystem. Since the atomic number N of such composite nanostructures
easily exceeds 10,000 or even 100,000 atoms, the traditional O(N^3^)
scaling DFT method cannot be used. Partly because of this, the internal
electric field problem remains one of the most outstanding unsolved
problems in colloidal nanoscience. For example, we don't even know
whether there is a large internal electric field in a simple quantum dot
consisted with dipolar semiconductors, and if there is one, what is the
cause of such internal electric field.

Here, we propose to use our newly developed LS3DF method to calculate
10,000--100,000 atom composite nanosystems using the local density
approximation (LDA) of the DFT. We will investigate nanostructures with
different geometries and heterostructure composites: A/B nanorod, A/B
core/shell nanowire, A/B/A dumbbell structures. We will also study 1)
the effect of different surface passivation and surface termination
layers; 2) the effect of a charge trapped in a surface dangling bond,
and 3) the effect of a single dopant in a nanostructure, and test the
bulk concepts of p-type and n-type semiconductors, and p-n conjunctions
in colloidal nanosystems. Our theoretical calculations could help
experimentalists design better solar cells using nanostructures, which
has been identified as one of the possible interruptive technologies
that can revolutionize the solar cell field.

## ALGORITHM

While planewave LDA codes, like Qbox, have demonstrated the capability
of using hundreds of thousands of processors on the BlueGene/L computer
\[1\], the *O(N*^3^) scaling for its floating point operation (and
*O(N*^2^) communication) may not necessarily be the most efficient way
to solve a given science problem. To reduce the time to solution, one
would ideally like to use a linear scaling method *ab initio* method
\[2\]. Most of these methods use localized orbitals, and minimize the
total energy as a function of these orbitals. Unfortunately the use of
localized orbitals can introduce extraneous local minima in the total
energy functional, which makes the total energy minimization more
difficult. On the computational side, it is difficult to effectively use
thousands of processors because localized orbitals can have strong
overlaps that make parallelization a nontrivial task. As a result of
these challenges and in spite of a decade of intense research, the
application of these types of linear scaling methods is still quite
limited. Another *O(N*) approach, the LSMS method \[3\], has been shown
to scale to thousands of processors. However this method can only be
applied to metals and has been mainly used to study metallic alloys and
magnetic systems. The LSMS method could not be used to study, for
example, the class of semiconductor nanostructures that we propose here.

Our LS3DF is based on the observation that the total energy of a system
can be decomposed into the quantum mechanical part (the wavefunction
kinetic energy and the exchange correlation energy), and the classical
electrostatic part. While the electrostatic energy (Coulomb energy) is
long range, the quantum mechanical energy is local in nature (short
sighted). While the long range Coulomb interaction can be solved
efficiently by the Poisson equation even for million-atom systems, the
quantum mechanical part is the most difficult one to be solved. For
this, we will take advantage of its locality by using the aforementioned
spatial decomposition divide-and-conquer method. While there are
previous methods \[4,5\] based on this divide-and-conquer concept, they
all rely on positive spatial partition functions to divide and patch the
spaces. There are intrinsic difficulties to use this positive partition
function technique, especially for dividing the kinetic energies. In
contrast, our novel division-patching method avoids these problems, thus
result in a much more accurate algorithm (with the accuracy essential
the same as the original full system LDA method).

Our LS3DF spatial division-patching technique is inspired by the
fragment molecular orbital (FMO) method proposed by Kitaura *et al*
\[6,7\] and combined with the ideas from our own charge patching
techniques \[8\]. FMO is used for organic chain like molecules. In FMO,
the long chain molecule is chopped into fragment pieces. Then each piece
and pairs of nearby pieces are calculated. The electron charge density
is then added up from all the pieces and their pairs, with positive sign
for the pair and negative sign for the pieces itself. The usage of pairs
and negative signs are innovative as this allows the calculation of the
energy of the artificial boundaries, which can subsequently be
subtracted from the total energy and charge density summation. Our LS3DF
method extends this technique to arbitrary 3 dimensional systems.
Instead of using pairs of pieces, we divide the system using overlapping
regions (pieces, fragments). More specifically, our division scheme is
illustrated in Fig.1 for a 2 dimensional system for clarity. Here, a
supercell is divided into m1x m2 grid points. From each grid point
corner (i~1~,i~2~), we can defined four pieces with dimension: 1x1, 1x2,
2x1, 2x2. Note that, they are overlapping pieces. Now, after all the
pieces at all the (i~1~,i~2~) corners are calculated, the total charge
density of the whole system can be patched together
as:$\rho_{\text{tot}}(r) = \sum_{(i1,i2),D}^{}{\text{sign}_{D}\rho_{(i1,i2),D}}(r)$,
here D denotes the dimension 1x1, 1x2, 2x1, 2x2, and the sign~D~ is +
for 1x1 and 2x2, and -- for 1x2, 2x1. The total energy (especially the
kinetic energy part) can be expressed in similar fashion using the
wavefunctions of each piece, although the electron-electron Coulomb
interaction is expressed based on the total charge density ρ~tot~(r). To
make sense out of the above formula, we can check each point inside a
piece (A point in Fig.1). Note that, each spatial point will be included
in 3^2^ pieces: 4 "2x2" pieces, 2 "2x1" pieces, 2 "1x2" pieces, and 1
"1x1" piece. After the above +/- cancellations, it will be covered by
only one piece, that is what we need. We can also check for each
boundary point. A boundary can be defined with a direction (i.e,
boundary from A to B is different than boundary from B to A, we has used
an arrow in Fig.1 to represent a directional boundary). A given
directional boundary is covered by 6 pieces, with equal numbers of
positive and negative signs. Since all these pieces have the same
(directional) boundary at that point, and given the near-sightedness,
their charge density will be the same near that point. As a result, the
boundary effect will be cancelled out. The same is true for the corner
effects. This division scheme can be extended to 3 dimension, where at
each corner point (i~1~,i~2~,i~3~), there will be eight pieces: 1x1x1,
1x1x2, 1x2x1, 2x1x1, 1x2x2, 2x1x2, 2x2x1, 2x2x2. In this case, each
spatial point will be covered by 3^3^ pieces. The sign in the formula is
positive for 2x2x2, 1x1x2, 1x2x1, 2x1x1, while negative for 2x2x1,
2x1x2, 1x2x2, 1x1x1.

![](media/image1.png){width="2.25in" height="2.0277777777777777in"}

The LS3DF result is extremely accurate when compared to the full-system
direct LDA calculation. For example, using a cubic 8 atom cell in a
diamond Si structure as our smallest 1×1×1 piece to calculate a Si box
passivated with H atoms, we found that the relative energy error is
8.E-6, which is smaller than the error introduced by other sources of
numerical approximations. The absolute energy error is less than 4
meV/atom = 0.1Kcal/mol and thus better than chemical accuracy, and the
total electron charge density has a relative error of 0.03%, which is
essentially the same as the directly calculated results. Finally, the
atomic force error is 6.4E-5 a.u, which is an order of magnitude smaller
than the typical stopping criterion used in *ab initio* atomic
relaxation. Thus, for all practical purposes, the result of the LS3DF
method can be considered the same as the LDA calculation. Also, note
that the final accuracy depends on the size of the small pieces, so that
for the same accuracy, a large system can be divided into similar sized
fragments, at the cost of a larger number of fragments. This ensures the
linear algorithmic scaling of the problem, good parallelization, and
acceptable accuracy.

## CODE AND SCALING

The LS3DF code is based on the planewave DFT PEtot code. The flow chart
of the LS3DF code is shown in Fig. 2 (right side), in comparison to the
original LDA code (left side). The LS3DF code consists of several
components: PEtot_F, which divides the number of processors into
processor groups, and the calculation of the fragment wavefunctions
ψ~F,I~ by each group for a given fragment potential V~F~(r). It also
calculates the fragment charge density ρ~F~(r) from the wavefunction
ψ~F,I~; Gen_dens patches together the fragment charge densities ρ~F~(r)
to generate the total charge density ρ~tot~(r) of the whole system. The
Poisson step generates the LDA total potential V~tot~(r) from the total
charge density ρ~tot~(r). This step solves the Poisson equation for the
whole system using a global FFT. It also uses the Pulay scheme to mix
the resulting LDA potential that is used in the next iteration. Finally,
Gen_VF generates the fragment potential V~F~(r) from the input total
potential V~tot~(r).

![](media/image2.jpeg){width="5.0in" height="2.817361111111111in"}

Fig. 2. Program flow chart for a conventional LDA method and the LS3DF
method

All the codes in Fig. 2 are parallelized. As a demonstration, we have
used this approach to calculate a \~3000 atom Si~2253~H~652~ quantum
dot. The calculated total charge density of this quantum dot is shown in
Fig. 3a. The calculation took approximately 2 hours on 1024 processors
of the NERSC Seaborg computer.

In Fig. 2, the PEtot_F step is the most time consuming part. The scaling
of this part for the 3000 quantum dot (QD) test case on Seaborg computer
is shown in Fig. 3b. As can be seen, this step scales well up to 1024
processors. We believe it should scale effectively to tens of thousands
of processors since each small processor group solves the fragments
independently.

![](media/image3.png){width="2.0in" height="1.4986111111111111in"}
![](media/image4.jpeg){width="1.6180555555555556in"
height="1.7222222222222223in"}

Fig. 3. a) the charge density isosurface (green) plot of a \~3000 atom
Si quantum dot passivated by H atoms. The pink color indicates the
bonds. b) the speedup (blue line) as a function of the number of
processors used in PEtot_F.

![](media/image5.jpeg){width="2.5347222222222223in"
height="1.6666666666666667in"}Note that there are about 4605 fragments
in this test case. Thus, the number of fragments is slightly larger than
the total number of atoms. The number of atoms in each fragment ranges
from 2 to 150 atoms. By using 1024 processors, we have 64 groups
(assuming 16 processors/group). Thus each group will solve on average 72
fragments. We estimate that, in order to keep a proper load balance,
each processor group at least will need to solve two fragments in
average. As a result, the maximum number of processors that we can
efficiently use will be in the order of 10 times the number of atoms in
the system. Thus for a 100,000 atom system, we should be able to
efficiently use a million processors. For our 3,000 Si atom QD, when we
used 1024 processors, the PEtot_F part of each self-consistent iteration
took about 10 minutes. The Poisson solver we used in this calculation
was based on FFTs. For the 3,000 atom Si QD case, the real space
numerical grid employed was 240×240×240 and it took about 1 minute to
solve the Poisson equation using 128 processors. It took about half a
minute each to finish the Gen_dens and Gen_VF program using 128
processors. Thus, in total, for the 1024 processors calculation above,
it took approximately 12 minutes to finish one self-consistent field
(SCF) step. For such a large QD, it typically takes 10 to 20 SCF steps
(the outer loop in Fig. 2) to converge the self-consistent iteration. As
a result, it will take about 2 to 4 hours to finish one self-consistent
calculation (for a fixed atomic position) for this type of QD system
using 1024 processors.

The comparison between the O(N) LS3DF with the conventional O(N^3^) LDA
method is shown in Fig.4 for the total number of floating point counts
for one SCF step. As we can see, the cross over based on this floating
point operation count is at about 500 atoms, which is similar to the
reported cross over for the localized orbital method \[10\].