CECAM/Psi-k Workshop

Local orbitals and linear-scaling ab initio calculations

CECAM, Lyon, 3-7 September 2001

Abstracts for Wednesday, 5 September 2001

We have developed an ab initio nearly O(N) method that combines an accurate optimized-orbital solution of the electronic structure problem with an efficient Green's function technique for evaluating the quantum conductance. The method has been implemented on massively parallel computers and enables studies of very large systems. This talk will describe both the methodological developments and applications to carbon nanotube-metal contacts. The electronic structure problem is solved using a generalization of the Galli-Parrinello approach [1] to real-space grids, where powerful multigrid-based convergence acceleration techniques can be used. Briefly [2], we minimize the total energy functional

=

for a set of nonorthogonal orbitals . Here, is the density matrix in the basis and is defined by . The electronic density is then given by

.

The conjugate orbitals, defined by where is the overlap matrix, satisfy the relation . It has been shown recently [3] that orbitals generated by this procedure decay faster than either orthogonal Wannier functions or the density matrix. If the orbitals are constrained to be zero outside of fixed localization radii, one obtains a variational O(N) expression for the total energy whose accuracy, compared to the full density-functional solution, depends on the size of the localization radii. Our orbitals are variationally optimized using multigrid preconditioning techniques until they accurately describe the ground state of the system. This procedure allows us to use a small number of orbitals per atom, much smaller than in LCAO or gaussian-based calculations, because the orbitals are optimized on the grid according to their environment. In order to ensure fast convergence and accuracy - even for metallic systems - we use both occupied and unoccupied orbitals. The scaling of the most expensive parts of the calculations is still O(N) due to localization, but there is a small O(N³) part, which is dealt with by parallelizing the necessary subspace diagonalization across many processors on a massively parallel supercomputer.

An optimized electron-density function for a carbon nanotube. Note that although the allowed localization region extends over 6 Å, this function is largely confined to one carbon-carbon bond. The plotting plane is along the surface of the nanotube.

The conductance of the full open system (infinite left lead, conductor, infinite right lead) is evaluated via the transmission function as

,

where is Green's function of the conductor (C) and are functions that describe the coupling of the conductor to the left (L) and right (R) leads. In a general non-orthogonal localized orbital scheme, the Green's function of the whole system can be explicitly written as

where and are the self-energy terms due to the semi-infinite leads, and and are the Hamiltonian and overlap matrices for the localized orbitals in the conductor. The expressions for the self-energies are derived using the formalism of principal layers in the framework of the surface Green's function matching theory. The coupling functions can be easily obtained once the self-energy functions are known. A full description of the method is given in Ref. 4.

The conductance calculations are carried out by expanding all quantities in the above basis of optimized, localized orbitals. It is then possible to efficiently evaluate the quantum conductance of a lead-conductor-lead system in O(N) steps, by dividing the system into principal layers that interact only with their nearest neighbors. Due to the small number of orbitals, the size of the matrices that enter the quantum conductance calculation and the computational cost of the whole procedure are thus minimized, while maintaining the high accuracy of ab initio calculations. This allows for ab initio treatment of the infinite leads in full atomistic detail, and for a complete and consistent description of the coupling of the conductor to the leads.

As an important illustrative example, we have investigated carbon nanotube-metal contacts and explained the anomalously large contact resistance observed in nanotube devices as due to the spatial separation of their conductance eigenchannels. The results for various contact geometries will be briefly described, and strategies for improving device performance will also be discussed.

1. G Galli and M Parrinello, Phys. Rev. Lett. 69, 3547 (1992).
2. J-L Fattebert and J Bernholc, Phys. Rev. B 62, 1713 (2000).
3. L He and D Vanderbilt, Phys. Rev. Lett. 86, 5341 (2001).
4. M Buongiorno Nardelli, J-L Fattebert and J Bernholc, Phys. Rev. B, submitted.

