noloco: An efficient implementation of van der Waals density functionals based on a Monte-Carlo integration technique.

The treatment of van der Waals interactions in density functional theory is an important field of ongoing research. Among different approaches developed recently to capture these non-local interactions, the van der Waals density functional (vdW-DF) developed in the groups of Langreth and Lundqvist is becoming increasingly popular. It does not rely on empirical parameters, and has been successfully applied to molecules, surface systems, and weakly-bound solids. As the vdW-DF requires the evaluation of a six-dimensional integral, it scales, however, unfavorably with system size. In this work, we present a numerically efficient implementation based on the Monte-Carlo technique for multi-dimensional integration. It can handle different versions of vdW-DF. Applications range from simple dimers to complex structures such as molecular crystals and organic molecules physisorbed on metal surfaces.

Introduction

Density functional theory (DFT) [1] is the most powerful and popular ab initio method for describing structural and electronic properties in a vast variety of materials. The only quantities which require an approximation are the expressions for exchange and correlation in the effective potential and the total-energy functional. The most widely used approximations during the last decades, the local-density approximation (LDA) and different flavors of the generalized gradient approximation (GGA), are based on the well-studied limit of the uniform electron gas. They have been very successful in the description of binding and structural properties of tightly-bound materials such as semiconductors or metals. However, for weakly-bound systems such as graphite, molecular crystals, or molecules physisorbed on surfaces, both types of functionals exhibit severe shortcomings as they cannot capture the long-range tail of the non-local van der Waals (vdW) forces [2,3]. While GGAʼs vastly overestimate equilibrium distances or do not lead to any binding, LDA performs reasonably well in case of overlapping densities, but cannot account for the proper asymptotic behavior. Therefore, the development of density functionals which correctly include vdW interactions is of prime interest.

The theoretical description of vdW interactions, and, in particular, a numerically efficient implementation, has been a challenging task for ab initio calculations. Recently, several approaches have been developed to include them within the framework of useful DFT approximations. They range from simple schemes utilizing a semi-empirical vdW correction term [4] to very sophisticated first-principles methods based on the adiabatic-connection formula [5,6]. Among the latter, the vdW-DF approach by Dion et al. [7] attracted considerable interest when demonstrating applicability for a number of vdW-bound systems, as recently reviewed by Langreth et al. [8]. In fact, it enabled the treatment of many cases where conventional approximate DFT functionals failed badly, including dimers [9,10], polymers [11], molecular crystals [12,13], and molecules adsorbed on surfaces [14-17].

vdW-DF theory provides an expression for the non-local correlation energy in terms of a six-dimensional integral based on the knowledge of the electron density. The numerical evaluation of this integral becomes cumbersome for large unit cells and, in particular, for periodic systems. Consequently, a number of recent publications deal with efficient schemes to evaluate the non-local correlation energy [18-20]. In this work, we present an alternative implementation of vdW-DF which is based on the Monte-Carlo method for multi-dimensional integration. The Monte-Carlo integration technique has a range of advantages over standard quadrature rules in cases where the integration is performed over large data grids, as required for the long-ranged vdW interactions in 2D or 3D periodic systems.

We provide a set of examples such as the Ar dimer and graphite to test our implementation in terms of accuracy and computational effort. Moreover, the program is applied to more challenging cases: We evaluate the cohesive energy of the molecular crystal naphthalene and present the non-local correlation energy of thiophene on Cu(110) as a typical system of an organic molecule/metal surface junction.

Theoretical background

The vdW density functional

The remarkable achievement of the vdW-DF theory is a simple and computationally attractive expression for the non-local correlation energy. This energy is given by an integral which describes the interaction between electron densities at two points r and [7]: The kernel ϕ is a non-local function of the coordinates, the electron densities, and their gradients, which, can be tabulated in advance in terms of the two dimensionless variables, D and δ, defined by where is, by construction [7], a function of the electron density and its gradient at a given point r. When both d and are large, a simplified asymptotic form of can be used with . Thus, the interaction energy has the correct dependence for large separation r.

Multi-dimensional integration

Numerical quadrature algorithms are known to be the best method for low-dimensional integrals [21,22]. However, their efficiency declines when the dimensionality of the integral increases. As an example, an integral of a d-dimensional function f over the hypercube , evaluated with the help of the trapezoidal rule is expressed as Here n is the number of grid points in each dimension and are the weight factors, depending on the integration scheme. To obtain the integral value, evaluations of the integrand function are required. Moreover, with increasing dimension d, the error estimate deteriorates drastically. In Monte-Carlo integration techniques, on the other hand, the error estimate scales as , independent of the number of dimensions. Hence, the Monte-Carlo method becomes more attractive for numerical integration in high dimensions.

