Proposal for new k-point generation scheme

[shyuep]

People who definitely need to look at this code and comment

• @computron @kristinpersson @hautierg @zbwang @johnson1228 @montoyjh @wmdrichards

I will leave this in the branch for three days. After that, I will merge if there are no comments.

Summary

• Today's kpoint generation scheme is highly inefficient, with rounding to integers occurring at multiple steps along the way.
• The major problem is that even for lattices where you have lattice lengths of 13 A or more, the code assigns a 2x2x2 grid with the MPRelaxSet. This makes calculations of such large cells extremely expensive. The following code gives

from pymatgen import Structure
from pymatgen.io.vasp.sets import MPRelaxSet
import numpy as np
for a in np.arange(3, 17, 1):
    s = Structure.from_spacegroup("Pm-3m", np.eye(3) * a, ["Fe"], [[0, 0, 0]])
    v = MPRelaxSet(s)
    kpts = v.kpoints.kpts[0]
    print("For a = %.1f, k = %d x %d x %d" % (a, *kpts))

For a = 3.0, k = 8 x 8 x 8
For a = 4.0, k = 6 x 6 x 6
For a = 5.0, k = 6 x 6 x 6
For a = 6.0, k = 4 x 4 x 4
For a = 7.0, k = 4 x 4 x 4
For a = 8.0, k = 4 x 4 x 4
For a = 9.0, k = 4 x 4 x 4
For a = 10.0, k = 2 x 2 x 2
For a = 11.0, k = 2 x 2 x 2
For a = 12.0, k = 2 x 2 x 2
For a = 13.0, k = 2 x 2 x 2
For a = 14.0, k = 2 x 2 x 2
For a = 15.0, k = 2 x 2 x 2
For a = 16.0, k = 1 x 1 x 1

Obviously, a 2x2x2 grid for a 15 A cell is ludicrous, especially for insulators (this is the main class of materials the MPRelaxSet handles - metals are a separate class we deal with separately). We have grids of 4x4x4 all the way up to a 9A cell.

New proposal

• We will fix this in several areas.
• In the Kpoints.automatic_density_by_vol method, kppa will no longer be rounded to an integer.
• In the Kpoints.automatic_density method, we will retain non-integer num_div, followed by an improved rounding scheme that rounds something like 2.9 x 2.9 x 2.9 to a MH 2x2x2 grid, a 3.2x3.2x3.2 to a MH 4x4x4 grid, and for n = 1 and n > 8, the current scheme of rounding to an odd Gamma grid is retained. Revised code is below.

num_div = [max(mult / l, 1) for l in lengths]
if all([k < 2 for k in num_div]):
return Kpoints(comment, 0, Kpoints.supported_modes.Gamma, [[1, 1, 1]], [0, 0, 0])

is_hexagonal = latt.is_hexagonal()

# VASP documentation recommends to use even grids for n <= 8 and odd
# grids for n > 8.
num_div = [int(round(i - i % 2 + 2 * math.floor(i % 2))) if i <= 8 else int(round(i - i % 2 + 1)) for i in num_div]

The new grids generated using the same test code is as follows:

For a = 3.0, k = 8 x 8 x 8
For a = 4.0, k = 6 x 6 x 6
For a = 5.0, k = 4 x 4 x 4
For a = 6.0, k = 4 x 4 x 4
For a = 7.0, k = 4 x 4 x 4
For a = 8.0, k = 2 x 2 x 2
For a = 9.0, k = 2 x 2 x 2
For a = 10.0, k = 2 x 2 x 2
For a = 11.0, k = 2 x 2 x 2
For a = 12.0, k = 1 x 1 x 1
For a = 13.0, k = 1 x 1 x 1
For a = 14.0, k = 1 x 1 x 1
For a = 15.0, k = 1 x 1 x 1
For a = 16.0, k = 1 x 1 x 1

Which is a lot more reasonable, especially for an 8x8x8 A sized cell (pretty large by most standards) onwards.

This new scheme has been implemented in the kpts branch.

[computron]

Hi Shyue

For the most part I am OK with the change.

The part I would hedge a bit on is the enforcement of even/odd mesh as per VASP documentation. This is unrelated to your change since it was there before, but I just would take the opportunity to revisit it.

Pretend there are two materials with almost the same lattice parameters. One of them, after doing the various operations, shows that we should have 2.9 kpoints per vector. The other shows 3.1 kpoints. These materials are almost the same, yet one is run with 2x2x2 and the other as 4x4x4 (8X as many reducible kpoints!) That kind of huge discontinuity doesn't feel right to me. This is also because I don't fully understand the reasoning behind VASP's "even for n<=8 rule", which makes me question how seriously we should be taking this rule.

[shyamd]

It might be useful to a have a legacy_scheme flag or something to that effect that retains the old generator method until the effect of this new kpt scheme can be determined. I wonder what it will do for things like NMR and Piezo which generally have very stringent convergence criteria.

[shyuep]

