How to perform calc to ecutwfc and ecutrho in Quantum ESPRESSO

Question

People I am using a script to calculate the best ecutwfc value for my system, however, what should be the value for ecutrho? My script should test different values for ecutrho too or I do not need to put this variable in my input file during tests?

Answer 1

The short answer is that you should test ecutrho as well, which is why Quantum Espresso (and other plane-wave DFT programs) allow you to set it as a semi-independent parameter. The range of good choices depends on other details of your calculation, in particular the:

• wavefunction cut-off energy, ecutwfc
• pseudopotential method
• exchange-correlation functional

As a rough guide, a reasonable compromise between accuracy and computational efficiency is often found for ecutrho between 4$\times$ecutwfc and 9$\times$ecutwfc.

Wavefunction cut-off energy

The so-called "soft" density, $\rho(r)$, is related to the Kohn-Sham wavefunctions by, $$ \rho(r) = \sum_{bk}f_{bk} \left\vert\psi_{bk}(r)\right\vert^2, \tag{1} $$ where $f_{bk}$ is the band occupancy (and I've ignored spin degeneracy and k-point weights, for simplicity).

Since the density depends on the wavefunctions to the power of two, the density has Fourier components up to twice the wavevector. We may write,

$$ G_\mathrm{cut}^\rho = 2 G_\mathrm{cut}^\psi \tag{2} $$ where the two $G_\mathrm{cut}$ scalars refer to the cut-off reciprocal lattice vectors for the density and wavefunctions.

The cut-off energy depends on the wavevector squared, $$ \begin{align} E_\mathrm{cut}^\psi &\propto \left\vert G_\mathrm{cut}^\psi\right\vert^2\\ E_\mathrm{cut}^\rho &\propto \left\vert G_\mathrm{cut}^\rho\right\vert^2\\ &= \left\vert 2G_\mathrm{cut}^\psi\right\vert^2\\ &= 4\left\vert G_\mathrm{cut}^\psi\right\vert^2\\ &= 4E_\mathrm{cut}^\psi. \tag{3} \end{align} $$ Thus we see that for the "soft" density, we require ecutrho ($E_\mathrm{cut}^\rho$) to be four times ecutwfc ($E_\mathrm{cut}^\psi$) in order to capture all of the contributions from $\psi$.

Short-cuts

At high cut-off energies, the plane-wave coefficients ($c_{Gbk}$) for the plane-waves with the highest kinetic energy are very small ($c_{Gbk}\rightarrow 0$ as $\vert G \vert\rightarrow \infty$). The Fourier components of the density are a convolution of the wavefunction's Fourier components, which means that the large wavevector components of the density go to zero even quicker than the wavefunction's plane-wave coefficients.

For sufficiently high $E_\mathrm{cut}^\psi$, the large $G$ components of the density may be neglected, in which case a value of $E_\mathrm{cut}^\rho$ less than $4E_\mathrm{cut}^\psi$ may be used safely. This is an approximation, but can speed up calculations on smaller systems where the Fourier transforms are a relatively large proportion of the computational effort, and reduce memory because the density grid storage is a significant fraction of the storage required.

Pseudopotential method

For conventional norm-conserving pseudopotentials, the "soft" density is the only contribution to the electronic density, $\rho(r)$; however, modern pseudopotentials often include charge augmentation (ultrasoft, PAW, even modern norm-conserving pseudopotentials) and so,

$$ \rho(r) = \rho^\mathrm{soft}(r) + \rho^\mathrm{aug}(r). \tag{4} $$ In contrast to $\rho^\mathrm{soft}(r)$, $\rho^\mathrm{aug}(r)$ does not principally depend on the wavefunctions $\psi$, but on the pseudoatomic augmentation charges $Q(r)$, which are properties of the pseudopotential. These augmentation charges are actually designed to contain the high-wavevector components of the charge density, and so they typically require significantly finer real-space grids (and hence larger Fourier wavevectors) to represent them accurately.

When performing calculations with these pseudopotentials, $E_\mathrm{cut}^\rho$ is typically a higher multiple of $E_\mathrm{cut}^\psi$ than it would be for conventional norm-conserving potentials. If your simulation has chemical elements with dense semicore states (e.g. $d-$ or $f-$ states), then the charge augmentation functions $Q(r)$ will have very large Fourier components, and you will probably need $E_\mathrm{cut}^\rho$ to be much larger than $4E_\mathrm{cut}^\psi$.

Note that $E_\mathrm{cut}^\rho$ is higher as a multiple of $E_\mathrm{cut}^\psi$, but it is not usually much larger in absolute terms, since the same physical charge density features should be present regardless of the pseudopotential formalism chosen. The advantage of the charge augmentation approach, however, is that by moving the high-wavevector components into $\rho^\mathrm{aug}$, they no longer have to be represented by $\rho^\mathrm{soft}$ and so a lower plane-wave cut-off energy may be used for the wavefunctions.)

Exchange-correlation functional

At first glance, the choice of exchange-correlation functional would not appear to have any bearing on $E_\mathrm{cut}^\rho$, beyond any effect on the pseudopotential augmentation charges (since we use pseudopotentials generated with the same exchange-correlation functional). However, $E_\mathrm{cut}^\rho$ defines the Fourier and real-space grids used not only for the densities, but also the potentials.

If we look at the Hartree potential $\widetilde{V}_H(G)$ in Fourier space, we see that, $$ \widetilde{V}_H(G) \propto \frac{\widetilde{\rho}(G)}{\vert G\vert^2} \tag{5} $$ and so the large $G$ components of the density are actually suppressed (because they are divided by $\vert G\vert^2$), and we could even use a lower cut-off energy for this contribution. (Side-note: the small $G$ components of the density are enhanced, which is a major contributor to charge-sloshing instabilities in SCF methods.)

If we look at the local density approximation to exchange-correlation, we have,

$$ V_{xc}^{LDA}(r) \propto \rho(r)^\frac{1}{3} \tag{6} $$ which also suppresses the contributions for large $G$, since the power is less than 1.

When we consider more sophisticated exchange-correlation functionals, they often include terms with the gradient of the density, $\nabla \rho$. In Fourier space, a real-space derivative multiplies each Fourier coefficient by $i$ times its corresponding wavevector ($G$), i.e. $$ \widetilde{\nabla \rho}(G) \propto iG \widetilde{\rho}(G), \tag{7} $$ where $i$ is the imaginary number. This actually enhances contributions from the coefficients for large $G$-vectors, and means that it is common for GGA exchange-correlation functionals to need larger values for $E_\mathrm{cut}^\rho$ than the corresponding LDA calculations. This is even more pronounced for meta-GGAs, as well as any functionals involving multiple density (or wavefunction) products, such as hybrid functionals.