Automated optimization of convergence parameters in plane wave density functional theory calculations via a tensor decomposition-based uncertainty quantification

As a first step towards an automated uncertainty quantification of the total energy surface and derived quantities we start with an analysis of DFT convergence. As model system we consider the energy-volume curve of bulk fcc-Al constructed by changing the lattice constant $a_{lat}$ of the cubic cell at $T=0$ K. Due to crystal symmetry all atomic forces exactly vanish, so that the only degree of freedom is the lattice constant or equivalently the volume of the primitive cell ($V=a_{lat}^{3}/4$). The actual calculations are performed using VASP kresse1993 ; kresse1996 ; kresse19962 . The results and conclusions are not limited to this specific code but can be directly transferred to any plane wave pseudpotential DFT code.

In this study we focus on the two most relevant convergence parameters of a plane wave (PW) DFT code: the PW energy cutoff $\epsilon$ and the k-point sampling $\kappa$. The integer value $\kappa$ defines the Monkhorst-Pack mesh $(\kappa\times\kappa\times\kappa)$ used to construct the k-point set.

To discuss and derive our approach we first compute and analyze the total energy surface $E(V,\epsilon,\kappa)$ as function of volume $V$, energy cutoff $\epsilon$ and k-point sampling $\kappa$. We map the energy surface on an equidistant set of $n_{V}$ volumes $V_{i}$ ranging over an interval $\pm 10\%$ around the equilibrium volume $V_{\rm eq}$ (see Sec. 2.3.2), $n_{\epsilon}$ energy cutoffs $\epsilon_{j}$ and $n_{\kappa}$ k-point samplings $\kappa_{k}$ (see Fig. 1). The generated shifted k-point mesh starts at a minimal k-point sampling of $3\times 3\times 3$ $(\kappa_{min}=3)$ and goes up to a maximum of $91\times 91\times 91$ $(\kappa_{max}=91)$ and an energy cutoff of $\epsilon_{min}=200$ eV up to $\epsilon_{max}=1200$ eV. All calculations are executed for primitive cubic cells with a single atom. The k-point sampling is increased in equidistant steps of $\Delta\kappa=2$ to always exclude the Gamma point. The energy cutoff is increased in equidistant steps of $\Delta\epsilon=20$ eV. For the following analysis we focus on fcc bulk aluminum. Extensive tests for other elements and pseudopotentials show qualitatively the same behavior and will be discussed in Sec. 3.5.

In total $n_{V}\times n_{\epsilon}\times n_{\kappa}=21\times 51\times 45=48195$ DFT calculations have been performed for a single pseudopotential. We note already here that such a large number of DFT calculations is not required for the final algorithm. The complete mapping of cutoff and k-point-sampling is used here only to derive and benchmark the asymptotic behavior for $E(V,\epsilon,\kappa)$ and $B_{eq}(\epsilon_{i},\kappa_{j})\mid_{V_{eq}}$ when going towards large (i.e. extremely well converged) parameters. To analyze the energy surface $E(V,\epsilon,\kappa)$ we define the convergence error with respect to the maximum energy cutoff $\epsilon_{max}$ and k-point sampling $\kappa_{max}$, respectively:

$$\Delta E_{\epsilon}(V,\epsilon,\kappa)=E(V,\epsilon,\kappa)-E(V,\epsilon_{max},\kappa)$$

(1)

$$\Delta E_{\kappa}(V,\epsilon,\kappa)=E(V,\epsilon,\kappa)-E(V,\epsilon,\kappa_{max})\quad.$$

(2)

We will first analyze the dependence of the energy over volume of these quantities since this dependence directly impacts convergence of equilibrium quantities such as bulk modulus, lattice constant etc. The energy volume dependence of the energy cutoff convergence $\Delta E_{\epsilon}$ is shown in Fig. 2(b) for a fixed $\kappa$. As can be seen, using a non-converged energy cutoff $\epsilon_{min}=260$ eV gives rise to a convergence error in the energy that strongly depends on the volume. In the shown example the error is largest for small volumes and monotonously decreases with increasing volume. As a consequence, the equilibrium volume will be shifted to a smaller value compared to the fully converged one.

Fig. 2(b) also reveals a remarkable and highly useful behavior of the energy cutoff convergence $\Delta E_{\epsilon}$: It is in first order independent of the k-point sampling. Taking the difference

