The ONETEP linear-scaling density functional theory program

A common starting point for performing density functional theory (DFT) calculations with linear-scaling computational cost is to use the one-particle density matrix

which provides a complete description of the fictitious Kohn–Sham (KS) system, where {ψ_i} are the KS orbitals and {f_i} are their occupancies. KS DFT calculation methods working directly with the density matrix have been developed such as the approach by Li, Nunes, and Vanderbilt.¹ The density matrix is a quantity that decays exponentially with the distance |r − r′| (for gapped systems at zero temperature and metals at a finite temperature), as formulated by Walter Kohn in his “nearsightedness of electronic matter” principle.^2,3 To take practical advantage of this principle and develop methods with reduced or linear-scaling computational cost, we need to truncate the exponentially decaying tail of the density matrix. As the system size (number of atoms) increases, we will eventually reach the point where the remaining amount of information increases linearly with the size of the system. This is not straightforward to implement if we use the KS orbitals directly. A more practical approach is to work with a set of localized orbitals {ϕ_α} that, in general (but not necessarily), are non-orthogonal as non-orthogonality can result in better localization.⁴ These localized orbitals can be considered to be connected to the KS orbitals via a linear transformation M as follows:

where we have used Greek indices for non-orthogonal components and Latin indices for orthogonal components. We also use the Einstein implicit summation convention for repeated greek indices. We will follow these conventions throughout this paper. It is easy to show that the density matrix can be expressed in terms of the local orbitals in the following form:^5,6

where the matrix K, which we call the density kernel, is the representation of the density matrix in the basis of the localized orbitals as follows:

Linear-scaling schemes can be developed by following this framework. Calculations with such schemes need to obey the idempotency property for the density kernel

which is equivalent to the necessary requirements of KS orbital orthonormality and integer occupancies of 1 or 0 (for calculations in materials with a gap).

One class of linear-scaling approach typically uses a fixed basis set of localized functions, such as Gaussian atomic orbitals or numerical atomic orbitals. However, it is difficult to systematically improve the quality of such basis sets toward convergence while simultaneously maintaining linear-scaling behavior, as this requires the inclusion of a large number of extended atomic basis functions, which can make it difficult to impose the necessary localization constraints on the density matrix.

To overcome this limitation, methods have emerged, which aim to achieve linear-scaling cost with near-complete basis set accuracy or “plane wave accuracy.” The main idea in these methods is to optimize both the non-orthogonal localized orbitals {ϕ_α} and the density kernel K, both subject to localization constraints, so as to achieve the variational freedom of a large basis set while retaining the localization and small matrix sizes of a minimal AO basis. onetep⁷ belongs to this class of methods; other codes in this category include conquest⁸ and bigdft.⁹ In onetep, the localized orbitals {ϕ_α} are called non-orthogonal generalized Wannier functions (NGWFs) and are optimized self-consistently during the calculation.¹⁰ 

In onetep, we aim to achieve controllable accuracy equivalent to that of plane-wave pseudopotential calculations, so the basis set we use to expand the local orbitals is a set of periodic sinc (psinc) basis functions {D_i},^10,11 which are unitary transformations of plane waves. The quality of the psinc set is controlled by a single parameter, the psinc kinetic energy cut-off, in analogy with plane wave codes. Thus, the psinc basis set provides a uniform description of space, similar to plane waves, and the kinetic energy cutoff controls the resolution of this description, providing a finer resolution with an increase in the cut-off value. Typical kinetic energy cutoff values in routine calculations are in the range of 500–1000 eV, although higher or lower values are possible and have been used in specialized cases.

Another advantage of using a psinc basis set is that many operations on the NGWFs that are needed for the construction of the Hamiltonian in the NGWF representation can be performed efficiently using Fast Fourier Transforms (FFTs), as in plane wave codes. This involves transforming the NGWFs to reciprocal space, applying the appropriate operator in diagonal form (e.g., the kinetic energy operator), and transforming back into real space. However, unlike in plane wave codes, these operations need to be performed with computational effort that is both relatively small and independent of the size of the system (i.e., independent of the number of atoms). To achieve this, we perform the FFTs on NGWFs in a miniature simulation cell, which we call the “FFT box.”¹² In order to preserve important properties of the integrals [such as hermiticity and uniformity of the quantum mechanical (QM) operators] computed with the FFT box, particular care must be taken in how this box is defined and used, as we have described in early publications.¹³ These publications also present the use of the FFT box for the efficient construction of the electronic density and the NGWF gradient.¹⁰ 

The ground state KS energy can now be considered as the minimum of a functional of the NGWFs and the density kernel

where the local potential $V^loc=V^ps,loc+V^H$ is the sum of the Hartree potential and the local pseudopotential, $V^nl$ is the non-local part of the pseudopotential in the Kleinman and Bylander representation, and $n=nr=ρr,r$ is the electronic density.

As the NGWFs are optimized to minimize the ground state energy, it turns out that they can provide a good description of the valence bands, but a very poor description of the conduction bands. For this reason, in order to provide a set of functions suitable for describing the conduction bands (e.g., for excited state calculations), we can optimize a second set of NGWFs {χ_α}, after the completion of the ground state calculation, which provide an accurate description of a set of low-lying conduction bands (see Sec. IV E 1).

onetep was designed from its inception as a parallel code,¹⁴ and the extent and efficiency of its parallelism have continued to improve during its development.^15,16 Originally created within the paradigm of the Message Passing Interface (MPI), the code is now able to extensively harness “hybrid” parallelism enabled by OpenMP multithreading. MPI divides the available computational resources into processes, each of which can then be subdivided by OpenMP into threads, each occupying a single physical core. Each of the principal data structures and computational tasks is distributed in a manner designed to best exploit the locality on which the linear-scaling formalism relies. Data and tasks associated with atoms are distributed according to an algorithm designed to balance the competing demands of locality, aimed at ensuring nearby atoms are on the same or nearby processes to minimize communication, and load-balance, aimed at ensuring all processes have nearly equal numbers of NGWFs. This algorithm first sorts all the atoms based on a Peano space-filling curve¹⁷ and then assigns them to processes based on a modified greedy algorithm,¹⁸ looping over sets of atoms with decreasing numbers of NGWFs so that the largest are assigned first.

Data associated with each atom, particularly the coefficients of each atom-centered NGWF, are distributed and held locally on the process of that atom. The functions of each atom are then ideal for a second level of parallelism, implemented using OpenMP threading, whereby each thread treats a single function at a time, for example, when applying an operator to the set of functions. As demonstrated in Ref. 16, this leads to efficient scaling with respect to both MPI processes and threads, enabling the code to take advantage of very large-scale HPC (high-performance computing) architectures, with efficient runs on as many as 16 384 parallel cores. Subsequent work has extended the threaded parallelism to the implicit solvent model and exact exchange functionalities, detailed in Secs. IV D 1 and IV B 1. However, more optimization is always possible, and a current limitation is that MPI communications within threaded regions are handled with CRITICAL regions that are only entered by one thread at a time, which causes whole-cell grid operations to scale relatively poorly with thread count beyond around 8–10 threads.

Figures 1 and 2 demonstrate the linear-scaling capabilities of onetep in two very different systems. Figure 1 presents results for amyloid fibril segments of varying sizes using several different functionals, demonstrating linear-scaling behavior at least up to 14 000 atoms.¹⁹ Figure 2 presents results for slabs of anatase-structure TiO₂ of varying sizes and shows linear-scaling behavior at least up to 3500 atoms.²⁰ 

In this section, we describe the various capabilities of onetep. Each feature has its own subsection, which includes a brief introduction to the functionality as well as relevant citations to other publications for those interested in either the theory or the application of the functionality. The capabilities are grouped thematically—the themes are “Improvements to core algorithms and functionality,” “Beyond-DFT methods,” “Dynamics and structure optimization,” “Including environmental effects,” “Excitations and spectroscopy,” and “Partitioning the density/energy.”

Blöchl’s projector augmented wave (PAW) approach enables all-electron (AE) calculations with the efficiency, systematic accuracy, and transferability of plane-wave pseudopotential calculations.²¹ PAW works by relating the soft pseudo-(PS-)wavefunctions, easily represented by plane-waves or on a grid, to AE wavefunctions, varying rapidly near the nuclei, by means of a linear transformation augmenting the soft part of the wavefunction with partial waves near each atom. PAW increases the accuracy feasible at a given computational cost, enabling efficient treatment of, for example, transition metal oxides, lanthanides, and actinides, which would be highly challenging for norm-conserving pseudopotentials (NCPPs). PAW also allows calculation of quantities requiring knowledge of wavefunctions, densities, and electric fields in the regions near to the atomic nuclei, which are not correctly represented by norm-conserving methods. Hine²⁰ introduced a novel adaptation of the PAW formalism suited to linear-scaling DFT within onetep. The nomenclature here follows that of Ref. 22, which also contains an introduction to the PAW method as applied to traditional DFT.

The fundamental PAW transformation relates AE orbitals |ψ_n⟩ to PS orbitals $|ψ̃n$⁠,

