Including dispersion interactions in the ONETEP program for linear-scaling density functional theory calculations

While density functional theory (DFT) allows accurate quantum mechanical simulations from first principles in molecules and solids, commonly used exchange-correlation density functionals provide a very incomplete description of dispersion interactions. One way to include such interactions is to augment the DFT energy expression by damped London energy expressions. Several variants of this have been developed for this task, which we discuss and compare in this paper. We have implemented these schemes in the ONETEP program, which is capable of DFT calculations with computational cost that increases linearly with the number of atoms. We have optimized all the parameters involved in our implementation of the dispersion correction, with the aim of simulating biomolecular systems. Our tests show that in cases where dispersion interactions are important this approach produces binding energies and molecular structures of a quality comparable with high-level wavefunction-based approaches.


Introduction
Calculations of properties and processes in materials from 'first principles' quantum mechanical approaches are widely used today in fields as diverse as materials science, biochemistry and engineering (Marzari 2006). These calculations have become established because of their ability to provide an accurate description at the atomic scale, where quantum rather than classical mechanics apply. One of the most widely used first principles methods is density functional theory (DFT; Hohenberg & Kohn 1964) as formulated by Kohn & Sham (1965).
While these calculations have found applications in diverse areas, it would be desirable in many cases to also use them to study nanoscale objects, such as semiconductor nanostructures that have potential applications in electronic devices or entire biological macromolecules. For example, in the case of biological macromolecules such calculations could be used to provide very accurate binding energies as DFT includes, by construction, the electronic changes that take place, such as charge transfer and polarization, which are omitted in the commonly used atomistic force field approaches. However, conventional firstprinciples methods have a computational cost that scales with the third (or greater) power of the number of atoms in the calculation. Owing to this, unfavourable scaling computations with these methods on more than a few hundred atoms are, in general, not feasible, even on modern supercomputers. To overcome this limitation, several novel reformulations of DFT have been developed, which aim to achieve linear-scaling computational cost (Goedecker 1999). The implementation of these reformulations into robust computational approaches has proved very difficult as challenging mathematical, algorithmic and theoretical problems needed to be overcome (Bowler et al. 2008). Nevertheless, intense research efforts spanning more than a decade have allowed satisfactory progress to be made in addressing the challenges involved with linear-scaling DFT and have led to the development of new codes for such calculations by several groups worldwide (Yang 1991;Ordejón et al. 2002;Skylaris et al. 2005;Bowler et al. 2006;Anglada et al. 2008).
Linear-scaling DFT approaches greatly extend the length scales that we can access with DFT calculations, from hundreds to many thousands of atoms, and are constructed with the intention of producing as closely as possible the same results as one would obtain with conventional approaches if it were possible to use them on such length scales. Consequently, the well-documented inability of common DFT functionals to describe dispersion interactions correctly remains. Dispersion (or London 1930) forces between atoms arise owing to the instantaneous dipoles, brought about by the fluctuations in the positions of the electrons. These induce dipoles in a nearby atom or molecule which then, in turn, interact with the dipole on the original atom or molecule. The energy of these attractive interactions can be approximated by the London formula, which, for a pair of atoms, has the following form: where I j is the ionization potential of atom j; a 0 j is its polarizability volume; and r ij is the distance between the atoms i and j. The potential between the two atoms that interact by dispersion forces can be modelled by the Lennard-Jones formula, The attractive r K6 term comes from the London (1930) formula, while the r K12 term represents a repulsive potential. This is required because at a closer range the electrons on each atom repel each other. The form of the repulsive term, however, is mainly a computational convenience as, for example, an exponential form e Kr=r 0 is a more accurate approximation, but also more costly (Lennard-Jones 1931). Dispersion forces are very weak (e.g. binding energies for noble gas atom pairs are less than 0.25 kcal mol K1 ) but can collectively be responsible for determining the geometry of many molecules and solids. Important cases include the stacking interactions between the p-electron systems, such as between graphene sheets in graphite, and base pairs in DNA. In biomolecular simulations, they are often described as 'hydrophobic' interactions and often play an important role in determining the structure and energetics; therefore, they need to be described as well as other non-covalent interactions (such as ion pairs or hydrogen bonding). The difficulty of describing dispersion with DFT is not an intrinsic failure of the theory as the exact exchange-correlation energy functional would be able to describe all such interactions correctly. However, its form is unknown and approximations are required (Kristyán & Pulay 1994). Commonly, these approximations are based on the local electron density and its gradient, and therefore give a poor description of interactions occurring outside the area of electronic overlap, which is the case with dispersion interactions (Pérez-Jordá & Becke 1995;Meijer & Sprik 1996;Zimmerli et al. 2004). DFT does provide an adequate description of the repulsive interactions at a closer range, where the electron densities overlap. The local density approximation will often appear to show binding (in, for example, a noble gas dimer; Pérez-Jordá & Becke 1995;Elstner et al. 2001;Wu & Yang 2002) but this binding is spurious as it results from the exchange part of the functional, whereas dispersion is a dynamical correlation effect (Kristyán & Pulay 1994;Rydberg et al. 2003). Gradient corrected functionals will usually show no binding at all, although basis set superposition error may often give rise to the appearance of weak binding (Jurečka et al. 2007).
Density functionals capable of explicitly including dispersion interactions are being developed by several groups (Langreth et al. 2004;Sato et al. 2005a,b;Zhang & Salahub 2007). While these functionals are promising, they have a considerably higher computational cost than conventional DFT functionals as they include non-local terms and their description of binding due to dispersion is not yet consistently comparable with high-level wavefunction-based methods, such as coupled-cluster approaches. More pragmatic efforts to improve the treatment of dispersion in DFT have instead focused on empirical corrections, such as the inclusion of a damped London term in the total energy expression. Following the pioneering application of such approaches by Böhm & Ahlrichs (1982) in Hartree-Fock calculations, which lack dispersion by definition, such schemes are implemented by summing the attractions between all distinct pairs of atoms, where f damp (r ij ) is a damping function that decays to 0 for small r ij and is 1 at large distances. This damping function is required because electronic structure calculations provide an adequate description of short-range attractions, and therefore the empirical correction becomes superfluous at small distances. If a damping function is not applied to the dispersion term then the total energy will be distorted, because of the resulting significant artificial strengthening of every covalent bond.
In this work, we describe the implementation of approaches for empirical dispersion corrections in the ONETEP linear-scaling DFT code, including the optimization of their parameters for use in biomolecular simulations. In §2a, a summary of the ONETEP method is given, and its connections with the conventional plane wave pseudopotential DFT approach are highlighted. Section 2b describes the dispersion correction schemes we have explored, and the procedure we have used for the optimization of the parameters follows in §2c. In §3, we present tests and comparisons of the methods on a variety of systems demonstrating that they produce dramatic improvements in binding energies and molecular structures in cases where dispersion interactions are important. Section 4 concludes the paper.

