In particular, I'm interested in the general process of selecting values for the planewave cutoff (e_cut) and k-point grid parameters. Although it's very common to read statements like "the calculation should be converged with respect to the number of k-points" in tutorials and documentation, I've spent quite some time investigating this and have found the literature to be surprisingly sparse on rigorous detail.

For the sake of what I'm about to write, I assume that parameters relate to the calculation as described in page 8 of these slides. To be more precise, the choice of k-point grid prescribes the terms of the outer "integral", in fact a summation (over k), and the choice of planewave cutoff e_cut defines the points g used for the inner summation.

As this is ultimately a quadrature approach, I would expect that increasing the fineness of the k-point grid without increasing the planewave cutoff should be "less important", in the sense that no matter how accurate the approximation to the integral of the sum of the planewave terms becomes, it will still be limited in "true" accuracy by the quality of that interior sum. This suggests to me that a sensible approach to obtaining a calculation accurate to some epsilon would be to first converge the planewave cutoff until the delta in total energy falls below that epsilon, and then increase the fineness of the k-point grid, again until the delta falls below the same epsilon. The experiments I've done seem to produce "sensible" results here, but I'm not entirely confident in the method. Does this make sense?

Nevertheless, having done this, I've noticed in particular that the convergence in increasing density of k-point grids seems to be quite dependent on the planewave cutoff. More specifically, calculations with low planewave cutoffs seem to converge more slowly in k-point grid size than do those with high planewave cutoffs. I find this to be a little surprising! Again, my intuition suggests that convergence rates for an increasing-in-fineness sequence of quadrature grids applied to functions (in this case, the sum of the planewave terms) with more-or-less the same smoothness should be very similar.

To illustrate my point, I've attached a dataset run, which generates total energies for bulk silicon for the same increasing sequence of k-point grid sizes (gamma-point centred Monkhorst Pack, points per side 1, 2, 4, 8, 16, 32, 64), taken at two different values of ecut (100 eV and 500 eV). I used the standard JTH PAW pseudopotentials, with double-grid cutoffs taken as twice the planewave cutoff, so 200 eV and 1000 eV respectively. The results fall out as:

Code: Select all

```
etotal11 -7.3274128604E+00
etotal12 -7.8649778973E+00
etotal13 -7.9518906162E+00
etotal14 -7.9582974662E+00
etotal15 -7.9584428547E+00
etotal16 -7.9584504688E+00
etotal17 -7.9584557575E+00
etotal21 -7.3539551058E+00
etotal22 -7.8809892419E+00
etotal23 -7.9666352650E+00
etotal24 -7.9732760505E+00
etotal25 -7.9734171428E+00
etotal26 -7.9734174432E+00
etotal27 -7.9734174422E+00
```

Is anybody able to offer an explanation for this? I have read, for example, the various canonical papers from Monkhorst and Pack and so on, and nothing in them jumps out to me as an immediate reason for this. Perhaps I'm missing something very obvious?

(Please note that I do understand that this calculation is dubious for real-world use, and could be made more efficient through the use of, e.g., different shifts for the k-point grid -- I'm presenting it as a minimal working example.)