$|ψ̃n$ are the projectors, and |φ_ν⟩ and $|φ̃ν$ are the AE and PS partial waves (see Refs. 22 and 20), respectively.

Within the linear-scaling DFT approach, we do not have direct access to the eigenstates, so we can instead use an equivalent transformation to relate the AE density matrix to the PS density matrix,

The NGWFs in onetep, $ϕαr$⁠, which are constructed out of psinc functions and thus equivalent to the plane waves of standard PAW, are an ideal means to represent this “soft” part of the density matrix (DM), so one can define

Notably, this allows one to retain the same basic forms of almost all the optimization algorithms for the kernel and NGWFs, having first redefined the overlap matrix between NGWFs to account for the PAW overlap operator,

Further complications associated with the representation of the compensation density, non-local potential energy terms, NGWF gradient, forces, and preconditioning must also be dealt with: these are discussed in detail in Ref. 20. The result is a method that has demonstrably superior convergence properties compared to NCPPs (see Fig. 2 of Ref. 20) and thus can achieve greater accuracy for a given computational cost. At fixed parameters such as NGWF radius and cutoff energy, the overhead of a PAW calculation is relatively low, at around 10%–20% compared to NCPPs. Linear-scaling computational cost remains just as achievable as for NCPPs, since the majority of the terms in PAW appear naturally on a per-atom basis. Reference 20 also shows excellent parallel scaling, demonstrating efficient simulation of 1824 atoms on 1920 cores.

Ensemble density functional theory (EDFT)^23–25 calculations are available in onetep. Within this formalism, the KS states have fractional occupancies, which can be determined in several different ways, including the Methfessel–Paxton,²⁶ Gaussian smearing,^27,28 and Fermi–Dirac²⁴ schemes. In onetep, the occupancies are determined by the Fermi–Dirac distribution, regulated by a constant electronic temperature. The macroscopic state function is the Helmholtz free energy, which obeys the variational principle of quantum mechanics. onetep directly minimizes the Helmholtz free energy functional to self-consistently find the KS states and their fractional occupancies.²⁹ 

The onetep approach requires a Hamiltonian diagonalization step in each SCF (self-consistent field) iteration, which makes the overall cost scale with the cube of the system size. However, the cost of the Hamiltonian diagonalization is kept to a minimum as a consequence of using a minimal set of NGWFs. In this way, the dimension of the Hamiltonian matrix is reduced to the minimum possible, and the prefactor due to diagonalization is reduced significantly. The remaining steps of the method can be performed at a linear-scaling cost due to the sparsity of the matrices required in the NGWF representation.

Additionally, parallel diagonalization techniques can be used to scale the diagonalization of the Hamiltonian to many computational cores and to distribute the associated memory requirements. As a result, simulations on metallic systems with thousands of atoms can be achieved using the onetep EDFT implementation.

In free-spin EDFT calculations, the charge and net spin of a system do not necessarily have to be integers, and the net spin can additionally be relaxed during a calculation. This makes the exploration of systems with unknown magnetic properties much easier.

As in linear-scaling methods for insulators, the idea behind the Fermi operator expansion (FOE)^30,31 is to be able to perform an SCF calculation without using a cubic-scaling diagonalization step. FOE, however, also works for metallic systems where we want to apply fractional occupancies to the KS states.

In FOE, the density matrix, with finite-temperature electronic occupancies, is constructed as a polynomial expansion of the matrix form of the Fermi–Dirac occupancy formula,

where I is the identity matrix, H is the Hamiltonian matrix, S is the overlap matrix, μ is the chemical potential, and β = 1/k_BT. Equation (11) cannot be applied directly as the condition number of the matrix to be inverted will be too large in general. Many options for expansions have been proposed over the last three decades.^32–38 In onetep, we have implemented the AQuA-FOE method described in detail in Ref. 39. AQuA-FOE takes advantage of the trigonometric mapping between the Fermi–Dirac occupancy formula and the hyperbolic tangent function. In this way, we can construct the density kernel at an electronic temperature many times that of the target temperature. The target temperature density kernel can then be recovered by repeatedly halving the temperature of the density kernel with the matrix analog of the hyperbolic tangent double angle formula. A fixed expansion length is always employed at the high temperature to generate the density matrix, requiring fewer terms to reach a given accuracy than at lower temperature.

To build the high temperature density kernel, we use a Chebyshev expansion of Eq. (11) using the fast re-summing algorithm of Liang et al.³³ As an alternative option, this algorithm can be used as the sole FOE method in onetep without the annealing and quenching steps.

The AQuA-FOE method in onetep can also find the chemical potential that conserves electron number by using a safeguarded Newton’s algorithm and the analytic derivative of the number of electrons with respect to the chemical potential. The chemical potential of the density kernel is then updated using the addition theorem of the hyperbolic tangent function so that the expansion need not be recalculated from scratch.

The AQuA-FOE method is presently invoked within the EDFT method (Sec. IV A 2) in onetep in order to produce finite temperature density kernels for each trial Hamiltonian matrix. In the future, this method could be used within density-mixing type or direct minimization algorithms.

onetep exploits the nearsightedness of the electronic structure (Sec. I) by using NGWFs that are strictly localized within spherical regions and which give rise to sparse matrix representations that may be truncated in real space in a systematically controllable manner and enable scalability of the method to very large systems. There are cases, however, in which it is advantageous to use a hybrid approach, in which the representation of electronic structure is localized along certain directions (Wannier-like) and delocalized along others (Bloch-like). Use cases of such a hybrid approach include layered systems and surface slabs, in which localization is applied only in the direction normal to the layer or slab, and nanowires, in which localization is applied in the two directions normal to the nanowire axis.

The concept of such so-called “hybrid Wannier functions” was first introduced by Sgiarovello et al.⁴⁰ This idea, within the formalism of maximally localized Wannier functions,⁴¹ was subsequently used, for example, to study the dielectric properties of layered materials^42,43 and topological insulators.^44,45 In these applications, the hybrid Wannier functions were obtained by post-processing the Bloch eigenstates from a traditional DFT calculation.^46,47 In onetep, we have introduced an approach in which hybrid NGWFs are optimized directly in situ, avoiding the system-size bottleneck associated with the O(N³) scaling of the traditional DFT calculation. We have also included the flexibility for the user to choose the lattice vector directions that are represented by a localized description and those that are extended, providing the full range of options in going from fully localized NGWFs to hybrid NGWFs that are localized in either two directions or one direction to a fully extended Bloch-like representation (see Fig. 3). Just as in the case of fully localized NGWFs, the hybrid NGWFs are constrained to be strictly zero outside of their localization regions.

Figure 4 shows the total energy per atom, calculated with hybrid NGWFs that are localized along the [001] direction (“2D-HNGWFs”) with different localization lengths r_cut, for a 2 × 2 × 10 slab supercell of an eight-atom cubic diamond carbon conventional unit cell. The slab is terminated with hydrogen atoms on its {001} surfaces, and periodic replicas in the [001] direction are separated by 15 Å vacuum regions. We use a kinetic energy cutoff of 1200 eV for the psinc basis, norm-conserving pseudopotentials for the C and H atoms, the local-density approximation (LDA) for exchange and correlation, and the Brillouin zone is sampled at the Γ-point only. The carbon diamond lattice parameter is set to its equilibrium bulk value within LDA-DFT, which we calculate to be 3.539 Å. The energies shown in Fig. 4 are referenced to an equivalent calculation with the traditional plane-wave code castep.⁴⁹ As can be seen, with hybrid NGWFs, it is possible to attain equivalence with a traditional plane-wave approach using a localization length of around 5 Å.

While the relaxation of the localization constraints results in a less favorable scaling of computational effort with system-size than using fully localized orbitals, the computational cost still only scales linearly with the size of the system along the directions in which the localization constraints are still applied (e.g., the layer/slab normal in a layered/slab system). Furthermore, in our experience, there are additional advantages that offset the nominal increase in computational cost per self-consistent iteration: the use of hybrid NGWFs has a tendency to improve the conditioning of the energy minimization, resulting in faster convergence; and the accuracy of the ionic forces tends to be better with a reduction in the “egg-box” effect^50–54 that is usually present in real-space local-orbital methods. For further details on these issues, refer to Ref. 48.

dl_mg⁵⁵ is a library that solves variants of the Poisson equation with the general form

in 3D over a cuboid domain with periodic, Dirichlet or mixed boundary conditions. A regular discretization grid is used over which the derivatives are approximated with finite differences of adjustable order. The library is developed in close coordination with onetep and is distributed as part of the package.

In Eq. (12), $εr$ is the relative permittivity, $vr$ is the electrostatic potential (ESP), and $ntotr$ is the fixed charge density. These terms feature in the generalized Poisson equation (GPE)

