Algorithms for Predicting the Physical Properties of Nanocrystals and Large Clusters

JAMES R. CHELIKOWSKY

Center for Computational Materials, Institute for Computational

Engineering and Sciences, University of Texas at Austin, Austin, Texas 78712, USA; Departments of Physics and Chemical Engineering, University of Texas at Austin, Austin, Texas 78712, USA

1.1 Introduction

Quantum confinement offers one the opportunity to alter physical properties of matter without changing the chemical composition. For example, quantum confinement can be used to tune the optical gap across the visible spectrum in CdSe nanocrystals or it can be used to modify the magnetic moment of ferromagnetic clusters.

However, constructing an efficacious algorithm for predicting the role of quantum confinement and its role in determining the properties of nanocrystals is a difficult task owing to the complexity of nanocrystals, which often contain thousands of atoms. Here, I will illustrate new algorithms for nanoscale systems. These algorithms are especially designed for highly parallel platforms. The nature of these algorithms will be discussed and applied to complex systems such as magnetic clusters, and doped semiconductor nanocrystals.

1.2 The Electronic Structure Problem

Many properties of materials can be predicted by a solution of the Kohn–Sham equation:

[MATHEMATICAL EXPRESSION OMITTED], (1.1)

where Vpion is an ionic pseudopotential, VH is the Hartree or Coulomb potential, and Vxc is the exchange-correlation potential. The Hartree and exchange-correlation potentials can be determined from the electronic charge density. The density is given by

[MATHEMATICAL EXPRESSION OMITTED]. (1.2)

The summation is over all occupied states. The Hartree potential is determined by

[MATHEMATICAL EXPRESSION OMITTED]. (1.3)

This term can be interpreted as the electrostatic interaction of a point charge with the charge density of the system.

A physical interpretation of the exchange-correlation potential is more difficult as it has no classical analogue. This potential can be evaluated using a local density approximation. The central tenet of this approximation is that the total exchange-correlation energy may be written as a universal functional of the density:

[MATHEMATICAL EXPRESSION OMITTED], (1.4)

where εxc is the exchange-correlation energy density. Exc and εxc are to be interpreted as depending solely on the charge density. The exchangecorrelation potential, Vxc, is then obtained as [MATHEMATICAL EXPRESSION OMITTED].

It is not difficult to solve the Kohn–Sham equation (eqn (1.1)) for an atom as the atomic charge density can be taken to be spherically symmetric. Thus, the Kohn–Sham problem reduces to solving a one-dimensional problem. The Hartree and exchange-correlation potentials can be iterated to form a self-consistent field. This atomic solution provides the input to construct a pseudo-potential representing the effect of the core electrons and nucleus. The "ion-core" pseudopotential, Vpion, can be transferred to other systems such as molecules and nanocrystals.

The Kohn–Sham equations represent a nonlinear, self-consistent eigenvalue problem. Typically, a solution is obtained by first approximating the Hartree and exchange-correlation potentials using a superposition of atomic charge densities. The Kohn–Sham equation is then solved using these approximate potentials. From the solution, new wavefunctions and charge densities are obtained and used to construct updated Hartree and exchange-correlation potentials. The process is repeated until the "input" and "output" potentials agree and a self-consistent solution is realised. At this point, the total electronic energy can be computed along with a variety of other electronic properties. Once the Kohn–Sham equation is solved, the total electronic energy, ET, of the system can be evaluated from

[MATHEMATICAL EXPRESSION OMITTED]. (1.5)

The structural energy can be obtained by adding the ion-core electrostatic terms. The forces can be obtained by taking the derivative of the energy with respect to position as from the Hellmann–Feynman theorem.

1.3 Algorithms for Solving the Kohn–Sham Equation

The Kohn–Sham equation as cast in eqn (1.1) can be solved using a variety of techniques. Often, the wavefunctions can be expanded in a basis such as plane waves or gaussians and the resulting secular equations can be solved using standard diagonalisation packages such as those found in VASP.