$$\Delta\Delta E_{\rm noise}=\Delta E_{\epsilon}(V,\epsilon,\kappa_{2})-\Delta E_{\epsilon}(V,\epsilon,\kappa_{1})$$

(3)

between any two k-point samplings $\kappa_{1}$ and $\kappa_{2}$ for a fixed energy cutoff $\epsilon$ results in a volume dependence that resembles random noise. This is shown exemplary in Fig. 2(d) when computing $\Delta\Delta E$ setting $\kappa_{1}$ and $\kappa_{2}$ to the minimum and maximum cutoff value: The average $\left<\Delta\Delta E(V)\right>_{V}$ is approximately zero, the distribution is Gauss-like and any smooth volume dependence is absent. Due to these characteristics we can therefore regard the variance of this contribution

$$\Delta\Delta E(\epsilon,\kappa)=\sqrt{\left<(\Delta\Delta E_{\rm noise}(V))^{2}\right>_{V}}$$

(4)

as a statistical error. Its origin are discretization errors arising from discontinuous (discrete) jumps whenever a continuous change in $\kappa$ or $\epsilon$ results in a discontinuous change in the integer number of k-vectors or plane waves.

Fig. 2(c) shows for the k-point convergence $\Delta E_{\kappa}(V,\epsilon,\kappa)$ an analogous behavior: It is in first order independent of $\epsilon$. The difference $\Delta\Delta E$, is shown in Fig. 2(d) and by construction (see Appendix A) identical to the one obtained from the cutoff convergence $\Delta E_{\epsilon}(V,\epsilon,\kappa)$. We can therefore conclude that the energy cutoff and k-point convergence can be separately considered, i.e.,

$$\Delta E_{\epsilon}(V,\epsilon,\kappa)\approx\Delta E_{\epsilon}(V,\epsilon,\kappa_{max})\pm\Delta\Delta E(\epsilon,\kappa)$$

(5)

and

$$\Delta E_{\kappa}(V,\epsilon,\kappa)\approx\Delta E_{\kappa}(V,\epsilon_{max},\kappa)\pm\Delta\Delta E(\epsilon,\kappa).$$

(6)

The above equations can be summarized in the following expression of the total energy surface:

$$\begin{split}E(V,\epsilon,\kappa)&\approx E(V,\epsilon_{max},\kappa_{max})\\ &+\Delta E_{\epsilon}(V,\epsilon,\kappa_{max})+\Delta E_{\kappa}(V,\epsilon_{max},\kappa)\\ &\pm\Delta\Delta E(\epsilon,\kappa).\end{split}$$

(7)

The above equation decomposes the original rank-three tensor, which would require to perform $n_{V}\times n_{\epsilon}\times n_{\kappa}$ DFT calculations into three rank-two tensors. The first two require only $n_{V}\times(n_{\epsilon}+n_{\kappa})$ DFT calculation. The last contribution, $\Delta\Delta E(\epsilon,\kappa)$ requires for each pair of $\epsilon,\kappa$ also the computation at all volumes, i.e., $n_{V}\times n_{\epsilon}\times n_{\kappa}$ DFT calculations. In the following Section we will analyze the statistical error and derive a computationally efficient approach to compute this error.

Having the full dataset of DFT energies as function of V, $\epsilon$ and $\kappa$ allows us to directly compute the statistical error as defined in Eq. (4). The results are summarized in Fig. 3(a). In the double-logarithmic plot the k-point convergence shows an almost linear dependence, with a vertical shift when changing the energy cutoff (color coded). The fact that the slope remains unchanged when changing the energy cutoff indicates that the k-point convergence of the statistical error is independent of the energy cutoff except for a proportionality factor. To verify this independence we plot the normalized statistical error $\sigma^{E}(\epsilon,\kappa)/\sigma^{E}(\epsilon,\kappa_{min})$ in Fig. 3(b). Having this insight we consider the normalized statistical error along the energy cutoff, $\sigma^{E}(\epsilon,\kappa)/\sigma^{E}(\epsilon_{min},\kappa)$. The fact that all curves coincide clearly shows that except for a proportionality factor the cutoff dependence is identical.

Using the above identified empirical relations we can approximate the statistical error by:

$$\sigma^{E}(\epsilon,\kappa)\approx\frac{\sigma^{E}(\epsilon,\kappa_{min})\cdot\sigma^{E}(\epsilon_{min},\kappa)}{\sigma^{E}(\epsilon_{min},\kappa_{min})}\quad.$$

(8)

The computation of $\sigma^{E}(\epsilon,\kappa_{min})$ and $\sigma^{E}(\epsilon_{min},\kappa)$ requires $n_{V}\times(n_{\epsilon}+n_{\kappa})$ additional DFT calculations. To obtain a sufficiently large magnitude of the statistical error these calculations are performed at the minimum of the convergence parameter set where the statistical noise is largest. As a consequence these extra calculations are computationally inexpensive. To compute the statistical error a second reference next to $\epsilon_{min}$ or $\kappa_{min}$ is needed (see Eq. 3). We find that $\epsilon_{max}$ and $\kappa_{max}$ provide an accurate estimate. These values do not require any additional DFT calculations since they are identical to the ones used to construct the systematic convergence errors $\Delta E_{\epsilon}$ and $\Delta E_{\kappa}$ in Eq. (7). The above formulation allows a highly efficient computation of the statistical error by reducing the computational effort from $(n_{V}\times n_{\epsilon}\times n_{\kappa})$ DFT calculations to $2n_{V}(n_{\epsilon}+n_{\kappa})$.

Eq. (7) together with Eq. (8) provide a powerful and computationally highly efficient approach to compute energy-volume curves $E(V,\epsilon,\kappa)$ for any set of convergence parameters $\epsilon$ and $\kappa$ and are a key result of the present study. They also form the basis for the automated approach to derive optimum convergence parameters that will be derived in the following. The underlying relations and assumption have been carefully validated by computing and analysing an extensive set of pseudopotentials and chemical elements, as presented below.

The compact representation of the energy surface $E(V,\epsilon,\kappa)$ as function of both physical/materials parameters (i.e. volume $V$) and DFT convergence parameters ($\epsilon$, $\kappa$) by Eq. (7) together with the set of converged hyper-parameters derived in Sec. 2.3.2 allows us to analyze the convergence of important materials parameters such as equilibrium bulk modulus, lattice constant, cohesive energy etc. with a very modest number of DFT calculations. To test the accuracy and predictive power of the approximate energy surface Eq. (7) we first compute the physical quantity from the full set (i.e. $V_{i}$, $\epsilon_{j}$, $\kappa_{k}$) of DFT data. Like in the previous section we focus on the bulk modulus, which is highly sensitive to even small errors.