The Monte-Carlo method gives the estimate E for the integral on the left-hand side of Eq. (6) using N random samples of the integrand at points [22,23] The law of large numbers ensures that the Monte-Carlo estimate converges to the true value of the integral, i.e., In contrast to quadrature rules, the accuracy of the Monte-Carlo integration has a probabilistic error bound. One can show that where is the standard deviation of the function . In other words, the Monte-Carlo error estimate is on the average and scales as independent of the dimension d. However, the convergence of the integral to the true value is rather slow as the number of sampling points increases. Fortunately, several techniques have been developed to improve the situation and make Monte-Carlo integration more efficient from a practical point of view [22-24]. The most important variance-reducing techniques are stratified sampling, importance sampling, globally adaptive subdivision, and quasi-random sequences. These are widely used in general, and, in particular, in the CUBA library [24] which is utilized in our implementation.

Implementation

Technical details

We proceed describing our implementation to evaluate the 6-dimensional integral in Eq. (1). In order to calculate the integrand one only needs information about the electron density distribution. The three-dimensional total electron density is provided on a regular three-dimensional grid, which can be generated by any DFT code. In the current implementation, the density, the density gradients (optional), and the information on the unit cell parameters are read in a format which is also used for visualization in the XCrySDen package [25].

The kernel has been pre-calculated and tabulated on a fine grid, i.e., points in the ranges and . As the kernel ϕ diverges for small D values, the starting value is chosen . For values , a constant value of is assumed. This choice of parameters has been checked to have no influence on the final result of the integral. For values , the asymptotic form of Eq. (5) is used. In Fig. 1, a plot of vs D for several values of δ is shown, which agrees perfectly with published data [26]. Since the Monte-Carlo integration requires the evaluation of the integrand at random points, we use a trilinear interpolation [27] to compute the values of the kernel, the electron density, and its gradients at those points.

Fig. 1

The vdW-DF non-local correlation kernel ϕ (multiplied by ) as a function of D for various values of δ.

noloco is a FORTRAN program and makes use of the CUBA library for multi-dimensional integration [24]. Among the four different Monte-Carlo integration algorithms implemented in the CUBA library, the routine DIVONNE has been chosen. This routine shows the best results in terms of numerical stability and computational time for our purpose.

Input

The noloco input file is presented in Fig. 2 for the example of the argon dimer. In line 2 one needs to specify the location of the tabulated vdW-DF kernel. The corresponding file kernel.dat can be found in the program source directory. Line 3 contains the information on the flavor of the vdW-DF functional to be used. ‘vdW-DF’ refers to the version by Dion et al. [7], ‘vdW-DF2’ to the latest release by Lee et al. [28], and ‘VV09’ stands for the version by Vydrov and Voorhis [29]. In line 4, the user specifies the size of a supercell which is used for periodic systems to evaluate . The supercell is created as shown in Fig. 3. Here, denote the number of translational images along the three lattice vectors, corresponding to a total number of unit cells

Fig. 2

Input file as used for the calculation of the non-local energy contribution for the argon dimer, Fig. 4.

Fig. 3

Construction of a supercell to evaluate the non-local energy, , for a periodic crystal according to Eq. (1). denote the numbers of translational images along the three lattice vectors. While the variable r is considered in the unit cell (gray area), needs to be taken into account in the whole supercell.

Setting equal to corresponds to the case of an isolated molecule, i.e., . In this case, the molecule must be placed at the center of a primitive cell which should be large enough to prevent wave-function overlap with periodic images. Note that noloco does not take into account the symmetry of the cell since typical systems in which vdW interactions play an important role (e.g., molecular crystals or molecule/surfaces interfaces), generally have rather low symmetry. In case of periodic structures, due to the translational symmetry of the electron density, the kernel, and hence the whole integrand, fulfills the condition . Therefore, the region of integration over the variable r runs over one unit cell only, while the integration region for is the supercell. The ultimate size of the supercell is determined by the convergence of with increasing number of translational images , depending on the individual system under investigation.