@shyamd I would say for convergence, MP has bigger issues to worry about than the kpoints. Most of the problems for NMR and Piezo has to do with proper force convergence. To me, the solution for those is for separately adapted input sets that crank up the density parameter for kpt generation, as well as add the necessary EDIFFG, etc.

[computron]

@shyuep I believe you've also brought up the suggestion in the past to just use VASP's fully automatic scheme, in which a single number is specified and VASP does all the various subdivisions.

https://cms.mpi.univie.ac.at/vasp/vasp/Automatic_k_mesh_generation.html

A couple of notes on this:

• Nothing in their own scheme seems to really enforce the even/odd rule. Thus, they don't seem to be taking that rule too seriously themselves.
• It's possible (although not guaranteed or tested) that the fully automatic method (i.e., single parameter automatic) might have fewer convergence problems than our semi-automatic (3 parameters) scheme. Our justification in the past for doing the latter was that our algorithm would take into account non-orthogonal lattice vectors better. But apart from possibly more stable behavior from VASP, another reason to do fully automatic might be because this generates "gamma-centered Monkhorst Pack" grids which might be better than the original MP scheme for estimating (rough) band gaps.

I guess I don't really have issues with your changes vs what was there before but still want to know if we can do better if we are going to change anyway.

[shyuep]

@computron I have no issues with moving to fully-automatic. My only concern is that such a move requires someone to benchmark the system to determine the proper grid parameter, and whether a robust value can be determined. I seem to recall we had some early attempts to do this benchmarking but can't remember the result.

My current changes is mainly to solve an immediate problem without causing too much impact on MP. For the most part, typical experimental cells are ~4-8 A, which will result in basically a 6x6x6 or 4x4x4 grid. Only at 8A does it transition to a 2x2x2 (which I think is reasonable). The immediate problem is for supercell calculations (e.g., defect studies, NEB calculations, AIMD simulations). There, something like 2x2x2 for a 15A cell is pretty ridiculous and makes what is already an expensive calculation much more expensive.

[computron]

Ok I won't let the perfect be the enemy of the good

But for the future I'd still think about:

• removing even/odd rule
• going fully automatic

[shyuep]

@computron Ok, I have done some very basic testing of the VASP full auto mode. They basically ignore the even/odd rule. So I have removed that rule. Now the corresponding grids are:

For a = 3.0, k = Monkhorst 8 x 8 x 8
For a = 4.0, k = Monkhorst 6 x 6 x 6
For a = 5.0, k = Gamma 5 x 5 x 5
For a = 6.0, k = Monkhorst 4 x 4 x 4
For a = 7.0, k = Gamma 3 x 3 x 3
For a = 8.0, k = Gamma 3 x 3 x 3
For a = 9.0, k = Gamma 3 x 3 x 3
For a = 10.0, k = Monkhorst 2 x 2 x 2
For a = 11.0, k = Monkhorst 2 x 2 x 2
For a = 12.0, k = Gamma 1 x 1 x 1
For a = 13.0, k = Gamma 1 x 1 x 1
For a = 14.0, k = Gamma 1 x 1 x 1
For a = 15.0, k = Gamma 1 x 1 x 1
For a = 16.0, k = Gamma 1 x 1 x 1

The main affected grids are the ones at 7-9A, which now gives a Gamma 3x3x3.

[computron]