The second term on the right-hand side of Eq. (12) appears in the Poisson–Boltzmann equation (P–BE), representing the charge density of mobile ions (e.g., electrolyte). In this term, z_i and $c¯i$ are the electric charge and the average bulk concentration (N_i/V) for the mobile ion of type i, respectively, β = 1/k_BT is the inverse temperature, and $0≤λr≤1$ is the ion accessibility function that accounts for the short range repulsion effect between mobile ions and fixed ionic cores. α and γ are constants that depend on the used unit systems.

In addition to solving the non-linear Poisson–Boltzmann equation [NLP–BE, Eq. (12)], dl_mg employs specialized algorithms to solve both the GPE [Eq. (13)] and linearized Poisson–Boltzmann equation (LP–BE),

In brief, dl_mg finds the solution for the GPE or P–BE in two stages:

1. A solution is found for the second order finite difference discretization using a multigrid algorithm for the GPE/LP–BE or a combination of the inexact Newton method and multigrid for the NLP–BE.

2. The solution corresponding to the chosen finite difference order (if larger than two) is found with a defect correction procedure that starts from the second order solution. The largest finite difference order implemented is 12.

See Ref. 55 for a detailed description of the theory and methods employed in dl_mg.

dl_mg underpins several functionalities of onetep, including implicit solvent (Sec. IV D 1) and open boundary conditions (Sec. IV D 4), which utilize the flexible Poisson-solving capabilities of the library.

On the well-known Jacob’s ladder of density functional approximations^56,57 (Fig. 5), orbital-dependent functionals occupy the higher rungs above the local density approximation (LDA) and the generalized gradient approximation (GGA). Dependence on the form of the KS orbitals is added on top of the usual dependence on density n in LDA and density gradient ∇n in GGA. This increases the flexibility available in the functional forms, enabling the satisfaction of additional physical constraints and/or better fitting of empirical data. The upshot of this is that improved accuracy is possible, though not guaranteed.

Two families of orbital-dependent functionals are currently supported in onetep: meta-GGAs and hybrids, which depend on the KS orbitals via the kinetic energy density, τ, and the Hartree–Fock exchange energy, $ExHF$⁠, respectively. The primary challenge of supporting both families in onetep is maintaining the onetep formalism’s $O(N)$ scaling when evaluating quantities with dependence on the KS orbitals.

The orbital dependence of the meta-GGAs is via τ, i.e.,

with

where the summation is over the occupied KS orbitals. In onetep, τ is evaluated in terms of the NGWFs and density kernel,

This quantity has the same structure as the charge density n, but with ∇ϕ_α in place of ϕ_α, and can be evaluated using the same procedure,¹⁴ subject to the subtle approximation that extra overlapping pairs of NGWFs αβ that would result from the delocalization of {ϕ_α} upon application of the gradient operator are neglected.¹⁹ 

Once τ is constructed, evaluation of $ExcMGGA$ is simple, since the semi-local exchange–correlation energy density ϵ_xc(r) depends only on n, ∇n, and τ at point r. The difficulty arises in evaluating the exchange–correlation potential,

and the quantities that depend on this, such as the gradient of E[{ϕ_α}, K] [Eq. (6)] with respect to K and the NGWFs. Since $ExcMGGA$ is a functional of τ, but the dependence of τ on n is unknown, evaluation of the functional derivative is not straightforward. onetep’s implementation of meta-GGAs sidesteps this issue by employing the “functional derivatives of τ-dependent functionals with respect to the orbitals” (FDO) method,^58,59 which yields a generalized Kohn–Sham approach⁶⁰ with an XC potential comprised of a local GGA-like part and a non-multiplicative orbital-specific τ-dependent part. This can be expressed in a “per-NGWF” form

which can be easily inserted wherever the product of the XC potential and an NGWF arises (e.g., in Hamiltonian matrix elements and energy gradient expressions).

In onetep v5.2, two meta-GGA XC functionals are available: PKZB^61,62 and B97M-rV.^63,64 The former is an early meta-GGA form (1999), implemented primarily for testing and validation. The latter is a modern and very accurate⁶⁵ semi-empirical functional developed through a combinatorial approach, where the τ-dependent functional form itself was obtained by fitting to empirical data. B97M-rV additionally accounts for dispersion effects via a contribution from the (revised) VV10 non-local correlation functional^66,67 (see also Sec. IV B 2).

Energies computed with meta-GGA XC in onetep using a minimal basis of NGWFs show excellent agreement with large basis set calculations in qchem, and the cost of these calculations scales linearly with system size, N.¹⁹ Figure 1 shows that single point energy calculations on amyloid fibril segments up to 13 696 atoms are only a factor of two to three times more costly for meta-GGAs compared to Perdew–Burke-Ernzerhof (PBE) (a GGA functional). Consequently, even very large systems (∼10⁴ atoms) are accessible with meta-GGA XC in onetep.

Hybrid functionals include some fraction of the Hartree–Fock exchange (HFx) energy, $ExHF$⁠. For example, the widely used Becke three-parameter hybrids^68,69 (e.g., B3LYP⁷⁰) have the form

mixing gradient-corrected exchange and correlation (from $ExLDA$⁠, $ΔExB88$⁠, $EcLDA$⁠, and $ΔEcGGA$⁠) and HFx in fractions determined by the following three fitted parameters: a₀, a_x, and a_c.

Whereas τ(r) is a local quantity, depending on the orbitals at r, HFx is non-local, i.e.,

with an exchange energy density that depends on the form of the orbitals over all space,

The non-locality of HFx presents a challenge to onetep’s linear-scaling computational framework, which is premised on the “nearsightedness” of electron–electron interactions (Sec. I). Nevertheless, onetep is able to perform hybrid functional DFT with $O(N)$ scaling computational cost.⁷¹ This is achieved through a combination of a distance-based exchange interaction cutoff and the use of an auxiliary basis of atom-centered truncated spherical waves (SWs).

The HFx energy can be expressed in terms of the NGWFs and density kernel, i.e.,

with exchange matrix elements

where, for the sake of brevity, we have introduced the notation

The distance-based cutoff simply ensures that the number of non-zero elements in X increases linearly with system size by setting X_αβ = 0 where the distance between the centers of NGWFs ϕ_α and ϕ_β is greater than the cutoff, r_X (Fig. 6).

The SW basis reduces the cost of evaluating the electron repulsion integrals (ERIs) in Eq. (24) in two ways. First, inserting a resolution-of-the-identity in the SW basis (SWRI) into Eq. (24) converts the four-index ERIs into products of more manageable two- and three-index objects, i.e.,

where {f_p} are SWs, V^pq is the inverse of the Coulomb metric matrix in the SW basis, and (ϕ_αϕ_δ|f_p) and (f_q|ϕ_γϕ_β) are electrostatic integrals over NGWFs and SWs. Second, the evaluation of the two- and three-index objects themselves is simplified by the availability of analytic expressions for the Coulomb potentials of the SWs (SWpots).⁷² 

The SWs are products of spherical Bessel functions, j_l(x), and real spherical harmonics, $Zlm(r^)$⁠, localized within spheres of radius a,

SWs are a convenient basis in which to expand quantities used in onetep, since, like NGWFs, they are both strictly localized within a sphere (radius a, set to the NGWF localization radius) and are systematically improvable with respect to a kinetic energy cutoff parameter.⁷³ Since SWs are angular-momentum-resolved, they have also found utility as a basis for onetep’s distributed multipole analysis (DMA, Sec. IV F 1) and projected density of states (p-DOS, Sec. IV E 4) functionality.

Linear-scaling Hartree–Fock exchange in onetep was originally described in 2013.⁷¹ Although the initial implementation of ERI evaluation was not parallelized, $O(N)$ scaling of computational cost was demonstrated for short polyethylene and polyacetylene chains (∼ 10 s of atoms).

Recent developments have accelerated $O(N)$ Hartree–Fock exchange in onetep, making hybrid functional calculations feasible on much larger systems, with ∼ 2000 atoms becoming accessible with routinely available supercomputing hardware (Fig. 7). This dramatic improvement in performance was achieved by resolving bottlenecks in two major components of the implementation. First, a novel mixed numerical–analytic approach was developed for evaluating elements of the Coulomb metric matrix in the SW basis [V in Eq. (26)], reducing the computational cost (memory and execution time) of this step by orders of magnitude.⁷⁵ Second, the parallel implementation of the evaluation of Hartree–Fock exchange was completely overhauled with efficient data distribution among MPI processes, extensive use of OpenMP threading, and intelligent caching of computed quantities in memory.⁷⁶ The excellent parallel scalability of the overhauled implementation is illustrated in Fig. 8.

While great strides have been made in reducing the computational effort associated with hybrid XC calculations in onetep, using semi-local functionals remains considerably less costly. onetep can perform semi-local XC calculations on systems of >10⁴ atoms with routinely available computational resources (see, e.g., Fig. 1). However, applying hybrid XC to systems of this size remains prohibitively expensive. A pragmatic solution to this issue is to employ embedded mean-field theory (EMFT, Sec. IV D 3), wherein hybrid XC may be applied in a small active region (<1000 atoms) that is embedded in a larger environment, treated with a semi-local XC functional.