The main parameter required for a Monte-Carlo integration is the requested accuracy of the final value of the integral. This value, which serves as the stopping criterion of the integration, can be chosen either as relative (, line 5) or absolute accuracy (, line 6). The integrator tries to find an estimate E (Eq. (7)) for the integral I which fulfills the relation [24].

The CUBA library provides the possibility to specify a variety of methods to obtain the integrand estimate: a Korobov [30] or Sobolʼ [31] quasi-random sampling, a Mersenne–Twister [32] pseudo-random sampling, or the cubature rules of Genz and Malik [33]. Here, we use the quasi-random Korobov [30] sequences choosing the number of sampling points in the order of , which is given in line 7. The maximum total number of the integrand evaluations can be specified in line 8. This parameter is only used if the desired convergence has not been reached already.

Finally, the name of the file containing the structural data and 3D electron density is specified in line 12. It depends on the code and the underlying band-structure method whether the full density (all-electron codes), the pseudo-density (pseudopotential codes), or the quasi-all-electron density (PAW-based codes) are read in. The flag in line 11 specifies whether the density is given in e/bohr3 or . In line 13, the user can specify a file which contains the density gradients on the same grid as the electron density. In case this data is not available, noloco calculates the density gradients using a three-point formula.

Applications

We illustrate our program with applications to several prototypical vdW-bound systems such as the Ar dimer and graphite, which already have been studied with vdW-DF [7,10,14]. In addition, we provide two more involved examples, the molecular crystal naphthalene [12] as well as thiophene on Cu(110) [16]. We test both reliability and computational efficiency of our implementation. The total energy calculations have been carried out using Quantum ESPRESSO [34]. The information regarding computational time refers to a desktop computer with an Opteron processor, 2.4 GHz and 2 GB RAM.

The computational procedure of vdW-DF is a post-scf scheme [7]. In the first step, the crystal or molecular electron density is calculated using the GGA functional according to Perdew, Burke and Ernzerhof (PBE) [35]. In the second step, the non-local correlation energy is evaluated using the self-consistent electron density obtained in the previous step. In the final step, the vdW-DF total energy is determined as

We use the vdW-DF flavor [7] in all the examples presented here. Other variants [28,29] are, however, also implemented in our package.

Argon dimer

A typical vdW-bound system, in which standard xc functionals are known to fail, is the Ar dimer. Since the Ar atom has a fully occupied 3p shell, the binding in the dimer is purely due to dispersion forces. Self-consistent calculations were performed by putting two Ar atoms in a cubic supercell of 30 bohr and using a plane-wave energy cutoff of 60 Ry.

We first consider the convergence of by choosing different values for the parameter which governs the absolute accuracy of the Monte-Carlo integration (line 6 of the input file). The results are presented in Fig. 4 and Table 1 which also contains the CPU time and the numbers of evaluation needed to achieve the corresponding accuracy. For comparison, the total number of evaluations which would be required by a quadrature rule, , is of order , where , , and are the numbers of grid points along the unit cell vectors where the electron density is represented. As we use 150 points in each direction in this example, is of the order 1013 which corresponds to at least three orders of magnitude higher computing time compared to our approach.

Fig. 4

Binding energy of the Ar dimer as a function of the interatomic distance, as obtained with the vdW-DF, for different parameters, , governing the accuracy of the Monte-Carlo integration.

Table 1

Non-local correlation energy of the Ar dimer for different parameters which govern the accuracy of the Monte-Carlo integration. The corresponding CPU times and numbers of integrand evaluations, , are included.

ϵ_abs [Ha]	Ecnl [Ha]	CPU time	N_eval
10−4	0.18098870	1 min 48 s	0.8⋅108
10−5	0.18097679	8 min 41 s	0.4⋅109
10−6	0.18097998	84 min 55 s	0.2⋅1010

Fig. 4 shows the non-local interaction energy for the Ar dimer as a function of the interatomic distance for different parameters . It turns out that a value of is sufficient to obtain converged results. These are in agreement with the corresponding curves from literature [7,36].

Graphite

Another prototypical vdW system is graphite, for which the interaction between adjacent graphene sheets is of vdW type.

For the self-consistent calculations, an energy cutoff of 35 Ry and an k-mesh are used. The in-plane lattice parameter is fixed to the experimental value of , while the inter-plane distance is varied. As mentioned in Section 3, for periodic systems one needs to check the long-range behavior of . The results of such a test is presented in Table 2, where, for simplicity, we only change the lateral supercell size. Note that there are two layers in the unit cell, such that the non-local interaction of the carbon atoms in one plane with those of the other plane is probed. Fortunately, there is almost no dependence of the CPU time on the number of translational images for values larger than (1 1 0). This demonstrates the efficiency of the chosen integration scheme when dealing with periodic systems.