Theory
(a ) The ONETEP program ONETEP (Skylaris et al. 2005) is a linear-scaling approach for DFT calculations, which is based on the reformulation of DFT in terms of the one-particle density matrix. In terms of Kohn-Sham orbitals, the density matrix is represented as where j n (r) is a Kohn-Sham orbital and f n is its occupancy. An equivalent representation is where {f a (r)} are localized functions (Hernández et al. 1996) and K ab , which is called the density kernel, is the representation of f n in the duals of these functions. Most commonly in linear-scaling approaches the density kernel is optimized while keeping {f a (r)} fixed in some suitable form (e.g. pseudoatomic orbitals). Linear scaling is achieved by truncating the density kernel, thus exploiting the exponential decay of the density matrix (in non-metallic systems; Kohn 1996). A particular characteristic of ONETEP is that the localized functions {f a (r)} are also optimized during the calculation, subject to a localization constraint, and are thus known as non-orthogonal generalized Wannier functions (NGWFs; Skylaris et al. 2002). The NGWFs are expanded in a basis set of periodic sinc (psinc) functions (Mostofi et al. 2003), which are equivalent to a plane wave basis as they are related by a unitary transformation. The fact that the NGWFs are optimized in situ allows us to achieve plane wave accuracy with only a minimal number of NGWFs (and hence the smallest possible sparse matrices); furthermore, as our basis set is independent of atomic positions and provides a uniform description of space, ONETEP calculations are not affected by basis set superposition error ). The code is parallelized and allows calculations to be performed on large systems containing thousands of atoms .

