Chris-Kriton Skylaris

Introducing ONETEP: Linear-scaling density functional simulations on parallel computers

We present ONETEP (order-N electronic total energy package), a density functional program for parallel computers whose computational cost scales linearly with the number of atoms and the number of processors. ONETEP is based on our reformulation of the plane wave pseudopotential method which exploits the electronic localization that is inherent in systems with a nonvanishing band gap. We summarize the theoretical developments that enable the direct optimization of strictly localized quantities expressed in terms of a delocalized plane wave basis. These same localized quantities lead us to a physical way of dividing the computational effort among many processors to allow calculations to be performed efficiently on parallel supercomputers. We show with examples that ONETEP achieves excellent speedups with increasing numbers of processors and confirm that the time taken by ONETEP as a function of increasing number of atoms for a given number of processors is indeed linear. What distinguishes our approach is that the localization is achieved in a controlled and mathematically consistent manner so that ONETEP obtains the same accuracy as conventional cubic-scaling plane wave approaches and offers fast and stable convergence. We expect that calculations with ONETEP have the potential to provide quantitative theoretical predictions for problems involving thousands of atoms such as those often encountered in nanoscience and biophysics.

Nonorthogonal generalized Wannier function pseudopotential plane-wave method

We present a reformulation of the plane-wave pseudopotential method for insulators. This new approach allows us to perform density-functional calculations by solving directly for “nonorthogonal generalized Wannier functions” rather than extended Bloch states. We outline the theory on which our method is based and present test calculations on a variety of systems. Comparison of our results with a standard plane-wave code shows that they are equivalent. Apart from the usual advantages of the plane-wave approach such as the applicability to any lattice symmetry and the high accuracy, our method also benefits from the localization properties of our functions in real space. The localization of all our functions greatly facilitates the future extension of our method to linear-scaling schemes or calculations of the electric polarization of crystalline insulators.

Preconditioned iterative minimization for linear-scaling electronic structure calculations

Linear-scaling electronic structure methods are essential for calculations on large systems. Some of these approaches use a systematic basis set, the completeness of which may be tuned with an adjustable parameter similar to the energy cut-off of plane-wave techniques. The search for the electronic ground state in such methods suffers from an ill-conditioning which is related to the kinetic contribution to the total energy and which results in unacceptably slow convergence. We present a general preconditioning scheme to overcome this ill-conditioning and implement it within our own first-principles linear-scaling density functional theory method. The scheme may be applied in either real space or reciprocal space with equal success. The rate of convergence is improved by an order of magnitude and is found to be almost independent of the size of the basis.

Accurate ionic forces and geometry optimization in linear-scaling density-functional theory with local orbitals

Linear scaling methods for density-functional theory (DFT) simulations are formulated in terms of localized orbitals in real space, rather than the delocalized eigenstates of conventional approaches. In local-orbital methods, relative to conventional DFT, desirable properties can be lost to some extent, such as the translational invariance of the total energy of a system with respect to small displacements and the smoothness of the potential-energy surface. This has repercussions for calculating accurate ionic forces and geometries. In this work we present results from onetep, our linear scaling method based on localized orbitals in real space. The use of psinc functions for the underlying basis set and on-the-fly optimization of the localized orbitals results in smooth potential-energy surfaces that are consistent with ionic forces calculated using the Hellmann-Feynman theorem. This enables accurate geometry optimization to be performed. Results for surface reconstructions in silicon are presented, along with three example systems demonstrating the performance of a quasi-Newton geometry optimization algorithm: an organic zwitterion, a point defect in an ionic crystal, and a semiconductor nanostructure.

Total-energy calculations on a real space grid with localized functions and a plane-wave basis

We present a novel real space formalism for ab initio electronic structure calculations. We use localized non-orthogonal functions that are expressed in terms of a basis set that is equivalent to a plane-wave basis. As a result, advantages of the plane-wave approach also apply to our method: its applicability to any lattice symmetry, and systematic basis set improvement via the kinetic energy cut-off parameter. The localization of our functions enables the use of fast Fourier transforms over small regions of the simulation cell to calculate the total energy with efficiency and accuracy. With just one further variational approximation, namely the truncation of the density matrix, the calculation may be performed with a cost that scales linearly with system size for insulating systems.