The deviation between the bulk modulus and its converged value $B_{0}(\epsilon_{max},\kappa_{max})$ as function of the convergence parameters is shown in Fig. 6(a). The color code shows the magnitude of the convergence error in a logarithmic scale. As expected, the error shows a general decrease when going towards higher convergence parameters. The actual dependence, however, is surprisingly complex showing a non-monotonous behavior and several local minima. Performing the same fitting approach on the approximate energy surface Eq. (7) (see Fig. 6(c) gives a convergence behavior that shows the same complexity and is virtually indistinguishable from the one shown in Fig. 6(a). Thus, Eq. (7) provides a highly accurate and computationally efficient approach for uncertainty quantification of DFT convergence parameters.

It may be tempting to identify the local minima in the error surface as optimum convergence parameters that combine low error with low computational effort. Unfortunately, these local minima are a spurious product of an oscillatory convergence behavior, where at the nodal points the value becomes close to the converged result. Since DFT convergence parameters should be robust against perturbations caused e.g. by changing the shape of the cell, by atomic displacements etc. the local minima are likely to shift, merge or disappear. Thus, selecting parameters based on such local minima would make these parameters suitable only for the exact structure for which the uncertainty quantification has been performed. We therefore construct in the following an envelope function that connects the local maxima. The envelope represents the amplitude of the oscillatory convergence behavior and is roughly independent on the exact position of the nodes (phase shifts).

Fig. 6(b) shows the resulting envelope function. It is much smoother than the original error surface (Fig. 6(a) and 6 (c)), decreases monotonously when increasing any of the convergence parameters and is free of any spurious local minima. For the further discussion and interpretation of convergence behavior and errors we will exclusively use the envelope function.

The energy expression given by Eq. (7) consists of two systematic contributions ($\Delta_{\epsilon}$ and $\Delta_{\kappa}$) that smoothly change with the convergence parameters and provide an absolute value. It also includes a statistical contribution $\Delta\Delta E$, which quantifies the magnitude of the fluctuations around this value. This decomposition into systematic and statistical contributions can be directly transferred to the physical quantities derived from the energy surface. To get the systematic contribution only the systematic part of the energy (i.e. the first three terms in Eq. (7)) are used as input for the fit. The statistical contribution is obtained using a Monte Carlo bootstrapping approach: The last term in Eq. (7) is replaced by a normal distribution $N(\mu,\sigma^{2}=\Delta\Delta E)$ and the fitting is performed over a large number of such distributions on the best converged surface $E(V,\epsilon_{max},\kappa_{max})+N(\mu=0,\sigma^{2}=\Delta\Delta E(\epsilon,\kappa)$. The inset in Fig. 3(b) shows the computed propagation of the statistical error in total energy $E$ to the statistical error in the bulk modulus. One can observe a linear relation between the error in the bulk modulus and the magnitude of the noise in the energy-volume curve. This is because the bulk modulus is obtained as the second derivative of the polynomial fitted to the DFT-computed energy-volume curve. Because both fitting and taking the second derivative are linear operations, the resulting bulk modulus is a linear functional of the energy values. Hence, the noise in the bulk modulus scales linearly with the simulated noise in the energy-volume curve.

Fig. 6 shows the error surface for the statistical error (Fig. 6(a) and (d)), the systematic error (Fig. 6(b) and (e)) as well as the total convergence error (Fig. 6(c) and (f)) for Al and Cu. The two elements have been chosen since they represent the two most different cases that we observe for all investigated potentials.

Generally, the statistical error becomes dominant for low convergence parameters (in the bottom-left region) while the systematic error dominates for medium and high convergence parameters. The red solid line in Fig. 6(i) shows the boundary where the magnitude of both errors becomes equal. For Al (Fig. 6(c)) this line is absent since the systematic error dominates over the entire region, i.e., even for the lowest considered convergence parameters. Since we have chosen the minimum at or only slightly below the VASP recommended settings the region captures convergence parameters that are commonly chosen. For Cu, in contrast, the statistical error dominates even when using convergence parameters that are well above the recommended ones, i.e., at lower k-point sampling even for an energy cutoff of more than $600$ eV.

The boundary (red line in Fig. 6(i)) that separates the regions where either the statistical or the systematic error dominate can be interpreted like a boundary in a phase diagram: It directly provides information about the dominating error (phase) for any given set of convergence parameters (state variables). Such diagrams can be thus regarded as error phase diagram. The knowledge of the error type has direct practical consequences. If the systematic error dominates, the convergence error becomes a simple linear superposition of each individual convergence error (see Eq. (7)).¹¹1Strictly speaking Eq. (7) applies only for the total energy. For derived quantities such as the bulk modulus we observe sizeable deviations. The reason is that next to an explicit dependence of e.g. the bulk modulus on the convergence parameters also an implicit one via the equilibrium volume occurs, i.e. $B_{0}(\epsilon,\kappa,V_{0}(\epsilon,\kappa))$ It thus allows for any set of ($\kappa$, $\epsilon$) values to determine the deviation from the converged result including its sign. Due to its additive nature the total systematic error will be always dominated by the least converged parameter.

In contrast, the total statistical error can be reduced to any target by just converging a single convergence parameter, which is a direct consequence of its multiplicative rather than additive nature (see Eq. (8)). It also is the reason why the statistical error decays faster than the systematic error when increasing both convergence parameters simultaneously.

Knowledge of the above introduced and constructed error phase diagram can be directly used to find convergence parameters that minimize computational resources for a given target accuracy. In the region where the statistical error dominates, the multiplicative nature, i.e., where the targeted accuracy can be achieved by converging only a single parameter, which in practice will be the computationally less expensive one. In typical cases this will be the k-point sampling since the necessary CPU time scales linearly with the number of k-points and the number of k-points decreases with increasing system size. In contrast, if highly converged calculations are desired one will be in the region of the error phase diagram where the systematic error dominates. Since in this region the errors of the individual convergence parameters are additive, converging one parameter better than the other would be a waste of CPU time. Thus, the availability of such error phase diagrams will allow us to provide a highly systematic and intuitive way of identifying optimum sets of DFT parameters. In a way, error phase diagrams may become what thermodynamic phase diagrams are for materials engineers today: Roadmaps for identifying optimum paths (convergence parameters) in materials design (DFT calculations).

Using the concepts outlined in the previous sections allows us to construct an easy to implement automated approach. This approach computes for a given chemical element and its pseudopotential a set of optimum convergence parameters. These optimized parameters guarantee that the error for a user-selected quantity (e.g. bulk modulus, lattice constant etc.) is below a user-defined target error. In the present implementation the developed tool accepts only cubic structures where the unit cell can be fully described by the lattice constant as single variable.

The key steps of the automated approach are as follows:

• •

  Determine the approximate lattice constant at $\epsilon_{max}$ and $\kappa_{max}$ using the experimental lattice constant, 21 volume points, a volume interval of $\pm 10\%$, and a polynomial fit of order $d=11$.

• •

  Perform DFT calculations around the computed equilibrium lattice constant using again $21$ volume points and a volume interval $\pm 10\%$. Compute on each of the volume points the energies along ($\epsilon_{min}$, $\kappa_{i}$), ($\epsilon_{max}$, $\kappa_{i}$), ($\epsilon_{i}$, $\kappa_{min}$) and ($\epsilon_{i}$, $\kappa_{max}$).

• •

  Use Eq. (7) to compute the total energy $E$ on the full 3d mesh (V,$\epsilon$,$\kappa$). No extra DFT calculations are needed.

• •

  Compute the systematic and the statistical error of the target quantity $A$ (e.g. bulk modulus) following the discussion given in Sec. 2.3. Construct the envelope function of the combined error.

• •

  Determine on the error surface $A(\epsilon,\kappa)$ the line $\Delta A(\epsilon_{opt},\kappa_{opt})=\Delta A_{target}$, i.e., sets of convergence parameters where the predicted error equals the target error.

• •

  Take the set where the curvature of the line is maximum (i.e., close to a 90^o angle.

The above algorithm has been implemented in the pyiron framework janssen2019 . The pyiron-based module requires as input only the pseudopotential and the selection of the target quantity and error. The setup of the DFT jobs, submission on a compute cluster, analysis etc. is done fully automatically without any user intervention.

To demonstrate the robustness and the accuracy of the proposed approach we have employed our automated tool to compute error surface and optimized convergence parameters for a wide range of chemical elements and pseudopotentials. Fig. 7 shows the convergence behavior for $9$ elements using as target quantity the bulk modulus. For each calculation the PBE-GGA pseudo potential recommended by the VASP manual VASP_manual is chosen. These are also the same pseudo potentials used in the Delta project lejaeghere2016 . To visualize the convergence contour lines at target errors of $0.1$ to $5$ GPa are shown on the error surface. Some of the elements such as e.g. Ca allow a convergence down to $0.1$ GPa even for rather modest convergence parameters. Others, such as e.g. Cu achieve this lower limit barely at the maximum set of convergence parameters studied here. In general, systematic convergence trends that would allow to extract some simple rules for finding optimum convergence parameters are lacking. This emphasizes the importance to provide automatized tools for this task.

For our minimum parameter set the statistical error dominates for all elements except for Al, Ca and Pb as indicated by the grey line in Fig. 7, which marks the boundary where the statistical error and systematic error are equal. Going to errors below 1 GPa the systematic error becomes the dominating contribution except for elements with a very high bulk modulus such as Ir and Pt. In the region dominated by the systematic error the cutoff and k-point related error contribution are additive (Eq. (7)). As consequence, the contour lines in Fig. 7 are either parallel to the $\kappa$-axis (i.e. the dominating error is due to the energy cutoff $\epsilon$) or to the $\epsilon$-axis (with $\kappa$ causing the leading error).

To ’benchmark’ the choices made by our automated tool against the choices made by human experts we include the parameters used in two large and well-established high-throughput studies: The Materials Project jain2011 ; jain2013 ; jain2015 (orange squares) and the delta project lejaeghere2016 (magenta squares). The delta project, which aims at high precision to allow a comparison between different DFT codes, systematically shows an error between $1$ and $5$ GPa, with a clear tendency towards the $1$ GPa limit. The Materials Project, where the focus is on computational efficiency and not the highest precision the error is close to the $5$ GPa limit, for several elements (e.g. Ir, Pt, Au) the error becomes $10$ GPa or larger.