The release of onetep (v5.2) at the time of writing supports global hybrid functionals in ground-state DFT calculations and conduction NGWF optimization (Sec. IV E 1). The core framework for applying global hybrids in linear-response time-dependent DFT (LR-TDDFT) (Sec. IV E 2; see also Sec. II.C of Ref. 77) has been developed and will be included as a user-accessible feature in an upcoming release. Support for range-separated hybrid functionals in onetep is also planned, though this is currently at an early stage of development and is thus a longer term aspiration.

Properties and processes in many scientifically interesting and technologically important systems depend on weak dispersion interactions, such as London or van der Waals (vdW) interactions. Examples include molecular crystals, physisorption of molecules on surfaces, and inter-layer interactions in layered two-dimensional materials and heterostructures. Large-scale electronic structure methods are well-suited for studying such structurally complex systems; however, standard density functionals such as the LDA and GGA do not account for these non-local and long-range interactions. onetep includes two methods designed to alleviate this problem.

A popular approach to remedying the problem of dispersion interactions is by including a damped London energy expression as part of the total DFT energy expression with perhaps the most well-known of these methods being the Grimme D2⁷⁸ and D3⁷⁹ variants.

In addition to the original D2 empirical dispersion correction of Grimme,⁷⁸ the damping functions of Elstner⁸⁰ and Wu and Yang⁸¹ have been implemented in the onetep code using new parameter sets reoptimized⁸² for biomolecular interactions. This reoptimization has proven critical for the accurate simulation of dispersion interactions, with the interaction energy root mean square error computed at the PBE level improving from 5.915 kcal/mol to 0.813 kcal/mol with respect to the literature values of the 60 complexes used to fit the parameters. Dispersion corrections have also been included in the calculation of forces, enabling accurate structures to be produced through geometry optimization (see Sec. IV C 1).

Dion et al.⁸³ developed a general non-local contribution to the correlation energy functional for describing vdW interactions within the framework of DFT, which takes the form

where ϕ(r₁, r₂) is a kernel that depends on the separation r₁₂ = |r₁ − r₂| and the electronic density ρ(r) at and in the neighborhood of r₁ and r₂ [see Eqs. (14)–(16) of Ref. 83 for the definitions of ϕ(r₁, r₂)]. Equation (28) is a double integral over the whole simulation cell, i.e., over six spatial dimensions, and computing it directly is prohibitive for large-scale systems. Román-Pérez and Soler⁸⁴ noted that the kernel ϕ may be written as

where q₁ and q₂ are the values of a universal function $q0ρ(r),|∇ρ(r)|$ (see Ref. 83 for details) at the points r₁ and r₂, respectively, and the kernel is then represented in terms of a set of interpolating polynomials (cubic spline functions) p_α(q) for a specific set of values of the interpolation points q_α for each r₁₂. In Ref. 84, it was found that ∼20 interpolation points on a logarithmic mesh was sufficient to accurately describe the kernel up to a cut-off value of q₀ ≤ q_c at which point q₀ is smoothly saturated. The benefit of using the expression in Eq. (29) for the kernel ϕ is that Eq. (28) can be evaluated using Fourier methods as follows:

θ_α(r) ≡ ρ(r)p_α[q₀(r)], ϕ_αβ(r) ≡ ϕ(q_α, q_β, r), and $θ̃α(k)$ and $ϕ̃αβ(k)$ are their respective Fourier transforms. In practice, ϕ_αβ(k) is computed and tabulated once and for all on a radial grid, and the main computational effort during a calculation is associated with the evaluation of the Fourier transforms $θ̃α(k)$⁠. Expressions for the non-local potential corresponding to the non-local energy (required for self-consistent calculations and computing Hellmann–Feynman forces) may also be evaluated (see Refs. 84 and 85 for details).

In the case of the original vdW functional of Ref. 83 (known as “vdW-DF”), the total exchange and correlation functional is given by

which is the sum of the exchange term from revPBE,⁸⁶ the correlation term from PW92 LDA,⁸⁷ and the non-local correlation discussed above. In onetep, we have also implemented a number of further vdW density functionals within the same implementational framework, including vdW-DF2,⁸⁸ optPBE-vdW, optB88-vdW and optPBE_κ=1-vdW,⁸⁹ optB86b-vdW,⁹⁰ and C09x-vdW.⁹¹ 

This work has been crucial in enabling applications to van der Waals heterostructures.^92–94 In that work, vdW-DF functionals, particularly optB88,⁸⁹ have been shown to give interlayer distances in good agreement with experimental results for 2D material bilayers and heterostructures.

DFT+U,^95–97 occasionally known as LDA+U or LSDA+U, is a technique that can be used to significantly improve on conventional DFT’s description of strongly correlated systems and also many systems not normally thought of as strongly correlated. Local and semi-local functionals often fail to qualitatively reproduce the physics associated with localized, partially filled electronic subspaces of 3d and 4f character, including the emergence of Mott–Hubbard insulating gaps and magnetic order. Challenging systems in this regard include the first-row transition metals and lanthanoids (and the oxides of both), but also certain magnetic semiconductors, organometallic molecules, and even biomolecules.⁹⁸ Other localized subspaces, including those with 2p and 4d character, or those that are not partially filled, spin-polarized, or spherically symmetrich or atom-centered, are increasingly being successfully targeted using flexible DFT+U implementations.

Several varieties and generalizations of DFT+U are implemented in onetep,^99–101 including the corresponding contributions to the ionic forces and phonon energies. Extensions to conduction band optimization and the excited state regime (TDDFT+U¹⁰²), as well as combinations with other functionalities such as PAW (see Sec. IV A 1), are also implemented. The DFT+U method preserves linear-scaling with respect to system size, and in its conventional fixed-projector form, it introduces only a small increase in the computational pre-factor.⁹⁹ A key strength of the onetep DFT+U code is its modular nature and the ability to test diverse population analysis schemes when defining the DFT+U target subspaces. It provides a convenient starting point for the implementation of constrained DFT¹⁰³ (Sec. IV B 4) and DMFT^104,105 (Sec. IV B 5) in onetep. The onetep DFT+U code has been used for fundamental developments in strong electronic correlation effects and self-interaction error correction,^106–108 and it has been extended to calculate Hubbard U and Hund’s J parameters from first principles¹⁰⁹ and to approximately enforce Koopman’s condition.¹¹⁰ 

DFT+U works by adding a corrective term to the total energy for each subspace selected.^95–97 This is commonly interpreted as a correction for the self-interaction error exhibited by the approximate XC functional,^111,112 applied only to the subspaces where it is thought to be most problematic. Available in onetep is the DFT+U+J energy correction functional,^97,111,113 given by

Here, U^I represents the net, occupancy–occupancy type Hartree + XC self-interaction of subspace I, taking account of the spin σ and screening by the remainder of the system appropriately.^109,J is a counterpart to U that quantifies subspace self-interaction of magnetization–magnetization type. We simply have DFT+U either when J = 0 or when the last, unlike-spin term is neglected (giving an effective U: U_eff = U − J.) Another J term can be optionally added (see Ref. 113), and separate Hubbard parameters for the linear (U₁) and quadratic (U₂) energy terms and spin channels (U^σ) can be used.¹¹⁰ $nmIσm′$⁠, the density matrix of subspace I and spin σ, is defined through the set of non-orthogonal projector functions ${|φmI}$ by

$ρ^σ$ is the global KS density matrix, and O^I is the inverse of the projector overlap matrix $Omm′I=⟨φmI|φm′I⟩$⁠.

Different Hubbard corrections can be applied to atoms of the same chemical species, which is helpful if their chemical environments are not equivalent. The suggested population analysis for DFT+U is based on the pseudo-atomic orbitals used to initialize the NGWFs at the beginning of a onetep calculation. Hydrogenic atomic orbitals of a chosen effective nuclear charge Z may alternatively be used, as can a subset of the optimized NGWFs from a previous calculation, in a projection update procedure, which can be iterated to self-consistency.¹⁰¹ The projection of the local density of states (LDOS) (Sec. IV E 4) onto the union of the DFT+U+J subspaces can also be calculated. With U and J set to zero, the functionality can still prove useful, e.g., for its pseudo-atomic population analysis or for breaking spin and spatial symmetries.^107,109

Constrained density functional theory (cDFT)^114–118 is an established variational^103,119 generalization of DFT that enables practical simulation of charge (spin) excitation energies and transfer, as well as exchange coupling parameters in magnetic systems. cDFT also provides a practical and computationally convenient way of correcting intrinsic deficiencies of standard DFT, such as self-interaction errors (SIEs) and static electronic correlation.^120–122 For an extensive review of the wide range of applications of cDFT, see Ref. 118 and references therein.

cDFT finds ground state solutions of a Kohn–Sham DFT-like Hamiltonian, augmented with an additional (external) potential: the constraining potential(s), V_c. The energy W of a DFT system (of energy E[n]) subject to a constraint (N_c) on the electronic occupation of a suitably defined subspace can be written as¹²³ 