Here, we focus on a different approach. We solve the Kohn–Sham equation without resort to an explicit basis. We solve for the wavefunctions on a grid with a fixed domain, which encompasses the physical system of interest. The grid need not be uniform, but it greatly simplifies the problem if it is. The wavefunctions outside of the domain are required to vanish for confined sys- tems or assume periodic boundary conditions for systems with translational symmetry. In contrast to methods employing an explicit basis, such boundary conditions are easily incorporated. In particular, real-space methods do not require the use of supercells for localised systems. As such, charged systems can easily be examined without considering any electrostatic divergences. The Hartree potential is also solved on the grid from eqn (1.3) using conjugate-gradient minimisation. For confined systems, the boundary condition for this Poisson equation can be obtained by first performing a multipole expansion of the electronic charge density, the Hartree potential on the boundary of the domain can then be evaluated using the multipoles. A few lowest-order multipoles are found to be sufficient to obtain an accurate Hartree potential.

Within a "real-space" approach, one can solve the eigenvalue problem using a finite element or finite difference approach. We use a higher-order finite difference approach owing to its simplicity in implementation. The Laplacian operator can be expressed using

[MATHEMATICAL EXPRESSION OMITTED], (1.6)

where h is the grid spacing, N is the number of nearest grid points, and Ci are the coefficients for evaluating the required derivatives. The error scales as O(h2N+2).

Once the secular equation is created, the eigenvalue problem can be solved using iterative methods. Typically, a method such as a preconditioned Davidson method can be used. This is a robust and efficient method, which never requires one to store the Hamiltonian matrix. In this chapter, we outline a method that avoids solving large-eigenvalue problems explicitly. The method utilises a damped Chebyshev polynomial filtered subspace iteration. In this approach, only the initial iteration requires solving an eigenvalue problem, which can be handled by means of any available efficient eigensolver. This step is used to provide a good initial subspace (or good initial approximation to the wavefunctions). Because the subspace dimension is slightly larger than the number of wanted eigenvalues, the method does not require as much memory as standard restarted eigensolvers such as ARPACK and TRLan (Thick – Restart, Lanczos). Moreover, the cost of orthogonalisation is much reduced as the filtering approach only requires a subspace with dimension slightly larger than the number of occupied states and orthogonalisation is performed only once per SCF iteration. In contrast, standard eigensolvers using restart usually require a subspace twice as large and the orthogonalisation and other costs related to updating the eigenvectors are much higher.

The main idea of the proposed method is to start with a good initial eigenbasis, {ψn}, corresponding to occupied states of the initial Hamiltonian, and then to improve adaptively the subspace by polynomial filtering. That is, at a given self-consistent step, a polynomial filter, Pm(H), of order m is constructed for the current Hamiltonian H. As the eigenbasis is updated, the polynomial will be different at each SCF step since H will change. The goal of the filter is to make the sub space spanned by [MATHEMATICAL EXPRESSION OMITTED] approximate the eigensubspace corre sponding to the occupied states of H. There is no need to make the new subspace, {[??]n}, approximate the wanted eigensubspace of H to high accuracy at inter mediate steps. Instead, the filtering is designed so that the new subspace obtained at each self-consistent iteration step will progressively approximate the wanted eigenspace of the final Hamiltonian when self-consistency is reached.

This can be efficiently achieved by exploiting the Chebyshev polynomials, Cm, for the polynomials Pm. Specifically, we wish to exploit the fast growth property outside of the [–1, 1] interval. To obtain a good filter for a given SCF step, one needs to provide a lower bound and an upper bound of the unoccupied spectrum of the current Hamiltonian H. The lower bound can be readily obtained from the Ritz values computed from the previous step, and the upper bound can be inexpensively obtained by a very small number of (e.g., 4 or 5) Lanczos steps. Hence, the main cost of the filtering at each iteration is in computing the polynomial operation.

One can use an affine mapping, l(t), to map interval [a, b] to the interval [–1, 1]:

[MATHEMATICAL EXPRESSION OMITTED]. (1.7)

The interval is chosen to encompass the energy interval that needs to be filtered out. The filtering operation can then be expressed as

[MATHEMATICAL EXPRESSION OMITTED]. (1.8)

This computation is accomplished by exploiting the convenient three-term recurrence property of Chebyshev polynomials:

[MATHEMATICAL EXPRESSION OMITTED]. (1.9)

An example of a damped Chebyshev polynomial as defined by eqns (1.7) and (1.9) is given in Figure 1.1, where we have taken the lower bound as a = 0.2 and the upper bound as b = 2. In this example, the filtering would enhance the eigenvalue components in the shaded region.

The filtering procedure for the self-consistent cycle is illustrated in Figure 1.2. Unlike traditional methods, the cycle only requires one explicit diagonalisation step. Instead of repeating this step again within the self-consistent loop, a filtering operation is used to create a new basis in which the desired eigensubspace is enhanced. After the new basis, {ψn}, is formed, the basis is orthogonalised. The orthogonalisation step scales as the cube of the number of occupied states and as such this method is not an "order-n" method, i.e. it does not scale linearly with the number of occupied states. However, the prefactor is sufficiently small that the method is much faster than previous implementations of real-space methods. As for the traditional methods, the self-consistent cycle is repeated until the "input" and "output" density is unchanged.

The algorithm can be parallelised by partitioning the domain that encloses the physical system into continuous subdomains of grid points as implemented in the PARSEC code. The number of subdomains is the same as the number of processors used in the calculation. A mapping can then be defined that assigns a subdomain to a processor. For optimal load balancing, different subdomains contain roughly the same number of grid points. The elements in the eigenvectors and potential arrays are indexed by grid points, hence a continuous block within the arrays can be distributed to different processors according to the defined mapping. As a result, the data required by the algorithm is completely distributed, and the Hamiltonian matrix is not stored explicitly. The Chebyshev filtering method requires matrix-vector operation between the Hamiltonian matrix and the eigenvectors. The operations of the local part of the ionic pseudopotential, Hartree and exchange-correlation potentials on the eigenvectors are local and therefore do not require communication. On the other hand, the Laplacian operator in eqn (1.6) and the nonlocal part of the ionic pseudopotential require communication with the processors that correspond to neighbouring blocks of grid points. The communication is done using the standard message passing interface (MPI) library.

In Table 1.1, we compare the timings using the Chebyshev filtering method along with explicit diagonalisation solvers using the TRLan and the ARPACK. These timing are for a modest sized nanocrystal: Si525H276. The Hamiltonian size is 292 584 × 292 584 and 1194 eigenvalues were determined. The numerical runs were performed on the SGI Altix 3700 cluster at the Minnesota Supercomputing Institute. The CPU type is a 1.3-GHz Intel Madison processor. Although the number of matrix-vector products and SCF iterations is similar, the total time with filtering is over an order of magnitude faster compared to ARPACK, and a factor of better than four when compared to TRLan. Such improved timings are not limited to this particular example. We have also explored GaAs nanocrystal and large Fe clusters, including spin-dependent density functionals. In Figure 1.3, we present the scaling behaviour of this algorithm on Lonestar, which is a TACC Dell Linux cluster containing 5840 computing cores. The scaling is studied using an iron cluster containing 43 iron atoms. The Hamiltonian matrix has a size of 110 245 and 258 eigenstates were calculated. We can see that the scaling starts linearly with the number of processors, but flattens off at around 128 processors for this particular system. The number of processors at which the saturation of the scaling occurs depends on the system size, a larger system size will need more processors to saturate the scaling.

1.4 Applications

1.4.1 The Electronic Properties of Si Nanocrystals

We illustrate the Chebyshev filtering method for a prototypical system: hydrogenated silicon nanocrystals. These nanocrystals, or quantum dots, are small fragments of the bulk in which the surface has been passivated by hydrogen atoms. In the case of silicon, the passivation is accomplished experimentally by capping the surface dangling bonds with hydrogen atoms. These systems exhibit interesting changes as one approaches the nanoregime. For example, nanocrystals of silicon are expected to be optically active, whereas bulk crystals of silicon are not.

The largest dot we examined contained over ten thousand atoms: Si9041H1860. This dot is approximately 7 nm in diameter. A ball and stick model for a much smaller quantum dot, approximately, 3 nm in diameter, is illustrated in Figure 1.4.