Table 2

Convergence of the non-local correlation energy (in Hartree) of graphite with respect to the number of translational images together with the corresponding CPU time. The parameter has been chosen as .

(n_xn_yn_z)	Ecnl [Ha]	CPU time
(0 0 0)	0.19687	1 min 31 s
(1 1 0)	0.12717	2 min 48 s
(2 2 0)	0.12332	2 min 56 s
(3 3 0)	0.12284	3 min 07 s
(4 4 0)	0.12268	2 min 53 s
(5 5 0)	0.12261	3 min 28 s

The binding energy of the graphene sheets in a graphite crystal as a function of the interlayer distance is presented in Fig. 5. is calculated using the integration parameters and . The curves obtained are fully consistent with those previously reported [14]. For comparison, results for some widely used approximate xc functionals are shown also.

Fig. 5

Binding energy of graphene sheets in graphite as a function of the inter-layer separation. For comparison, we also include results for different xc potentials and experimental data (Exp. 1 [37], Exp. 2 [38]).

Crystalline naphthalene

Naphthalene is the simplest representative of the oligoacene family. Its crystal structure is visualized in the center of Fig. 6. The dependence of the non-local energy on the parameter is displayed in the left panel together with the corresponding CPU time for the molecular crystal as well as the isolated molecule computed within the same unit cell. are fixed to . One can see that a value of is sufficient to obtain the non-local energy within 1 mHa. The calculation at this accuracy level only takes around 5 min on a desktop machine. As the computing time, however, increases exponentially with this input parameter, needs to be chosen carefully in case of such complex systems.

Fig. 6

Convergence tests for naphthalene. The crystal structure (center) is depicted together with the non-local energy and CPU time as a function of the parameter (left) and the supercell size (right). In the left panel, we used . Results are shown for the crystal as well as for the isolated molecule. In the right panel, data are depicted for and .

The right panel is dedicated to the convergence with respect to the number of periodic images. The results for the naphthalene crystal show quick convergence of the non-local energy as a function of the integration volume. It is interesting to note that the characteristic fluctuations of the integral values are of the order of the requested integration accuracy. As in the case of graphite, there is almost no dependence of the computing time on the integration volume.

Thiophene on Cu(110)

As the last example we consider the adsorption of an organic molecule on a metal surface. Thiophene (C4H4S1) on Cu(110) has been studied in detail in Ref. [16] where it turned out that vdW interaction is a substantial contribution to the binding mechanism. We consider a geometry as depicted in the central panel of Fig. 7, with a molecule–metal distance of 2.6 Å. As in the previous case, we test the convergence behavior as a function of the parameter (left panel) and the supercell size (right panel). is required to obtain the non-local energy within 0.5 mHa. is found to be sufficient for the same accuracy. Again, the computing time hardly increases with the number of periodic images, but increases exponentially with . Since the electron density is small in the region between the molecule and the surface, a dense grid is required. That would render a geometry relaxation with standard integration techniques for such systems hardly doable.

Fig. 7

Convergence tests for thiophene on Cu(110). The adsorption geometry (center) is depicted together with the non-local energy and CPU time as a function of the parameter (left) and the supercell size (right). In the left panel, the integration volume is confined to one unit cell. In the right panel, the data correspond to .

Conclusions

We have presented a numerically efficient implementation of the Langreth–Lundqvist vdW-DF which is shown to be very useful, in particular, when dealing with periodic systems. The computational efficiency is achieved through the usage of the Monte-Carlo integration technique as implemented in the CUBA library [24]. noloco is freely distributed and can be downloaded from http://amadm.unileoben.ac.at/software.html or http://exciting-code.org. Our code presented here is carefully tested, has different flavors of vdW density functionals included, and can be used on top of any DFT code. So far, it has been used in combination with Quantum Espresso [12,39], Blöchelʼs PAW code [16], VASP [17], exciting [40], Wien2k, and FHI-AIMS [41] being applied to a variety of vdW bound systems, ranging from molecular dimers to surface systems [16,17], molecular crystals and their surface energies [12] as well as for investigating the role of non-local correlations in the coinage metals [41].