where $ρ^$ is the electronic density operator and $P^$ is the projection operator for the cDFT-subspace(s). Maximization of W with respect to V_c determines the V_c needed to enforce the targeted electronic occupation (N_c) in the cDFT-subspace.^103,119 Atomic forces can then be calculated,^100,123 allowing geometry optimization and molecular dynamics (MD) simulations.

The cDFT implementation in onetep rests on the use of non-orthogonal, atom-centered projectors (|φ_n⟩) to calculate the electronic occupation of the given cDFT-subspace, needed in Eq. (35), as

where O^mn is an element of the matrix O⁻¹ and O_mn = ⟨φ_m|φ_n⟩. The mathematical framework underlying the tensorially invariant approach to cDFT-subspace occupancies, energies, and forces is extensively discussed in Refs. 100 and 123.

In onetep, the cDFT problem is solved via nested conjugated gradients,¹⁴ with the W maximization loop nested between the NGWF (ϕ_α) and density kernel (K) loops. Figure 9(a) illustrates this procedure for a nitrogen molecule, where an electronic occupancy difference of one electron between the two atoms is enforced, demonstrating excellent convergence within 10 NGWF iterations.

In analogy with DFT+U^99–101 (Sec. IV B 3), the cDFT problem can be solved self-consistently by using the NGWFs (ϕ_α) optimized for a given cDFT constraint as cDFT-projectors (φ_n) for a successive cDFT-optimization, as shown in Figs. 9(b) and 9(c), providing rapid convergence, as shown in Fig. 9(d). In a cDFT geometry optimization, a ground state geometry can be obtained initially using cDFT-projectors obtained self-consistently for the initial geometry. New cDFT-projectors can then be found for the current geometry, and the process iterated until W is self-consistent. Figure 9(d) shows that this converges extremely rapidly—within one or two steps. Figures 9(e) and 9(f) report the rate of convergence of this fully self-consistent W and the optimized N–N bond-distance (d_NN), respectively, with respect to the localization radius of the cDFT projectors and NGWFs for two different cDFT-subspaces. W and d_NN are different for the two subspaces for small radii, but as the radii grow, increasing the variational freedom of the problem, the values of W_SC and d_NN for the subspaces converge. This indicates that the self-consistent projector implementation in onetep is relatively insensitive to the choice of the cDFT-subspace representation, demonstrating its robustness.

cDFT in onetep has been applied to photovoltaic materials,^124,125 to interfaces in graphene-based composite materials,¹²³ and to the limitations of cDFT in constraining SIE-free solutions in standard DFT exchange–correlation functionals.¹¹⁰ Going forward, we expect applications of cDFT to large-scale systems (e.g., biomolecules, nanoparticles, nanotubes, etc.) currently not amenable to simulation by standard DFT or beyond DFT approaches.

Dynamical mean field theory (DMFT) is a Green’s function method tailored for treating systems with strong local electronic correlation. Formally, correlation can be defined as the physics/chemistry due to multi-determinantal wavefunctions (that is, beyond Hartree–Fock theory). The effects of correlation encompass dispersion interactions, particle lifetimes, magnetism, satellites, collective excitations, the Kondo effect, Mott insulators, and metal–insulator transitions. Many (but not all) of these phenomena appear in systems containing transition metals or rare-earth atoms whose 3d- or 4f-electron shells are partially filled. Electrons in these shells are in especially close proximity with one another, and thus, their interaction is too pronounced to be adequately described by semi-local DFT, which can even provide qualitatively incorrect descriptions of the electronic structure.

DMFT maps the many-body electronic problem to an Anderson impurity model (AIM) Hamiltonian with a self-consistency condition.^126,127 This model Hamiltonian includes Hubbard and Hund’s-like terms, and local quantum fluctuations are fully taken into account —that is, the self-energy is dynamic but is assumed to be local. Consequently, DMFT can be expected to yield improved results where strong correlation is expected to be (a) important and (b) localized (cf. GW, which can account for non-local self-energies but only includes correlation effects perturbatively). Like DFT+U (Sec. IV B 3), DMFT is typically used to treat localized regions where correlation is important, as defined by some atom-centered Hubbard projectors.¹²⁸ In the case of DFT+U, each correlated subspace is subjected to an additional potential; in DFT+DMFT, these subspaces are subjected to a full Green’s function treatment.

The AIM that the system is mapped on to contains both impurity and bath sites. The impurity sites typically correspond to the 3d/4f orbitals, while the bath sites collectively represent the rest of the physical system. The AIM Hamiltonian is parameterized by inter-site hopping terms, which are determined via a fitting procedure, and Hubbard U and Hund’s J parameters, which are external user-defined parameters that can be determined using linear response or cRPA.^109,111,129 Alternatively, they can be treated as variational parameters.

Once the AIM Hamiltonian has been constructed, it is solved using exact diagonalization, continuous-time quantum Monte Carlo, or some other method, thereby obtaining the impurity Green’s function and self-energy. This is typically the most computationally expensive step in DMFT and places a severe restriction on the total number of sites (impurity plus bath) one can use. The self-energies from each correlated subspace are then “upfolded” using these localized quantities to construct the total self-energy. During this upfolding, one must correct for the exchange and correlation already accounted for at the DFT level. The whole process is then repeated until the total self-energy is converged.

DMFT has proved capable of capturing complex electronic behavior where semi-local DFT fails. For example, DMFT corrects the lattice constants of strontium vanadate, iron oxide, elemental cerium, and delta-phase plutonium, which semi-local DFT substantially underestimates.^130–132 Similar improvements have been observed for spectral properties such as densities of states and optical absorption.¹³³ DMFT also captures important dynamic properties that are enhanced by strong correlation, such as satellite peaks in the photoemission spectra of cerium and nickel.^134,135 It can capture magnetic behavior both above and below the Curie temperature,¹³⁵ has been used to study high-temperature superconductivity in cuprates^136,137 and other materials,^138,139 and can predict mass renormalization in heavy fermion materials.^140–144 

The implementation of DFT + DMFT using onetep involves an interface with the toscam code toolbox.¹⁰⁵ onetep and toscam are responsible for separate sections of the DMFT loop: onetep performs much of the initialization and mappings, while toscam is responsible for solving the AIM. As the calculation proceeds, control alternates between these two programs, with the entire procedure being driven by an overarching script. toscam will be freely available for download in the near future. For more details, refer to Ref. 105. onetep + toscam has already been used to explain the insulating M₁ phase of vanadium dioxide,¹⁰⁴ to demonstrate the importance of Hund’s coupling in the binding energetics of myoglobin,^145,146 and to study the super-exchange mechanism in the dicopper oxo-bridge of hemocyanin and tyrosinase.¹⁴⁷ 

onetep can relax ionic positions with respect to forces either in Cartesian coordinates using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (following the original castep implementation⁴⁹) or in the limited memory l-BFGS version;¹⁴⁸ or in delocalized internal coordinates using the approach of Andzelm et al.¹⁴⁹ 

The forces in onetep are calculated as the sum of local and non-local pseudopotential terms with Ewald and, where applicable, non-linear core correction terms.⁵³ Since onetep’s NGWF basis is atom centered and local, the forces also include Pulay terms due to the change in basis functions with respect to atomic positions.⁵⁴ The total forces in onetep are calculated as

Pre-conditioners are used to speed up the geometry optimizers in onetep. By default, (l-)BFGS uses the physics-based pre-conditioner suggested by Pfrommer,¹⁵⁰ which takes an estimate for the average optical phonon frequency at the Γ-point as input. The algorithm also estimates this quantity throughout the calculation in order that restarts and future calculations may be sped up. onetep’s geometry optimization performance may be improved dramatically¹⁵¹ by using the force-field based pre-conditioners of Packwood et al.¹⁵² These methods take simple interatomic potentials to construct initial sparse approximations to the Hessian matrix. These pre-conditioners are available with the l-BFGS geometry optimizer.

onetep does not perform any extrapolation of NGWFs or density kernels between geometry optimizer iterations; this helps prevent biasing the end result. For computational expediency, however, onetep instead takes the final NGWFs from the previous geometry optimizer iteration as initial NGWFs to the next iteration. The NGWFs are shifted with respect to the simulation grid to keep the ions at their center. On restarted calculations onetep can read the Hessian matrix or update vectors along with the final NGWFs and density kernel from the previous run.

onetep supports constraints on individual atoms in all geometry optimizer modes. Constraints may be to fix an atom in place, to keep it on a vector, or to keep it on a plane. Relaxation of the lattice parameters, however, is not possible within onetep due to the limited demand for the majority of applications of onetep. Although onetep does not calculate stresses for periodic systems, pressure can be applied to finite systems via the electronic enthalpy method, demonstrated in Sec. V C.

The calculation of phonon frequencies is implemented through the well-established finite-displacement methodology (FDM);^153,154 in our case, the dynamical matrix is calculated using a central difference formula with a choice of either a three- or five-point stencil in each Cartesian direction (i.e., 6N/12N displacements in total for an N-atom system).

