>[!WARNING] Example of this Warning >``` > WARNING: PSMAXN for non-local potential too small > PSMAXN should be > 13.5545592237991 NTYP= 1 >``` >[!Note] Acknowledgment > I wrote this note because Dan Criveanu (University of Nottingham) asked me about the dangers of ignoring `PSMAXN` warnings, and I realized I didn't have a good answer. While the answer in this note isn't entirely satisfying to me, it has improved my understanding of the PAW method and how VASP implements it under the hood. > > I'd also like to thank Griffin group members for helpful discussions. # Solution Summary This warning is frequently encountered when the cutoffs of the different PAW sets in your `POTCAR` are quite different (i.e., a few are soft, and one or more are much harder). Not an uncommon error in oxides and nitrides since Oxygen/Nitrogen pseudopotentials/PAW sets tend to be harder than many other elements. Sometimes, the calculations will run to completion (possibly suffering from numerical problems), and sometimes, they will just crash down the line with errors. The solution is to decrease your cutoff (I suggest steps of 10 eV) until the warning goes away (see below to figure out the exact number). In many cases, you only need to reduce it slightly, unless it was already absurdly high. The other option is to use softer pseudopotentials for your hardest elements (e.g, switch out `O` for `O_s`, if it doesn't affect your physics/chemistry), or vice versa. There is some old (ca. 2004) advice on the VASP forums about reordering your POTCAR to have the softest element first. I do not believe this works in modern versions of VASP, based on my reading of `pseudo.F`, but I could be wrong. # Explanation > [!Hint] Suggested Reading > This note isn't about the PAW method, so I'm not going to explain it, but you need to be familiar with it if you want to understand how this warning comes about, what it means, and what it can lead to. I strongly recommend Professor Qijing Zheng's excellent [post](http://staff.ustc.edu.cn/~zqj/posts/VASP-All-Electron-WFC/) on obtaining all-electron wave functions from VASP as a starting point if you're at least vaguely familiar with the PAW method and pseudopotential theory. I will use the notation of that note and follow it closely in the preliminaries section, so refer to it if you need clarification or more details. > > Otherwise, you can start with [Blöchl's seminal paper](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.50.17953) on the method. Still, it's not particularly helpful to read if you don't have a basic understanding of pseudopotential theory, the precursor of the PAW method. ## Preliminaries The PAW projectors $\ket{\tilde{p}_j}$ for some given atom $j$ at location $\boldsymbol{\tau}_j$ are given by $ \left\langle\mathbf{r} \mid \tilde{p}_j\right\rangle=\tilde{p}_j\left(\mathbf{r}-\mathbf{R}_L-\boldsymbol{\tau}_j\right) ; \quad\left|\mathbf{r}-\mathbf{R}_L-\boldsymbol{\tau}_j\right|<r_j^c $ where $r^c_j$ is the PAW sphere cutoff radius and $\boldsymbol{R}_L$ is the cell vector of the $L^\text{th}$ cell within the Born-von Karman supercell. Everything that follows assumes we're working in the "home cell" $\boldsymbol{R}_0 = (0,0,0)$ for simplicity. This can be decomposed into radial and angular parts $ \tilde{p}_j(\mathbf{r}-\boldsymbol{\tau}_j)=f_j(|\mathbf{r}-\boldsymbol{\tau}_j|) \cdot Y_{l m}(\widehat{\mathbf{r}-\boldsymbol{\tau}_j}) $ where $f_j(r)$ is a smooth, real function and $Y_{lm}$ are real spherical harmonics in the direction $\widehat{\boldsymbol{r} - \boldsymbol{\tau}_j}$. The choice of real functions is for computational efficiency reasons. We can compute this projector in real space $ \begin{aligned} \mathcal{F}\left[\tilde{p}_j(\mathbf{r}-\boldsymbol{\tau})\right](\mathbf{k}) & =\frac{1}{\sqrt{V}} \int_V \tilde{p}_j(\mathbf{r}-\boldsymbol{\tau}) e^{-i \mathbf{k} \cdot \mathbf{r}} \mathrm{d} \mathbf{r} \\ & =\frac{i^{-\ell}}{\sqrt{V}} \cdot g_j(k) \cdot Y_{l m}(\hat{\mathbf{k}}) \cdot e^{-i \mathbf{k} \cdot \boldsymbol{\tau}} \end{aligned} $ where $ g_j(k)=4 \pi \int_0^{\infty} f_j(r) j_l(k r) r^2 \mathrm{~d} r $ is the spherical Bessel transform of $f_j(r)$, and it is also smooth and real by construction. The `POTCAR` file only contains $g_j(k)$ and $f_j(r)$ and doesn't include the spherical harmonics and phase terms. ## Evaluating the Projectors When computing the nonlocal part of the potential, VASP needs to evaluate quantities like $\braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}}$, where $\ket{\tilde{p}_{j \boldsymbol{k}}} = e^{- \boldsymbol{k} \cdot (r - \tau_j)}\ket{\tilde{p}_{j}(\boldsymbol{r} - \boldsymbol{\tau}_j)}$ and $\ket{\tilde{u}_{n \boldsymbol{k}}}$ is the periodic part of the pseudo-wave. There are two primary ways to evaluate this expression: reciprocal space (`LREAL = .FALSE.`, the default option in VASP) or real space. ### Evaluating Projectors in Reciprocal Space In reciprocal space, $\braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}}$ works out to $ \braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}} =\frac{1}{\sqrt{V}} \sum_{\mathbf{G}} i^l g_j(|\mathbf{G}+\mathbf{k}|) Y_{l m}(\widehat{\mathbf{G}+\mathbf{k}}) e^{i \mathbf{G} \cdot \tau_j} \cdot C_{n k}(\mathbf{G}) \tag{1} $ where $g_j(k)$ is the spherical Bessel transform of the radial part of the projector, $Y_{l m}(\widehat{\mathbf{G}+\mathbf{k}})$ is a real spherical harmonic in the direction of $\boldsymbol{G} + \boldsymbol{k}$, and $C_{nk}(\boldsymbol{G})$ is the Fourier coefficient of the (periodic part of the) pseudo wavefunction. Since $f_j(k)$ is a highly localized, smooth radial function and is always chosen to be zero outside the PAW sphere with a radius $r_j^c$ , then its spherical Bessel transform $g_j(k)$ by construction is also localized and smooth in k-space and drops to zero after some $k = G_j^c$. On the other hand, the pseudo wavefunctions are expanded only up to the cutoff $G_\text{cut} = \sqrt{E_\text{cut}}$ (in Rydberg atomic units, i.e., $G_\text{cut}$ is in $a_0^{-1}$ and $E_\text{cut}$ is `ENCUT` but in Ry) so that $C_{nk}(\boldsymbol{G})$ is identically zero for all $|\boldsymbol{G}| > G_\text{cut}$. Thus, the terms in this sum are essentially zero after $\max\limits_G (G_\text{cut}, G_j^c)$ . Since the PAW method (with a frozen core) is intended to get rid of the high-frequency components of the wave functions in the core so we can use a smaller cutoff, naturally $G_j^c$ is higher than reasonable choices of $G_\text{cut}$ (e.g., the suggested `ENMAX` from the `POTCAR`) under normal circumstances. As a result, the sum normally includes all non-zero plane wave coefficients up to $G_\text{cut}$, which is precisely what we want. However, problems can arise in the uncommon situation where $G_\text{cut} > G_j^c$, where the sum in Eq. (1) would be truncated, and thus, we lose some of the non-zero plane wave coefficients $C_{nk}(\boldsymbol{G})$ with $G_j^c < |G| < G_\text{cut}$. The only possible resolution is to either increase $G_j^c$ (i.e., use a harder PAW data set) or decrease $G_\text{cut}$ (i.e., decrease `ENCUT` in your `INCAR`). This computation is performed in the subroutine `PROJ1` in `nonl.F`. ### Evaluating the Projectors in Real Space In this subsection, I will utilize the very simplistic approach that VASP uses when you set `LREAL = .TRUE.` purely for demonstrative purposes. Note that this is not recommended and can introduce significant errors; see below for further discussion. If we decide to evaluate $\braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}}$ in real space, we end up with $ \braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}} = \frac{V}{N_{\mathrm{fft}}} \sum_{\left|\mathbf{r}_k-\boldsymbol{\tau}_j\right|<r_j^c} e^{i \mathbf{k}_k \cdot\left(\mathbf{r}-\boldsymbol{\tau}_j\right)} \cdot f_j\left(\left|\mathbf{r}_k-\boldsymbol{\tau}_j\right|\right) \cdot Y_{l m}\left(\widehat{\mathbf{r}_k-\boldsymbol{\tau}_j}\right) \cdot \tilde{u}_{n \mathbf{k}}\left(\mathbf{r}_k\right) \cdot $ where $V$ is the volume of the unit cell and $N_\text{fft}$ is the number of points on the FFT grid. As mentioned above, this scheme is unrealistic anyway, but let's pretend for the sake of the argument that the Fourier transform is performed on an FFT mesh inside the PAW sphere, including Fourier components up to $G_j^c$ (which would, in almost all cases, be an extremely dense mesh compared to the coarse wave function grid). Thus, if $G_\text{cut} > G_j^c$, we will end up with wrap-around (aliasing) errors, where high-frequency components of the wave functions are reflected into the smaller sphere of radius $G_j^c$ and represented as lower-frequency components. Aliasing is a well-studied problem with discrete Fourier transforms, leading to well-known DFT errors with the charge density and local potentials. I recommend checking the [VASP wiki entry on aliasing](https://www.vasp.at/wiki/index.php/Wrap-around_errors). >[!Tip] The fine FFT mesh In more realistic approaches to evaluating real space projectors, you would use a second "fine FFT mesh" mesh for this evaluation (usually twice as big in each direction as the coarse FFT mesh used for the wave functions). In most circumstances, this mesh would be too small to include every component of the projectors (i.e., up to $G_j^c$), leading to aliasing errors. See further discussion in the section [[#What happens if I ignore this warning?]] ## Numerical Implementation in VASP The quantity $G_j^c$ is stored internally in `pseudo.F` as `PSMAXN`. Now that we've defined these quantities, we can now generalize the error cited in the [[#`PSMAXN` for non-local potential too small|first section]]: ``` WARNING: PSMAXN for non-local potential too small PSMAXN should be > {QMAXNL} NTYP= {NI} ``` which means that $G_j^c$ for the ion of species `NI` (the ordering in `POSCAR/POTCAR` gives the indices) is larger than `QMAXNL` = $G_\text{cut}$. See `pseudo.F` for more details. (Around line 1269 of the VASP 6.4.1 source code.) `PSMAXN` is calculated for each PAW data set when it is generated and can be found in each `POTCAR` file on the line above the *first* (of multiple, depending on how many nonlocal projectors) occurrence of the words "Non local part," followed by the letter "T". For example, `PSNMAX` (in units of $\mathring{A}^{-1}$) is on line 873 of the `Y_sv` POTCAR, `PSMAXN = 13.41705...` ``` 13.4170545189104 T Non local Part 0 2 1.51485304432685 ``` >[!Hint] What about the local part of the potential? > The local potential has nothing to do with the entire discussion of this note, but I wanted to mention it anyway. The largest $k$ in the local part of the potential is stored in the variable `PSGMAX`, and you can get a similar warning. > ``` > WARNING: PSGMAX for local potential too small > PSGMAX should be > {QMAXL} NTYP= {NI} > ``` > where now `QMAXL` = 2 $G_\text{cut}$ (the factor of 2 because the local potential is evaluated on the fine grid, not the coarse wave function grid). `PSGMAX` is usually quite large; I have never seen this error before. > `PSGMAX` (also in $\mathring{A}^{-1}$) can be found on the line just above `local part` (line 66 of `Y_sv`) > ``` > END of PSCTR-controll parameters > local part > 66.1404096002623 > ``` # The dangers of ignoring this warning The $\braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}}$ term goes into computing the all-electron wave function $ \left|\psi_{n \mathbf{k}}\right\rangle=\left|\tilde{\psi}_{n \mathbf{k}}\right\rangle+\sum_j\left[\left|\phi_j\right\rangle-\left|\tilde{\phi}_j\right\rangle\right] \cdot\left\langle\tilde{p}_j \mid \tilde{\psi}_{n \mathbf{k}}\right\rangle, $ and thus also the charge (and spin) densities. Since all observables in DFT are functionals of the density (whether explicit or implicit), essentially everything computed by VASP will be affected by errors in the evaluation of $\braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}}$. However, it is difficult to say precisely how much of an effect such errors would have and on which quantities without systematic testing. In the case of real space projectors, since high-frequency components get folded in because of the aliasing (see above), you essentially end up with noise in your calculation of $\braket{\tilde{p}_{j \boldsymbol{k}}|\tilde{u}_{n \boldsymbol{k}}}$. For simplicity, the approach discussed in the [[#Evaluating the Projectors in Real Space|section on real space projectors]] is extremely basic (based on `LREAL = .TRUE.`). However, this approach is not recommended precisely due to aliasing errors, and VASP Wiki recommends `LREAL = Auto`. That scheme uses an optimization algorithm (by Georg Kresse) to deal with some of the high-frequency components and minimize aliasing errors. Unfortunately, the scheme is unpublished, so I can't make any in-depth assessments about its efficiency in handling the aliasing errors introduced by the $G_\text{cut} > G_j^c$ situation, but in general, it's supposed to be better than [King-Smith, Payne and Lin](https://doi.org/10.1103/PhysRevB.44.13063). The scheme's code can be found in `optreal.F` in the `OPTREAL_NEW` subroutine, I will update this note if I get the chance to read it in full. >[!Tip] Fun Fact >Georg Kresse is still actively contributing to the VASP code base himself, based on a comment signed `gK` from 2019 in the `optreal.F` source. From a more practical perspective, combining $G_\text{cut} > G_j^c$ with real space, projectors can often lead to various crashes in VASP (e.g., `REAL_OPT: internal ERROR`). Kresse's scheme uses the `INCAR` tag `ROPT` to tune the FFT grid used to compute the real space projectors, and setting `ROPT` to small enough values for the species with softer PAW sets can help resolve crashes and/or numerical errors. I have not seen any crashes when using reciprocal space projectors, but that doesn't necessarily mean proceeding without more careful testing is a good idea. %%==I should probably say something extra here==%% >[!Summary] Conclusion > The safest and easiest way of dealing with this warning is to reduce your cutoff (or change your `POTCAR`). To find how much you need to reduce it by, locate all the `PSMAXN` values in your `POTCAR` (see above), find the smallest one (probably corresponds to the softest element, but not always), and convert it to eV as `ENCUT = PSMAXN**2 * AUTOANG**2 * RYTOEV`, where `AUTOANG` is the Bohr radius to Angstrom conversion factor, and `RYTOEV` is the Ry to eV conversion factor. Any cutoff below this value will eliminate the warning and possible associated errors.