Fix breaking change to ISMEAR=-5 handling for NKPT<4 (attempt 2)

[Andrew-S-Rosen]

Closes #335.

[esoteric-ephemera]

After some discussion with @Andrew-S-Rosen, we decided on the following:

1. For "tet" errors, always default to Gaussian smearing. We don't want to assume incorrectly why the user is using too low k-point density
2. For "dentet" errors, preserve most of the original logic. That is, if a calculation is using an automatic k-points scheme (KSPACING or the auto-generated grid), simply increase the k-point density slightly.

There's a new function _increase_k_point_density that does this more systematically, and permits a minimum number of k-points (useful for tetrahedron integration where that is 4).

If the calculation uses explicit k-points, or if the k-point density has already been increased, fall back to Gaussian smearing.

3. Split off the true tet errors from what we're now calling ksymm errors per #338. There's not much documentation on the VASP side for these errors, but it seems like the solutions are the same as for the bravais errors.

Only remaining task is to generate a test file for the tet error handler which I'll do soon

[Andrew-S-Rosen]

Thanks, @esoteric-ephemera! A few comments from me in return based on your changes. I fully approve the "tet" error handler changes but have some remaining questions about the "dentet" error handler.

1. On the VASP forum, for the DENTET error it says " (mind that Gamma has to be included in the k-mesh)". Did you want to include that in the swap as well?

2. If min_kpoints=0 then that will be ignored and turned into 1.

custodian/custodian/vasp/handlers.py

Line 1900 in d8a7896

min_kpoints = min_kpoints or 1

3. I fear the automated k-point handling here might be considered by some as "over-engineered" (😅) but I am not inherently opposed. It might warrant specific tests though.

4. For the DENTET error handling, it's worth thinking about what would happen when trying to model a very large structure (e.g. a MOF or zeolite) where you have a large lattice such that 1x1x1 k-points is more than sufficient. Is it going to try to increase the k-points by 20% up to 100 times? If so, that doesn't seem wise to me because that is surely going to be far too expensive to actually run. If we are going to iteratively update the k-points, I think it should be stop at a much smaller number of attempts.

5. For the DENTET error handling, if the max_inc is hit, the k-points are not updated (if I understand correctly). But then no change is made at all. It should be that if max_inc is hit, we instead resort to the ISMEAR/SIGMA swap.

[Andrew-S-Rosen]

@esoteric-ephemera: For the large structure test, you may wish to use https://next-gen.materialsproject.org/materials/mp-1200292 and/or https://next-gen.materialsproject.org/mofs/qmof-d475ee7 as representative examples. The former is a zeolite on the Materials Explorer, and the latter is a MOF from the MOF Explorer.

[esoteric-ephemera]

Thanks @Andrew-S-Rosen!

1. On the VASP forum, for the DENTET error it says " (mind that Gamma has to be included in the k-mesh)". Did you want to include that in the swap as well?

Good point - the first check is now ensuring the grid is Gamma-centered. If it is, then the k-point density is increased at most twice. Then the smearing is changed. Warnings are raised each time

2. If min_kpoints=0 then that will be ignored and turned into 1.

That's intentional - there should always be at least one k-point used in the calculation. Note these are the actual number of k-points used, a zero here doesn't mean automatic grid generation

3. I fear the automated k-point handling here might be considered by some as "over-engineered" (😅) but I am not inherently opposed. It might warrant specific tests though.

Yeah it may be, feedback on this is welcome. It's been moved to a separate custodian.vasp.utils file with tests

4. For the DENTET error handling, it's worth thinking about...

No this just bumps up the k-point density until it finds a grid which has more k-points than the original grid, and has at least the number of min_kpoints. I've added specific tests for mp-1200292 with an input 1x1x1 Gamma-centered grid, and it just bumps this to a 2x2x1 grid

5. For the DENTET error handling, if the max_inc is hit, the k-points are not updated (if I understand correctly). But then no change is made at all. It should be that if max_inc is hit, we instead resort to the ISMEAR/SIGMA swap.

Added as suggested - this is a good fallback

[Andrew-S-Rosen]

Thanks! My only remaining question is if bumping up to 2x2x1 in an example like the zeolite is wise. That's likely a very painful increase in cost, and so the question I'm still wrestling with is it worth it to do this or instead fall back to the ISMEAR change? And if so, where do we define that breakeven point? Do we try to avoid changing gamma-point only calculations at all? (NB: remember that switching to the standard version of VASP not only slows down execution but also means you can't rely on an old WAVECAR anymore and have to do SCF from scratch)

These are tough, subjective questions without a clear answer so let's not let it hold up this PR. But just something to think on.

In any case, your changes make sense to me. Thank you!

[shyuep]

I would suggest a small mod - add a max_natoms as a arg for the handler. This will be used to determine if any supercell operation should be done. A sensible default would be something like max_natoms = 100. If the original cell is already 50 atoms, we will not do a 2x1x1 even and just modify ISMEAR.

[esoteric-ephemera]

Thanks @Andrew-S-Rosen and @shyuep! I've added a max_natoms kwarg to VaspErrorHandler and set it to 50 to disable the k-point density increase

There may be discussion needed as to whether the k-point density increase is worthwhile keeping. There are clear cases (slabs/surfaces, wires) where increasing the k-point density isn't physically a good idea. And many examples that @Andrew-S-Rosen has brought up for practical reasons

It's at least there if we choose to stick with it, but the simplest and best solution may just be to change the smearing and warn the user

[Andrew-S-Rosen]

It's likely worth adding documentation here because it's not immediately obvious to the end-user what this is being used for.

With numpy docstrings, you can do

Attributes
----------
is_monitor
    blah blah
max_natoms
    blah blah blah

Also, I think max_natoms as a name might be misleading because the user might think that they should (for some reason) be specifying the maximum number of atoms across all the systems they intend to run with Custodian when in reality, it is simply a tolerance used elsewhere. Admittedly, I am not sure of a better name (maybe natoms_tol? no...).

[esoteric-ephemera]

Fleshed out the docstr a bit and renamed this to num_sites_kpoint_cutoff. Not the most appealing name but hopefully clearer

[Andrew-S-Rosen]

Thank you, @esoteric-ephemera! I appreciate your hard work on this. I have left one comment but otherwise am happy to defer to you and @shyuep for the remainder of this PR.

[shyuep]

This looks good to me. Is it ready to be merged?

[esoteric-ephemera]

Sorry! I was on vacation without internet access. It's ready to go from my perspective. Just for second opinions: @Andrew-S-Rosen and @mkhorton, do you have any concerns about the new logic/handling?

[Andrew-S-Rosen]

No comments from me, except thank you! Also, good catch about the bug in the VASP manual.