For crystalline materials, the FDM requires a supercell geometry to capture a sufficiently dense grid of q-points in the vibrational Brillouin zone; this is a natural choice for a linear-scaling code such as onetep. If a crystalline material is specified, only the atoms in the unit cell are displaced, and either a real-space truncation of the force constant matrix or a Slater–Koster style interpolation is used to obtain the phonon frequencies at arbitrary q-points. In order to significantly speed up the calculation, the converged self-consistent density kernel and NGWFs from the unperturbed system are taken as the starting guess for each displacement; it is also possible to limit the degrees of freedom to reduce the number of displacements required and to use a task-farming approach to run the displacements in batches as independent onetep calculations that can be run in parallel.

In addition to the phonon band structure and DOS, the code calculates a number of thermodynamic properties: the zero-point energy, and the free energy, entropy, internal energy, and specific heat in a given temperature range. These are calculated on a regular grid of q-points, which can be chosen independently of the original supercell repetitions. onetep does not calculate the effect of the long-range electric fields associated with longitudinal phonons near the Γ point (LO–TO splitting). Tests on bulk silicon have shown that well-converged quantities can be obtained using a 5 × 5 × 5 cubic supercell (1000 atoms) and a 10 × 10 × 10 q-point grid for the thermodynamic properties, while the phonon DOS requires a significantly denser q-point grid.¹⁵⁵ 

First principles molecular dynamics (MD) simulations constitute a powerful tool to study systems that are inherently dynamical, such as in conformational analysis of (bio)molecules at finite temperatures, bond breaking in chemical reactions, adsorption, and solvation phenomena. MD simulations are also employed to search for equilibrium structures under different external conditions^156–158 as well as for computing vibrational, IR, and Raman spectra¹⁵⁹ (see Sec. IV E 6). Two main paradigms exist to perform MD simulations within the DFT framework: Born–Oppenheimer MD (BOMD) and Car–Parrinello MD (CPMD).^160,161 The Car–Parrinello scheme involves direct propagation of the electronic degrees of freedom along with the nuclear degrees of freedom. On the other hand, in the Born–Oppenheimer MD (BOMD) scheme, one solves the static electronic problem to find the KS energy self-consistently [E_SCF({R_J})] at every MD step. Hence, the electronic wavefunctions are never propagated, since electrons are assumed to instantaneously adjust to the motion of the nuclei¹⁶² (Born–Oppenheimer approximation). This makes BOMD more suitable to linear-scaling approaches since the localization constraints are more easily imposed during the KS energy minimization. The other main advantage of BOMD over CPMD is that a larger time step can be chosen, which, in turn, results in fewer MD steps and reduced simulation times. Moreover, no extra ambiguous parameters, such as fictitious electronic masses, are required to stabilize the dynamics. On the other hand, in BOMD, the stability of the generated trajectory is quite sensitive to the accuracy of the ionic forces and more generally to the degree of convergence of the SCF loop. This is due to the well-known fact that the error in the forces is first order with respect to the error in the incompletely converged wave-function, even though the error in the Kohn–Sham energy is second order.¹⁶³ 

Within the onetep framework, the BOMD equations assume the following form:

where R_I is the position of the Ith ion at time t and M_I is its mass. $FISCF$ is the force acting on ion I (for a detailed description of the computation of ionic forces in onetep, see Ref. 53).

The BOMD equations of motion are integrated using a finite difference scheme (which usually preserves the symplectic structure of phase space)—in onetep, we employ the Velocity-Verlet (VV) algorithm.¹⁶⁴ 

It is possible to couple different thermostats to the BOMD equations in Eq. (38) in order to sample ensembles other than the microcanonical (also known as NVE). In particular, the following thermostats are available in onetep: Andersen,¹⁶⁵ Langevin,¹⁶⁶ Nosé–Hoover chains,¹⁶⁷ Berendsen,¹⁶⁸ and Bussi.¹⁶⁹ 

For each MD step, the dominant computational effort is associated with solving the KS equations self-consistently, whereas the computation of the ionic forces and the integration of the BOMD equations amount to a small fraction of the walltime. The number of SCF cycles required to reach a given level of self-consistency can be substantially reduced by using improved initial guesses for the electronic degrees of freedom. Various approaches have been implemented, which exploit the information included in the MD history to build more efficient initial guesses and speed up the SCF loop. These methods can be grouped into two major categories: extrapolation methods and propagation methods. Various flavors of linear and polynomial extrapolation are available in onetep. Here, the electronic degrees of freedom are represented as polynomials with respect to time t. All these schemes assume that the NGWFs and the density kernel can be treated as independent degrees of freedom and therefore can be independently extrapolated. However, the density kernel can be considered as the representation of the density matrix in the NGWF basis. Therefore, one can also first extrapolate the NGWFs and then project the density kernel onto this new basis. Accordingly, two density kernel projection methods have been implemented: naïve projection and tensorial-compliant projection.¹⁷⁰ The latter takes into account the tensorial nature of the density kernel in the NGWF basis.

For the propagation schemes, we have implemented the approach developed by Niklasson et al. based on an extended Lagrangian (EL) propagation.¹⁷¹ These schemes have shown superior performances in terms of energy conservation compared to simpler schemes (by alleviating the issue of the breaking of time-reversal symmetry). In practice, three flavors of EL propagation can be chosen in onetep: the naïve approach (deprecated), EL with dissipation (dEL),¹⁷¹ and inertial EL with a thermostat (iEL).¹⁷² The latter two methods show similar performance in terms of energy conservation and speed-up of the SCF cycle.¹⁷³ 

It should be stressed that the self-consistent solution is independent of the initial guess only at exact convergence. In practice, some breaking of time-reversal symmetry is to be expected, resulting from memory effects due to the failure of the Hellmann–Feynman theorem, as the electronic system is never exactly converged to its ground state.¹⁶³ 

Given decent initial guesses, geometry optimization (Sec. IV C 1) is a staple in locating useful minimum-energy structures. Transition states, often just as important as the endpoints in chemical reactions, diffusion problems, and solid state transitions, are less straightforward to locate. They are the maximum-energy points on the minimum-energy paths connecting structures. Several approaches exist to approximate these saddle points, with the nudged elastic band (NEB)¹⁷⁴ approach being among the most robust. NEB attempts to minimize several intermediate structures, or “beads,” onto the minimum-energy path from which the maximal point can be approximated or directly measured. The need for a strong transition state searching functionality in linear-scaling DFT packages is highlighted in a recent review of the opportunities and challenges facing this community.¹⁷⁵ To this end, a robust and scalable implementation of NEB has been designed for onetep.

NEB approximates the minimum energy path by connecting several intermediate structures, known as beads, typically initially guessed from linear interpolation, to their neighbors with springs of zero natural length. In order to ensure the beads lie equally distributed and on the minimum energy path, certain components of both the real and spring forces are projected out. This chain-of-states is then evolved until the “nudged” force on each bead is below some tolerance. Given a perfect approximation of the path tangent at each point on the path, the beads will lie on the minimum energy path, regardless of bead count. In practice, the path tangent at each bead’s position is approximated,¹⁷⁶ requiring a few beads to resolve minimum energy paths.

The implementation of NEB in onetep handles each bead in parallel with an internal driver approach. Based on a comparison of NEB minimizers,¹⁷⁷ both FIRE and the chain-global l-BFGS (GL-BFGS) minimizers were selected for the implementation in onetep. A climbing-image mode¹⁷⁸ is available in which the highest-energy bead moves in a modified potential energy surface (PES) in which the saddle point becomes an energy minimum, allowing the bead to climb the path to sample the transition state exactly.

NEB is not a global search—the path it finds depends on the initial guess path. For this reason, significant control is given to the user over the initial definition of the path; an optional intermediate structure can be defined in the input file or each bead’s position can be directly read from an input XYZ file.

onetep’s $O(N)$ formalism and scalable implementation make it an excellent tool for simulating large and complicated systems, such as biomolecules, nanomaterials, and surfaces. The environment in which these types of systems exist—often some kind of liquid solvent—can intimately affect their behavior and properties. Some means of representing solvent environments is therefore vital for realistic and informative electronic structure calculations on these systems.

A direct approach to solvent modeling would be to simply include a number of solvent molecules in the electronic structure simulation of the system of interest. This explicit quantum mechanical approach is superficially very simple, as it does not require the system to be partitioned into solute and solvent regions and treats internal solute and solute–solvent interactions on an equal footing. However, this approach is generally prohibitively expensive, even with an $O(N)$ method. In addition, even if sufficient resources are available to run an electronic structure calculation on the entire solute–solvent system, realistic results require many calculations in order to statistically average over solvent degrees of freedom.

Explicit atomistic solvent models can be made cheaper by treating the solvent environment at a lower level of theory, as in embedded mean-field theory (EMFT, Sec. IV D 3) and QM/MM, Sec. IV D 2). A further simplification is to represent the solvent as an unstructured dielectric medium. In this case, solvent–solvent interactions and solvent degrees of freedom are accounted for implicitly, and we only need to be concerned with the atomistic details of the solute.