Minimal parameter implicit solvent model for ab initio electronic-structure calculations

We present an implicit solvent model for ab initio electronic-structure calculations which is fully self-consistent and is based on direct solution of the nonhomogeneous Poisson equation. The solute cavity is naturally defined in terms of an isosurface of the electronic density according to the formula of Fattebert and Gygi (J. Comput. Chem., 23 (2002) 662). While this model depends on only two parameters, we demonstrate that by using appropriate boundary conditions and dispersion-repulsion contributions, solvation energies obtained for an extensive test set including neutral and charged molecules show dramatic improvement compared to existing models. Our approach is implemented in, but not restricted to, a linear-scaling density functional theory (DFT) framework, opening the path for self-consistent implicit solvent DFT calculations on systems of unprecedented size, which we demonstrate with calculations on a 2615-atom protein-ligand complex.

Electrostatic interactions in finite systems treated with periodic boundary conditions: application to linear-scaling density functional theory

We present a comparison of methods for treating the electrostatic interactions of finite, isolated systems within periodic boundary conditions (PBCs), within density functional theory (DFT), with particular emphasis on linear-scaling (LS) DFT. Often, PBCs are not physically realistic but are an unavoidable consequence of the choice of basis set and the efficacy of using Fourier transforms to compute the Hartree potential. In such cases the effects of PBCs on the calculations need to be avoided, so that the results obtained represent the open rather than the periodic boundary. The very large systems encountered in LS-DFT make the demands of the supercell approximation for isolated systems more difficult to manage, and we show cases where the open boundary (infinite cell) result cannot be obtained from extrapolation of calculations from periodic cells of increasing size. We discuss, implement, and test three very different approaches for overcoming or circumventing the effects of PBCs: truncation of the Coulomb interaction combined with padding of the simulation cell, approaches based on the minimum image convention, and the explicit use of open boundary conditions (OBCs). We have implemented these approaches in the ONETEP LS-DFT program and applied them to a range of systems, including a polar nanorod and a protein. We compare their accuracy, complexity, and rate of convergence with simulation cell size. We demonstrate that corrective approaches within PBCs can achieve the OBC result more efficiently and accurately than pure OBC approaches.

Electrostatic embedding in large-scale first principles quantum mechanical calculations on biomolecules.

Biomolecular simulations with atomistic detail are often required to describe interactions with chemical accuracy for applications such as the calculation of free energies of binding or chemical reactions in enzymes. Force fields are typically used for this task but these rely on extensive parameterisation which in cases can lead to limited accuracy and transferability, for example for ligands with unusual functional groups. These limitations can be overcome with first principles calculations with methods such as density functional theory (DFT) but at a much higher computational cost. The use of electrostatic embedding can significantly reduce this cost by representing a portion of the simulated system in terms of highly localised charge distributions. These classical charge distributions are electrostatically coupled with the quantum system and represent the effect of the environment in which the quantum system is embedded. In this paper we describe and evaluate such an embedding scheme in which the polarisation of the electronic density by the embedding charges occurs self-consistently during the calculation of the density. We have implemented this scheme in a linear-scaling DFT program as our aim is to treat with DFT entire biomolecules (such as proteins) and large portions of the solvent. We test this approach in the calculation of interaction energies of ligands with biomolecules and solvent and investigate under what conditions these can be obtained with the same level of accuracy as when the entire system is described by DFT, for a variety of neutral and charged species.

On the resolution of identity Coulomb energy approximation in density functional theory

The Resolution of the Identity approximation for the Coulomb (RI-J) energy in Density Functional Theory improves the computational efficiency of large-scale calculations but requires the use of a second, or “auxiliary” basis set. We examine the performance of some of the existing auxiliary basis sets with a variety of basis sets and molecules. We determine the accuracy of the RI-J approximation for these basis sets and suggest criteria for the selection of combinations of basis set and auxiliary basis set.

Accurate kinetic energy evaluation in electronic structure calculations with localized functions on real space grids

We present a method for calculating the kinetic energy of localized functions represented on a regular real space grid. This method uses fast Fourier transforms applied to restricted regions commensurate with the simulation cell and is applicable to grids of any symmetry. In the limit of large systems it scales linearly with system size. Comparison with the finite difference approach shows that our method offers significant improvements in accuracy without loss of efficiency.