(b ) Dispersion correction
The various dispersion correction schemes available differ in the form of the damping function f damp (r ij ) that they employ. Two major forms for this function have been widely used. The first form is the one introduced by Mooij et al. (1999) and later generalized by Elstner et al. (2001), ð2:3Þ and the second is a Fermi-like function introduced by Wu & Yang (2002), where c damp is a damping constant and R 0,ij is determined by the range of the overlap of atoms i and j (Elstner et al. 2001 (2.3), and this combination will henceforth be referred to as damping function 2 (DF2). Wu and Yang's damping function in equation (2.4) will be labelled in what follows as damping function 3 (DF3). It is possible to obtain heteroatomic R 0,ij from the homoatomic values using the following expression (Elstner et al. 2001): The homoatomic R 0,i can be estimated from atomic van der Waals radii (Mooij et al. 1999;Wu & Yang 2002). Figure 1 shows the three damping functions; while they have similar shapes, the range of r for which each damping function has values in the interval [0.01,0.99] varies. DF2 has a notably more gentle decay to zero, while DF3 decays particularly abruptly.
The C 6,ij coefficients can be calculated from the homoatomic C 6,i coefficients, which can, in turn, be obtained from experimental work or calculated from the atomic polarizabilities a i using the expression where N eff,i is the effective number of electrons (Elstner et al. 2001). Homoatomic C 6,i coefficients can be combined to give heteroatomic C 6,ij coefficients using one of the two equivalent forms (Elstner et al. 2001;Wu & Yang 2002) of the Slater-Kirkwood combination rule (Slater & Kirkwood 1931), Our aim is to use the empirical dispersion correction schemes to improve the description of biomolecular systems in large-scale DFT calculations. We have therefore implemented and tested in ONETEP the schemes described in §2b. Our approach also involves the optimization of the parameters involved for each exchange-correlation functional in ONETEP. In order to optimize the parameters, we required a benchmark set of complexes with dispersion interactions where the binding energies are known for high accuracy. Subsets of the JSCH-2005 and S22 sets (Jurečka et al. 2006) were chosen for this task. The S22 set is a set of 22 complexes designed to be used as a training set for the inclusion of dispersion corrections and consists of seven hydrogen-bonded complexes, eight complexes with predominant dispersion contribution and seven complexes with significant contributions from both dispersion and hydrogen bonding to the binding. The reference binding energies have been calculated by a combination of MP2 and CCSD(T) methods and extrapolated to the complete basis set limit of CCSD(T). The geometries of the S22 set were obtained by geometry optimizations using MP2 (using a cc-pVTZ basis set and applying a counterpoise correction) for the larger complexes and CCSD(T) (using cc-pVTZ or cc-pVQZ basis sets) for the smaller complexes (Jurečka et al. 2006). The JSCH-2005 set provides similarly high-quality binding energies and reference geometries for sets of base pairs and amino acid pairs. The geometries of the complexes in the subset of the JSCH-2005 set that is used in this work were obtained by hydrogen-only geometry optimizations of geometries obtained experimentally. A subset of the stacked base pairs and amino acid pairs from the JSCH-2005 set and the non-hydrogen-bonded complexes from the S22 set were used for our benchmarks. In addition to these, six sulphur-containing complexes from Morgado et al. (2007), with binding energies calculated predominantly by MP2, were included so that the parameters for sulphur could also be optimized. The geometries of these complexes were obtained by BLYP-D (with a TZV basis set) optimization (Morgado et al. 2007). In total, 60 complexes were chosen for the optimization of 11 parameters. The inclusion of further base pairs from the JSCH-2005 set was deemed undesirable as this could unbalance our chosen training set by giving a bias to base pairs. Also, the hydrogen-bonded complexes in the above sets were omitted, as the empirical dispersion corrections are not designed to describe hydrogen bonding; therefore, optimizing the parameters in a way that causes them to do so (by imitating the CCSD(T) description of the hydrogen bonds) could compromise their description of dispersion interactions.
Binding energies for the chosen set of complexes were obtained with all of the currently available GGA functionals in ONETEP: PBE (Perdew et al. 1996); PW91 (Perdew 1991); revPBE (Zhang & Yang 1998); and RPBE (Hammer et al. 1999), as differences of the single-point energies of the bound complex and the two monomers. The geometries were used as provided by the literature and were not modified in these calculations. Subtracting from the reference 'exact' binding energies gave the error in the binding energy (and ideal dispersion energy correction) for each complex. The goal of the optimization was to adjust the parameters in the dispersion formula (1.4) to minimize the difference between the value of the dispersion energy and the error in the binding energy for each complex. The parameters optimized were the C 6,i coefficients, the R 0,i and the c damp coefficients. Our optimization strategy involved the minimization of the following object function: where the index A runs over the complexes; DE disp,A is the current dispersion energy contribution; E bind uncorr;A is the pure ONETEP DFT binding energy (without dispersion); and E bind lit;A is the literature CCSD(T) or MP2 binding energy. For the optimization, the N eff and initial C 6,i parameters were taken from Wu & Yang (2002) for carbon, hydrogen, nitrogen and oxygen, and from Halgren (1992) for sulphur. The initial c damp parameters used were: 3.0 for DF1, as used by Elstner et al. (2001); 3.54 for DF2, the value Wu & Yang (2002) proposed , rather than Mooij's value of 7.19 (Mooij et al. 1999); and 23.0 for DF3 following Wu & Yang (2002). We used R 0,i values from Elstner et al. (2001). The optimization was considered converged when either of the following two criteria was satisfied.
-The largest change in any parameter from its initial value exceeds 20 per cent.
Since the initial parameters are derived from physical quantities, the optimized parameters should not vary considerably in order to preserve transferability and avoid over-optimization to the fitting set. -An iteration satisfied the following inequality: maximum percentage change in a parameter in the current step percentage change in Err in the current step ! 0:5; ð2:10Þ which ensures that the parameters were varied only when this led to a significant reduction in the object function.
The c damp parameter was not restricted by the former criterion as it is completely empirical; for example, Mooij's proposed c damp for DF2 is double Wu and Yang's proposed value. The parameters for sulphur were further optimized by starting from the parameters obtained with the entire set (of 60 complexes) and optimizing only the sulphur C 6 coefficient and R 0 with the set of sulphurcontaining complexes. In this case, the maximum parameter change was limited to 15 per cent, with the latter convergence criterion the same as above. To eliminate possible effects of the basis set, the single-point energy calculations with ONETEP were performed with a large kinetic energy cut-off of 1200 eV, giving a nearcomplete psinc basis set. Also, large NGWF radii of 8.0a 0 were used for all elements (except hydrogen that had NGWF radii of 7.0a 0 ).

(d ) Atomic forces and geometry optimization
ONETEP is able to compute atomic forces (as analytic derivatives of the total energy) and use these to perform geometry optimizations. We have included in the forces the contribution from the dispersion interactions so that their effect on determining molecular structure can be taken into account during geometry optimizations. As dispersion interactions dominate only in very weakly bound complexes, a very accurate calculation of all the forces is required. This is possible as the NGWFs are essentially expressed in a plane wave basis, and therefore the 'egg box' effect (Tafipolsky & Schmid 2006) of energy variation with respect to the real-space grid, which is typically observed in real-space techniques, is negligible in ONETEP.

Results and discussion
The 60 complexes we used for the fitting of the parameters are presented in table 1, in which the dispersion-including binding energies as obtained with ONETEP are given using the optimized parameters for the three damping functions (DF1, DF2 and DF3) and the PBE (Hammer et al. 1999) functional. Table 1 also contains the ONETEP binding energies that are obtained when no dispersion contribution is included. The binding energies are compared with the accurate ab initio benchmark binding energies for these complexes, which are subsets of the JSCH-2005, S22 and Morgado et al. sets of complexes (Jurečka et al. 2006;Morgado et al. 2007). We can observe from table 1 that the inclusion of the dispersion contribution dramatically improves the binding energies, in most cases leading to an agreement with the literature results that is better than 1 kcal mol K1 . The optimization of the parameters has been necessary to obtain this good agreement as, for example, for DF1 with the PBE functional the value of Err (defined in equation (2.9)) was reduced by 78 per cent in the initial parameter optimization, and in the subsequent sulphur parameter optimization of the value of Err (for the subset of sulphur complexes) was further reduced by 32 per cent. After optimization, DF1 (with PBE) produced binding energies with the lowest root mean square (r.m.s.) difference from the literature values of 0.813 kcal mol K1 . DF3 had a r.m.s. difference only slightly higher, 0.820 kcal mol K1 ; however, DF2 was notably worse with a r.m.s. of 0.926 kcal mol K1 , and a very similar trend was observed for the standard deviations. So, for the PBE functional, DF1 is expected to be the most accurate and consistent.
We argued that a key concern of our approach was to retain the transferability of the parameters. To check if this goal has been achieved, we performed validation calculations on a set of complexes that were not included in the fitting set. These complexes are presented in table 2. They are grouped into four categories: interstrand base pairs; stacked base pairs; hydrogen-bonded base pairs; and other hydrogen-bonded complexes (from the remainder of the S22 set). These systems Table 1. Binding energies (in kcal mol K1 ) for the complexes used in the fitting of the parameters. (The ONETEP results are given with and without dispersion interactions with the optimized parameters for the PBE functional and are compared with the 'exact' values from the literature.) literature (Jurečka et al. (2006) and Morgado et al. (2007)) were chosen as they represent a wide range of typical biomolecular environments and also because accurate binding energies are available for these structures in the literature (Jurečka et al. 2006). For the interstrand base pairs and the stacked base pairs the dispersion interaction is the dominant interaction; our results show the same dramatic improvement in the binding energies as in the complexes of table 1. Furthermore, the level of improvement in the binding energies of the stacked and the interstrand base pairs is similar even though only stacked base pairs were included in the fitting set, indicating the generality of the empirical dispersion correction. For the hydrogen-bonded base pairs and other hydrogen-bonded complexes, where the binding is mainly due to hydrogen bonds, the inclusion of the empirical dispersion contribution is not as successful. In a few cases, such as the water dimer, for example, the uncorrected ONETEP binding energy is already too large, and the dispersion correction leads to further overbinding. DF2 gave a significant overbinding for every hydrogen-bonded complex (r.m.s. difference 5.777 kcal mol K1 ), so this function is less applicable to systems with significant hydrogen bonding, which is the norm for many biological molecules. DF3 performed better than DF1 for all but two of the hydrogen-bonded complexes; the r.m.s. differences were 1.225 and 1.367 kcal mol K1 , respectively. For the nonhydrogen-bonded complexes, all the damping functions produced binding energies We have also investigated the effect of the dispersion contribution on the atomic forces by examining the molecular structures obtained during geometry optimization. We have performed full (unconstrained) geometry optimizations on four systems: a benzene dimer; a methane dimer; a methane-benzene complex; and an indole-benzene complex. All the calculations were performed with DF1 and the PBE functional and a rather tight maximum absolute force convergence threshold of 0.001 E h /a 0 was used as the forces due to dispersion are obviously very weak.
For the case of the benzene dimer, the optimization with dispersion contributions resulted in the equilibrium structure shown in figure 2, where the two benzene molecules are in a conformation with their planes parallel, at a separation of 3.7 Å. This is in close agreement with the value of 3.9 Å that has been obtained with CCSD(T) calculations with a near-complete basis set by Sinnokrot & Sherrill (2004). When we do not include dispersion interactions, the pair of benzene molecules experiences only the repulsive potential of the PBE functional and no binding is observed, but the geometry optimization is simply completed when their separation is 5.1 Å, as at this distance the forces are smaller than the set threshold. For the methane dimer, when dispersion interactions are included ONETEP was able to reproduce the geometry obtained by MP2 calculations using a large Gaussian basis set (cc_pVTZ; Dunning 1989). The final structure obtained with ONETEP is shown in figure 2. The structure has the correct symmetry, and the distance between the hydrogen atoms of the two molecules (3.08 Å) is in agreement with the MP2 value (3.07 Å). When the empirical dispersion contributions are omitted from the ONETEP geometry optimization, the methane molecules end up much further from each other and their orientation is very different from that obtained using the MP2 approach.
In the case of the benzene-methane complex, the benzene-methane distance we obtain after optimization when using our empirical dispersion contribution is 3.15 Å (figure 2), which is in close agreement with the value of 2.98 Å from an accurate MP2 geometry (Jurečka et al. 2006). Omitting the dispersion contribution results in a distance of 3.62 Å. For the indole-benzene complex, the indole-benzene distance increased from 3.40 Å in an MP2-optimized geometry (Jurečka et al. 2006) to 3.52 Å when optimized with our empirical dispersion and 3.92 Å when optimized without it. Clearly, in all cases, the inclusion of the empirical dispersion contribution has significantly improved the geometries obtained.

Conclusions
We have presented an implementation of the empirical dispersion contributions for the ONETEP code, including optimization of parameters for use in biological simulations. Optimization of the parameters significantly improved the obtained binding energies for weakly bound complexes. Further validation calculations, which compare with the literature results from explicitly correlated wavefunction methods, show that the inclusion of dispersion interactions provides an adequate description of the binding energies of weakly bound complexes. We found that the damping functions DF1 and DF3 produced the best binding energies, with DF3 being superior for hydrogen-bonded complexes. Inferior corrected binding energies