In onetep’s implicit solvent model, the solute is embedded in a cavity within the solvent, which is represented as a polarizable dielectric medium. The polarization of the dielectric medium by the solute charge induces a reaction potential that is incorporated into the Kohn–Sham Hamiltonian. The result is a self-consistent reaction field (SCRF) method in which the electrostatic interactions between the solute and implicit solvent are captured within the SCF procedure.

onetep’s implicit solvation model shares the essential features of the continuum dielectric SCRF method with other widely adopted models, such as COnductor-like Screening MOdel (COSMO)¹⁷⁹ and PCM,¹⁸⁰ but differs in the specifics of its implementation. In particular, the model has the following key features:

• The generalized Poisson equation [Eq. (13) with α = −4π, since onetep uses atomic units] is solved for the total electrostatic potential, $v$(r), in 3D using a high-order defect-corrected multigrid method [using the dl_mg library (see Sec. IV A 5)].⁵⁵ 

• The dielectric permittivity, ε(r), is defined in terms of the electron density, n_elec(r), and smoothly varies between vacuum [ε(r) = 1] and bulk solvent [ε(r) = ε_∞] values based on the form proposed by Fattebert and Gygi,¹⁸¹ 

which has only two parameters: n₀, the electron density value at the center of the solute–solvent interface, and β, which determines the width of the solute–solvent transition [between ε(r) = 1 and ε_∞].

• Non-electrostatic contributions to the free energy of solvation are computed by scaling the surface area of the vacuum cavity surrounding the solute, S, by a surface tension, γ, as used by Scherlis et al.¹⁸² to extend Fattebert and Gygi’s electrostatic model, i.e.,

S[n_elec] is evaluated by integration over a thin film centered on n₀,^182,183 and γ_eff is an effective surface tension, scaled to take account of solute–solvent dispersion–repulsion effects.¹⁸⁴,onetep’s minimal parameter implicit solvent model (MPSM) enables the computation of free energies of solvation of small molecules with accuracy comparable to or better than other widely used implicit solvent models, but with only two empirical parameters (e.g., Fig. 10—see Ref. 184 for further details). The model has also proven to be very useful for large-scale DFT calculations, where screening due to the polarizable dielectric medium corrects for spurious electrostatic effects that arise at artificial surface–vacuum interfaces. This effect has been very convincingly demonstrated for proteins and water clusters, where unscreened dipole moments across the system in vacuo can artificially close the HOMO–LUMO gap (e.g., Fig. 11—see Ref. 192).

onetep now also includes the functionality to calculate energies in electrolyte solutions in both open boundary conditions (OBCs) and periodic boundary conditions (PBCs).¹⁹³ The source term in the GPE now also includes the charge density of mobile electrolyte ions, i.e.,

with n_ions[$v$](r) dependent on the electrostatic potential, $v$⁠. The charge density of mobile electrolyte ions can be written in terms of their space-dependent concentration, c_i[$v$](r), as

where z_i is the charge on ion i. The concentrations can be described in terms of the bulk concentrations of ions (⁠$ci∞$⁠) as

where $β=1kBT$ and λ(r) is the electrolyte accessibility function, which smoothly varies from zero in the solute to one in the bulk electrolyte solution and prevents electrolyte charge accumulation near the solute. The expression in OBCs is the well-known Boltzmann concentration of ions. The expression in PBCs includes a renormalization factor in order to conserve the ions. onetep also includes the functionality for a linear approximation of the above concentrations when the electrostatic potential is small or when z_i$v$ ≪ k_BT. The above set of equations is solved self consistently within dl_mg (Sec. IV A 5).

Due to the long-ranged nature of Coulombic interactions, it is often necessary to include hundreds of atoms, if not more, in DFT calculations before the relevant properties are sufficiently converged.¹⁹⁴ Suitable truncation of larger systems to computationally affordable sizes usually proves to be far from trivial.¹⁹² Hybrid (QM/MM) approaches offer one way of addressing this problem, embedding a quantum mechanical (QM) calculation within an environment described at the molecular mechanical (MM) level of theory. This is well justified in systems whose properties of interest are spatially localized, with the quality of obtained results depending on the physical soundness of the interface and the accuracy of the MM description.

onetep implements a fully self-consistent interface to the MM suite tinker,¹⁹⁵ which implements the state-of-the-art polarizable force-field AMOEBA^196,197 and standard fixed-point-charge force-fields, such as GAFF.¹⁹⁸ In the case of AMOEBA, the two subsystems (QM and MM) interact via multipolar electrostatics and are fully mutually polarizable with direct energy minimization encompassing the entire (QM + MM) system. Distributed multipole analysis (DMA,^199–201 Sec. IV F 1) is used to represent the QM charge density in the form of point multipoles (up to a quadrupole) that are compatible with the AMOEBA formalism. In order to prevent unphysical charge spilling from the QM subsystem onto the point charges, which is an artifact of the point multipole model, damping of electrostatic interactions is used. Additionally, we employ a simple model for Pauli repulsion between QM and MM subsystems, yielding a suitable repulsive potential that ensures a physically sound localization of the QM charge density.

Our implementation^202,203 is the only one to date enabling in situ minimization of localized orbitals in the presence of a polarizable MM embedding and the only one where the cost of the QM calculation is linear-scaling. This model has been demonstrated to yield remarkably good energetics across the QM/MM interface (as compared to fully QM calculations) even for minimal QM subsystems.

In many potential applications of ab initio atomistic modeling, the important physics or chemistry of the system of interest is associated with a small “active region,” which would ideally be described using a high accuracy QM method. However, the rest of the system, or the “environment,” is still significant through its interaction with the active region. QM/MM (Sec. IV D 2) and implicit solvent (Sec. IV D 1) methods are two ways of approaching this problem; however, if the environment also requires QM treatment, a quantum embedding approach is required. Quantum embedding schemes enable an expensive high accuracy method to be used for the active region and a cheaper lower accuracy method for the environment.^204,205

The quantum embedding formalism implemented within onetep is the embedded mean field theory (EMFT) scheme, first proposed by Fornace et al.^206,207 In general, the density matrix or kernel is partitioned into regional sub-blocks,

where A labels the active region and B labels the environment. Each atom and its associated NGWFs are assigned to a region.

For DFT-in-DFT embedding, the only difference between the lower and higher accuracy methods will be in the exchange–correlation functional—E_xc in Eq. (6). The other terms in Eq. (6), which we will collectively call E₀[{ϕ_α}, K] here for brevity, are the same at both levels of theory. In EMFT, the total energy for the system is then written as

where $Exclow$ and $Exchigh$ are the exchange–correlation functionals of the lower and higher accuracy methods, respectively, and n and n_AA are the total and active region electron density, respectively.

The block orthogonalization procedure used in previous implementations of EMFT to avoid charge spilling²⁰⁸ interferes strongly with the ability to optimize the basis functions (the NGWFs) in onetep. To avoid this, all NGWFs are optimized at the lower level of theory and then fixed for the optimization of K within EMFT. Results on a variety of test systems show that the error associated with this approximation is 1% or less of the difference between the full high- and low-level calculations.²⁰⁷ Embedding a hybrid functional calculation within a semi-local DFT environment enables us to extend the use of HFx to periodic structures containing a subsidiary active region. Our implementation of EMFT has been successfully applied to the calculation of excited state and encapsulation energies for a variety of systems consisting of several thousand atoms²⁰⁷ and is capable of performing hybrid-in-semi-local ground state total energy calculations. Forces are also available within the EMFT framework. Future developments will seek to make this compatible with time-dependent DFT (TDDFT) (Sec. IV E 2), enabling the use of time-dependent EMFT (TD-EMFT).²⁰⁹ 

The psinc basis functions employed in onetep are unitary transformations of plane waves (cf. Sec. II). This enables onetep to leverage Fourier transforms for parts of the calculation that are cumbersome in real space, such as evaluating the Hartree potential, but also implies an underlying periodicity of the system. This can become problematic for systems that are not genuinely periodic, particularly when they are not charge-neutral or possess a significant dipole moment—in such cases, non-negligible artificial interactions between the periodic images arise. The traditional approach of vacuum padding and correcting²¹⁰ for the dependence of the energy on the size of the simulation cell has since been improved upon with a number of techniques for implementing open boundary conditions (OBCs) (see Ref. 211 and references therein). In onetep, the artificial periodicity can be abrogated (thus modeling OBCs) using three different approaches: cut-off Coulomb (CC),^212,213 Martyna–Tuckerman (MT),²¹⁴ and direct real-space calculations using a multigrid solver (MG)⁵⁵ to solve the Poisson equation.

The standard periodic approach evaluates the Coulomb integral over all space, assuming an infinite periodic repetition of the simulation cell. This leads to significant errors in isolated systems with net charges or large multipole moments, as the periodic images artificially interact with one another via long-range electrostatics. The cut-off Coulomb approach counters this by only allowing Coulombic interactions within a specified region, thereby emulating the electrostatics of an isolated system with a periodic representation of the charge density.²¹² By selecting an appropriate region, $D$⁠, one can eliminate spurious electrostatic interactions between charges in the periodic replicas and the home cell while also retaining the correct description of the potential between charges in the isolated system,