Mixing the constituent functions of a band over k-points affords a set of orthonormal Wannier functions. If also a linear transformation between different bands is applied, the result is a set of Nonorthogonal Generalised Wannier functions (NGWFs). Various kinds of Wannier functions have been calculated in a post-processing context in pseudopotential plane-wave calculations [1]. The desirable property of these functions is their spatial localisation. The charge density can be written in terms of the NGWFs and the density kernel as

The are expanded in a plane-wave basis

(1)

(2)

When the NGWFs decay exponentialy we expect that the restriction of the sum of equation (2) only to delta functions centred within a spherical region would be a good approximation. In this case the number of delta functions used for each NGWF will be independent of the size of the system. The same can be achieved with the number of plane waves in (1) with our FFT box method [2].

A new code we are developing performs density functional calculations by optimising both the density kernel and the NGWFs using the above formalism. It is thus closely related to the traditional plane-wave pseudopotential approach with familiar advantages such as an orthonormal basis specified by a kinetic energy cutoff, calculation of energy terms in real or in reciprocal space and natural applicability to any lattice symmetry. The optimisation of the NGWFs is essential in order to obtain plane-wave accuracy.

H s-fireball	Li s-fireball
H NGWF	Li NGWF

Example: NGWFs of LiH. Top: Initial guesses. Bottom: After the completion of a total energy calculation with this method.

1. N Marzari and D Vanderbilt, Phys. Rev. B 56, 12847 (1997).
2. C-K Skylaris, A A Mostofi, P D Haynes, C J Pickard and M C Payne, Comput. Phys. Commun. (in press).

This talk presents the application of the finite-element (FE) method to large-scale ab initio electronic-structure calculations in solids[1,2]. Conventional electronic structure calculations involve the solution of two differential equations - the Schrodinger equation and Poisson's equation. The FE method is used to solve differential equations in a vast range of engineering applications, and as a result has an extensive set of algorithms and codes for a wide range of problems. We have used a Galerkin approach, which casts the FE method into a variational form, to solve the Schrodinger equation and Poisson's equation for the Bloch-periodic boundary conditions that arise in electronic structure calculations on solids. The approach can be understood as an expansion method that uses a strictly local, piecewise-polynomial basis. Because the basis is polynomial in nature, the method is completely general and its convergence can be controlled systematically. Because the basis is strictly local, it produces sparse, structured matrices, enabling the effective use of iterative solution methods. The result is a variational real-space method that requires no Fourier transforms and is well suited for parallelization. The method thus combines the significant advantages of both basis-oriented and grid-based approaches.

We discuss the construction and properties of the required FE bases, and show how they can be used for the solution of the Schrodinger and Poisson equations in solid-state ab initio electronic structure calculations. We present results for the Schrodinger equation demonstrating the optimal and variational convergence of the method in electronic band structure calculations, and results for the Poisson equation demonstrating its optimal convergence and linear scaling with system size. We indicate the steps that are currently under way to develop a fully self-consistent program based on this approach. We conclude with recent applications of the method to large-scale ab initio positron distribution and lifetime calculations[3] for systems of over 4000 atoms.

This work was performed under the auspices of the U.S. Department of Energy by the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.

1. J E Pask, B M Klein, C Y Fong and P A Sterne, Phys. Rev. B 59, 12352 (1999).
2. J E Pask, B M Klein, P A Sterne and C Y Fong, Computer Physics Commun. 135, 1 (2001).
3. P A Sterne, J E Pask and B M Klein, Applied Surface Science 149, 238 (1999).

Wavelets are a new and versatile basis set in electronic structure calculations. Being localized in both real and Fourier space, wavelets combine the advantages of conventional localized basis sets such as Gaussians or finite elements with the advantages of plane waves. Wavelets are therefore an ideal basis set for the localized orbitals that are the basic quantity in O(N) schemes. Due to the data compression properties of wavelets localized orbitals and density matrices can be represented in a very compact way.