Ok - I like the smoothness of the changes better; not sure if switching between Gamma and Monkhorst may lead to issues (hence the VASP fully automatic saying something about "Gamma-centered Monkhorst-Pack" for everything, which I don't fully understand since I never really looked into k-point integration schemes). In any case I am OK with it unless someone else sees a problem.

[shyuep]

@computron In general, the so-called Gamma-Center MonkHorst just simply means you specify Gamma for all grids. Personally, I think the VASP people are just too lazy to figure out whether MH or Gamma and just use Gamma, since that is the most robust scheme (you do not care about whether the cell is hexagonal or not).

[shyamd]

So at some point, we have to consider the fact that this is for the MPRelaxSet, which suggests it reflects MP as is.

[computron]

@shyuep FYI - this can change grids quite a bit

Si structure relaxation was 8x8x8 (old) and 6x6x6 (new). Note that the new scheme seems to more closely resemble the desired grid (6.02 -> 6 instead 6.02 -> 8)

Si uniform run was 220 explicit kpoints (old) and is now 145 explicit kpoints - this should have been the "per volume" scheme rather than "per atom". This is a huge reduction(!)

The new algorithm is probably more stable, but I don't know whether the tests of k-point convergence need to be re-run and some of the default k-point densities cranked up. e.g., if we were previously having 220 kpoints for the Si uniform run, I am not sure that 145 is going to cut it

[computron]

Also the NEB runs in atomate tests used to 8x8x8. Now they are 3x3x3.

These are pretty big differences; I'm not sure this new k-mesh is well tested enough for atomate. Previously I was going to update all the k-point tests in atomate for new pymatgen, but for now I am just going to roll back the req pymatgen to 4.7.5 in atomate until we sort this out.

[shyuep]

Which NEB test are you referring to?

Si is a semiconductor. I am pretty sure a 6x6x6 grid is more than enough, even for the primitive cell. The number of kpoints, even without taking into account symmetry, drops from 512 to 256. So I am not surprised the uniform unique kpts drop from 220 to 145.

[computron]

The NEB test in atomate:

atomate/vasp/workflows/tests/test_neb_workflow.py:40

which appears to be for Li2O

[shyuep]

Well, it is for a 2x2x2 of the primitive cell of Li2O, which puts it at 6.7 A in each direction. I highly doubt you need a 8x8x8 kpoint grid for that cell. A 3x3x3 grid sounds reasonable to me.

@johnson1228 @HanmeiTang Comment?

[computron]

Sorry - slight modification (some better, some worse)

For the NEB workflow, the endpoint calcs look to have been 4x4x4 in the past and 3x3x3 revised (=no big deal).

For the actual image calcs (i.e., multiple POSCAR, same INCAR/KPOINTS/etc) - the calcs look to be 8x8x8 in the past and 1x1x1 (!!) in the revised scheme. This is the major difference - i.e., not in the endpoints but rather the actual images.

[shyuep]

I see. I am pretty sure the image calculations are wrong. You cannot possibly be using a 8x8x8 for image calculations. That NEB calculation will never finish. The image calculations should be <= the grid of end points.

[computron]

I might be misinterpreting the NEB tests then. Looks like there are multiple "configurations" being tested and I can't make sense of it.

The reference files for the old kpoints are here: atomate/atomate/vasp/test_files/neb_wf

There are 5 directories in there. I thought some of these were endpoints and others were the images, but it looks likes these are 5 different configurations of NEB runs as documented in:

atomate/vasp/workflows/tests/test_neb_workflow.py

The first 3 of these configurations show a rough match between previous and current kpoints (4x4x4 vs 3x3x3). The next two (configs 4 + 5) have 8x8x8 in the KPOINTS file - which change to 1x1x1 in the new scheme.

[shyuep]

Ok. @HanmeiTang @johnson1228 Pls fix the tests in atomate PROPERLY.

[computron]

Note that for configs 4 and 5, the KPOINTS file itself documents that the previously the grid points per atom were 1000, whereas in the new file it says 100 (an order of magnitude less).

OLD

pymatgen generated KPOINTS with grid density = 1000 / atom 0 Monkhorst 8 8 8

NEW

pymatgen 4.7.6+ generated KPOINTS with grid density = 100 / atom 0 Gamma 1 1 1

[HanmeiTang]

Okay, we will fix the unittests soon.

[HanmeiTang]

We have fixed the unittests issue.

Basically, the tests that we wrote before only check if file transfer using the RunNEBVaspFake class works properly in the NEB workflow, i.e. actual NEB calculations are not performed. The different KPOINTS settings from 1000/atom to 100/atom are just to make sure that KPOINTS file can be modified by the users.

In the KPOINTS file for perfect cell relaxations as well as endpoint relaxations, we have changed the k-mesh to Gamma 3x3x3, and Gamma 1x1x1 for NEB calculations.

[computron]

Ok I am going to close this again for now. At some point I might try to do a test of "new vs old" k-point meshes. I am still a little nervous because we benchmarked the old scheme, and now it seems we are almost always reducing the k-point mesh from that.

[computron]

also @HanmeiTang please push a PR to atomate when you have the chance

[HanmeiTang]

@computron Will do.

[shyuep]

@computron I have already done tests for Li2O and NaCl kpoint meshes. Generally, large band gap insulators are not going to be an issue. These meshes are almost certainly too small for metals. It is more the intermediate semiconductor cases that might be of concern.

Certainly the new scheme will always reduce the mesh. The idea was that the rounding and even/odd rules of the old mesh algo always rounds up, and sometimes by a significant amount. I highly doubt there is a major effect for structural relaxations.

For static and uniform calculations, I would certainly recommend upping the grid by 1.5x or even 2x if you feel uncomfortable. Those are cheap single SCF calculations.

Quite frankly, we should not be worrying about kpoint meshes for relaxations when MP is still using EDIFF * # of atoms, which sometimes result in EDIFF > 0.005 eV, and correspondingly large EDIFFGs. The error of those washes out the small errors you are likely to get with the kpoints. Many structures in MP, especially those with a large # of atoms, are actually poorly converged in terms of structure.

[shyuep]

@computron Final update - I have cranked up the reciprocal density for relaxation to 64 instead of 50. This results in a 7x7x7 grid for Si.

Note that the VASP examples themselves show a grid of 11x11x11 for fcc Si. (https://cms.mpi.univie.ac.at/wiki/index.php/Fcc_Si) Recp density = 50 gives a 9x9x9 grid, which is not far off. Recp density = 64 gives 10x10x10, which I think is pretty reasonable for HT.

The correct behavior at high lattice lengths is maintained, i.e., after about 12 A, the gamma grid is obtained. I still think that 2x2x2 mesh for a 10Ax10Ax10A structure is a bit too much, but I will live with this.