The modified Coulomb interaction term can then be convolved with the charge density to give the Hartree potential

Here, point r experiences Coulomb potentials from all n(r′) in region $D$⁠, and potentials from charge densities outside this region are set to 0. $VHCC(r)$ is evaluated in reciprocal space through a Fourier transform of (48). The form of the expression for $VHCC(G)$ is determined by the shape of the region $D$⁠, which also reflects the periodicity of the system. The simplest of these is the 3D Coulomb cut-off type, which truncates the Coulomb interaction to a sphere of radius R_c, thereby eliminating periodicity in the system. This can be extended to 2D (slab) and 1D (wire) geometries, maintaining periodic behavior in a desired plane or axis. The modified derivations of $VHCC(G)$ for each cut-off type are given in the work of Rozzi et al.²¹³ and Sohier et al.²¹⁵ 

To ensure sufficient space between the charge densities of the home cell and those of the periodic images, onetep places the unit cell within a larger, padded cell where n(r) = 0 for all r. Calculations of $VHCC(r)$ on a grid are performed by placing the charge density of the unit cell within the padded simulation cell. This is Fourier transformed to give ñ(G) and multiplied with the appropriate values of $v$^CC(G) to give $VHCC(G)$⁠. These values are subsequently reverse Fourier transformed, and the real-space Hartree potential is extracted from the padded grid to give $VHCC(r)$⁠. A similar approach is used for the local pseudopotential energy term, while for the core–core interaction, the usual Ewald sum is replaced with a direct Coulombic sum of interactions computed in real space. For more details about the implementation in onetep, refer to Ref. 211.

This approach retains Fast Fourier Transforms (FFTs) for the calculation of Hartree and local pseudopotential energy terms. The Coulomb operator is still periodic but is modified in such a way that the interactions between neighboring cells are eliminated.

Under PBCs, the real-space representation of a function f(r) is

where $f̃(G)$ is the Fourier transform of f(r), and the reciprocal space vectors G are chosen to be commensurate with the simulation cell. The resultant function f_per(r) has the periodicity of the cell and is a superposition of periodically repeated functions f(r), one in each cell. This is illustrated on the example of a Gaussian function in Fig. 12 (top).

The Martyna–Tuckerman approach²¹⁴ replaces f_per(r) with a different function f_MIC(r),

which is also periodic but is not constructed as a superposition of periodically repeated functions; rather, the minimum image convention (MIC) is used to construct it from a single function residing in the simulation cell while still satisfying periodicity. This is illustrated in the bottom panel of Fig. 12.

A detailed description of the form of $f¯(G)$ required to accomplish the above, and of the treatment of the G = 0 term, can be found in Ref. 211 and in the original paper by Martyna and Tuckerman.²¹⁴ The modified Hartree potential $V¯H(G)$ and the modified local pseudopotential $V¯locps(G)$ are obtained according to Eqs. (22)–(25) of Ref. 211. The core–core interaction is evaluated using a direct Coulombic sum rather than via Ewald summation.

An important consequence of the MIC approach is that the length of the simulation cell must be at least twice that of the molecule in any dimension; otherwise, unphysical interactions occur between the molecule and the nearest periodic image.

In this approach, FFTs are not used in the evaluation of the Hartree, local pseudopotential, and core–core energy terms. Instead, these calculations are performed in real space.

The Hartree potential under OBCs is obtained using a multigrid solver⁵⁵ (also see Sec. IV D 1) to solve the Poisson equation,

on a real-space grid. Defect correction²¹⁶ is used to improve the accuracy of the second-order solution to a higher order (up to 12th). Dirichlet boundary conditions of the form

must be provided on the faces of the simulation cell ∂Ω. To mitigate the computational cost of evaluating Eq. (52), charge coarse-graining and bilinear interpolation are used.

The local pseudopotential term is obtained directly in real space according to Eq. (B4) of Ref. 216. As in our implementation of the Martyna–Tuckerman method, the core–core interaction is evaluated using a direct Coulombic sum rather than via Ewald summation.

One important limitation of this technique is that the multigrid solver used in onetep only supports orthorhombic (cuboid) grids.

As mentioned in Sec. II, while the NGWF basis is optimized to accurately represent the occupied states, the unoccupied states are typically over-estimated in energy and are in some cases missing.^217,218 As described in Ref. 218, a second set of ‘conduction NGWFs’ may therefore be optimized to represent a select number of low energy conduction states using a projection operator-based approach. Aside from being non-self-consistent, the calculation proceeds along similar lines to a valence onetep calculation. The optimized conduction NGWFs are then combined with the valence NGWFs to form a joint NGWF basis. Provided that one uses sufficiently large conduction NGWF radii (typically larger than valence NGWF radii), it is possible to generate a joint NGWF basis that is capable of accurately representing both the occupied states and a given number of (bound) unoccupied states.^218,219 One can use this basis to calculate optical absorption spectra using Fermi’s golden rule,^218,219 which has been employed to study energy transport in the Fenna–Matthews–Olson light-harvesting pigment–protein complex,^220,221 or TDDFT^77,222 (Sec. IV E 2), or to calculate electron energy loss spectra²²³ (EELS) (Sec. IV E 5). Alternatively, the joint NGWF basis may be directly used to calculate band structures.

As described in Ref. 224, two possible approaches to calculating band structure are implemented in onetep. The first approach, denoted the tight-binding (TB) style approach, is related to both TB and Wannier interpolation, while the second approach, denoted the k · p style approach, is related to k · p perturbation theory. In the TB approach, k-dependent phase factors are directly built into the construction of the Hamiltonian matrix at a given k-point, while in the k · p approach, a k-dependent Hamiltonian operator is derived, with the Hamiltonian matrix being the NGWF matrix elements of this operator.²²⁵ For both approaches, the Hamiltonian is constructed in the joint basis at each of the requested k-points and diagonalized to get the corresponding energies. A method for unfolding supercell band structures has also been implemented.

For a well converged NGWF basis, the two methods give very similar band structures. However in cases where the NGWFs are not well converged, e.g., due to too small localization radii or too loose convergence criteria, the k · p approach tends to result in unphysical gaps at the Brillouin zone boundary with corresponding discontinuities in the unfolded band structure.²²⁴ It is therefore generally recommended to use the TB style approach. As an example, Fig. 13 shows the calculated band structure for a carbon nanotube (CNT) using the TB style approach. Although for this system the low energy conduction states are already relatively well represented, it is necessary to use a joint basis in order to accurately reproduce the band structure over the energy range considered.

Many potential applications of linear-scaling electronic structure methods involve the prediction of optical properties of large-scale systems, ranging from nanostructured materials to pigment–protein complexes. Often, as in the prediction of UV–vis absorption spectra, one is mainly interested in a small number of low-lying excited states. These states can be efficiently obtained in the framework of linear-response time-dependent DFT (LR-TDDFT), where the problem of calculating excitation energies is recast into solving for the lowest eigenvalues of an effective two-particle Hamiltonian.^226,227 In onetep, the resulting non-Hermitian eigenvalue problem is solved by exploiting an equivalence of the problem to that of a set of classical harmonic oscillators.²²⁸ 

To achieve linear-scaling computations of excitation energies, in onetep, the transition density matrix ρ^{1}(r, r′) describing the lowest excitation of the system is expanded in terms of valence and conduction NGWFs optimized in a ground-state DFT calculation

where K^{X} and K^{Y} are the excitation and de-excitation contributions to the transition density matrix, respectively. The lowest excitation energy of the system can then be written as a functional of K^{X} and K^{Y}, and higher excitations can be computed by enforcing an effective orthogonality constraint between the transition density matrices of individual excitations. Both the full LR-TDDFT method and the simplified Tamm–Dancoff approximation,²²⁹ corresponding to setting K^{Y} to zero, are implemented in the onetep code.^77,222

Linear-scaling computational effort with system size for calculating a fixed number of low-lying excited states is then achieved by exploiting sparsity in K^{X} and K^{Y}. However, unlike for a ground-state DFT calculation, where the sparsity of the density kernel is defined by a simple cutoff radius, the sparsity patterns of K^{X} and K^{Y} can be chosen in a physically motivated way to restrict the space of possible solutions to the TDDFT eigenvalue problem. By dividing the system into smaller subsystems, and setting K^{X}αβ and K^{Y}αβ to zero when NGWFs χ_α(r) and ϕ_β(r) are associated with different subsystems, it is possible to suppress all contributions to ρ^{1}(r, r′) arising from charge transfer between subsystems^230,231 (see Fig. 14). Since TDDFT systematically underestimates excitation energies of charge-transfer states,²³² especially when using semi-local functionals, the subsystem approach is an efficient way to remove unphysical low-lying excited states. The approach has been used to study excitations of dyes in explicit solvent clusters^231,233,234 and the computation of excitons in multi-chromophore pigment–protein complexes.²³⁰