How to Calculate Grid Density for Given KSPACING?

[hszhao.cn@gmail.com]

Hi here,

I'm currently working on setting up my VASP calculations and trying to understand how to derive the k-point grid mess size based on the KSPACING parameter. Could you guide me on how to calculate or estimate the grid density from a given KSPACING value?

Say, I set the following in INCAR:

Code: Select all

KSPACING = 0.02  in \AA ^-1

Then, I want to know the corresponding n1 * n2 * n3 mesh when using a KPOINTS file.

Any advice or pointers to relevant documentation would be greatly appreciated.

See here for the related discussion.

Thanks,
Zhao

[martin.schlipf]

If you want to compute the grid by hand, I recommend studying the description of the INCAR tag KSPACING. There it is explained how the grid density is selected based on the value that is set and whether KGAMMA is set or not.

Many times it may be easier to just run VASP to look at the generated mesh though. You may want to use the dry-run mode to skip over the expensive parts of the calculation

Code: Select all

mpirun -np 1 /path/to/vasp_std --dry-run

Then VASP will finish after the k-point mesh has been generated and you will find an output like

Code: Select all

Automatic generation of k-mesh.
 Grid dimensions derived from KSPACING:
 generate k-points for:   10   10   10

in the OUTCAR file that contains the information you need.

[hszhao.cn@gmail.com]

Dear martin.schlipf,

Based on the specification described in INCAR tag KSPACING and source code file mkpoints.F, in this case, NKPX, NKPY, and NKPZ is derived using the lattice norms (BNORM) and the spacing (KSPACING). The optional multiplier (MULTIPL) is also considered if provided. Here's the implementation in Python:

Code: Select all

import math

def calculate_grid_dimensions(bnorm, spacing, multipl=None):
    """
    Calculate the Monkhorst-Pack grid dimensions based on lattice norms and spacing.

    :param bnorm: Lattice norms (list or tuple of three floats).
    :param spacing: KPOINTS spacing (float).
    :param multipl: Optional multiplier (list or tuple of three integers).
    :return: Tuple of grid dimensions (NKPX, NKPY, NKPZ).
    """
    # Calculate initial grid dimensions based on spacing
    nkpx = max(1, math.ceil(bnorm[0] * math.pi * 2 / spacing))
    nkpy = max(1, math.ceil(bnorm[1] * math.pi * 2 / spacing))
    nkpz = max(1, math.ceil(bnorm[2] * math.pi * 2 / spacing))

    # Apply multiplier if provided
    if multipl:
        nkpx *= multipl[0]
        nkpy *= multipl[1]
        nkpz *= multipl[2]

    return nkpx, nkpy, nkpz

# Example usage
bnorm = [1.0, 2.0, 3.0]  # Example lattice norms
spacing = 0.5  # Example KPOINTS spacing
multipl = [2, 2, 2]  # Example multiplier

nkpx, nkpy, nkpz = calculate_grid_dimensions(bnorm, spacing, multipl)
print("Grid dimensions:", nkpx, nkpy, nkpz)

This Python function calculate_grid_dimensions mirrors the Fortran logic for determining the Monkhorst-Pack grid dimensions (NKPX, NKPY, NKPZ) based on the lattice norms (BNORM), the Monkhorst-Pack spacing (KSPACING), and an optional multiplier (MULTIPL).

However, one of my doubts is that the parameter MULTIPL does not seem to be mentioned in the input files of vasp, so I do not know how to provide it to the above function for use. If there are any other logical problems, welcome further communication and correction.

Regards,
Zhao

[martin.schlipf]

I think the only point where the multiplier is used is for cRPA calculations or do you have a calculation which runs multiple times through this routine? In that particular case the multiplier seems to be hard coded looping over 1..36.

[hszhao.cn@gmail.com]

Dear martin.schlipf,

I don't seem to have encountered any such situation before.

Regards,
Zhao

[hszhao.cn@gmail.com]

Dear martin.schlipf,

Many times it may be easier to just run VASP to look at the generated mesh though. You may want to use the dry-run mode to skip over the expensive parts of the calculation

Code: Select all

mpirun -np 1 /path/to/vasp_std --dry-run

Then VASP will finish after the k-point mesh has been generated and you will find an output like

Code: Select all

Automatic generation of k-mesh.
 Grid dimensions derived from KSPACING:
 generate k-points for:   10   10   10

in the OUTCAR file that contains the information you need.

Using vasp.6.4.3, I tried your above method, but the result information is not exactly the same as you showed:

Code: Select all

werner@x13dai-t:~/Desktop/Jinyu-Qiao/Bulk_Systems_-_Tutorial/Fcc_Si$ ug -i kspa
vasprun.xml:    <i name="KSPACING">      0.50000000</i>
werner@x13dai-t:~/Desktop/Jinyu-Qiao/Bulk_Systems_-_Tutorial/Fcc_Si$ ug -A2 -B3 -i 'grid dim'
OUTCAR- KPOINTS: pymatgen with grid density = 787 / numbe
OUTCAR-
OUTCAR-Automatic generation of k-mesh.
OUTCAR: Grid dimensions read from file:
OUTCAR- generate k-points for:    7    7    7
OUTCAR-

More specifically, as you can see, in my case, there are two points that need to be clarified:

1. KSPACING info only exists in vasprun.xml.
2. In OUTCAR, I see the following line:

Code: Select all

Grid dimensions read from file:

But you showed the following line:

Code: Select all

Grid dimensions derived from KSPACING:

Therefore, it seems that the storage position of the output results of these related information has been adjusted in the latest version of VASP.

Regards,
Zhao

[hszhao.cn@gmail.com]

In addition, I have another questione that are troubling me regarding the results I have obtained:

In the above test, I generate the input files using pymatgen, as a result, there is the following relevant information appears in KPOINTS and OUTCAR:

Code: Select all

(datasci) werner@x13dai-t:~/Desktop/Jinyu-Qiao/Fcc_Si$ ug -i pymatgen
KPOINTS:pymatgen with grid density = 787 / number of atoms
Binary file vaspout.h5 matches
OUTCAR: KPOINTS: pymatgen with grid density = 787 / numbe
OUTCAR: k-points in units of 2pi/SCALE and weight: pymatgen with grid density = 787 / numbe
OUTCAR: k-points in reciprocal lattice and weights: pymatgen with grid density = 787 / number

(datasci) werner@x13dai-t:~/Desktop/Jinyu-Qiao/Fcc_Si$ cat KPOINTS 
pymatgen with grid density = 787 / number of atoms
0
Gamma
7 7 7

But I wonder what's the relationship between pymatgen grid density and VASP KSPACING?

Regards,
Zhao

[martin.schlipf]

Can you provide the input files you used? Looking just at the output it appears that you still provided a KPOINTS file. KSPACING is ignored if you use a KPOINTS file.

I do not know how pymatgen sets up there k-point mesh so I cannot comment on the differences to VASP.

[hszhao.cn@gmail.com]

Dear martin.schlipf,

Looking just at the output it appears that you still provided a KPOINTS file. KSPACING is ignored if you use a KPOINTS file.

Thank you for pointing this out. Now, I use the following code to generate input files and confirm all your previous assertions:

Code: Select all

from mp_api.client import MPRester
from pymatgen.io.vasp.inputs import Kpoints
from pymatgen.io.vasp.sets import MPRelaxSet

# Use Materials Project API to get Si structure with Fd-3m space group
with MPRester() as mpr:
    results = mpr.materials.summary.search(
        formula="Si",
        spacegroup_symbol="Fd-3m"
    )

    if results:
        result = results[0]
        mp_id = result.material_id
        structure = mpr.get_structure_by_material_id(mp_id)
    else:
        print("No matching materials found.")
        exit()

# Set k-point density
kppa = 787

# Generate k-points using pymatgen's method
kpoints = Kpoints.automatic_density(structure, kppa)

# Get the actual k-point grid
grid = kpoints.kpts[0]

# User-defined INCAR settings
user_incar_settings = {
    'ALGO': 'N',
    'GGA': 're',
    'KSPACING': 0.3,
#    'KGAMMA': '.TRUE.',
}

# Create VASP input set
vasp_input = MPRelaxSet(structure, user_potcar_functional="PBE_64", user_incar_settings=user_incar_settings)
vasp_input.write_input("./Fcc_Si")

# Print results
print(f"Material ID: {mp_id}")
print(f"kppa (grid density) = {kppa}")
print(f"Generated k-point grid: {grid[0]} x {grid[1]} x {grid[2]}")

The run results of the above python code are as follows:

Code: Select all

Material ID: mp-149
kppa (grid density) = 787
Generated k-point grid: 7 x 7 x 7

The vasp dry-run of the input files generated above will give the following results:

Code: Select all

$ ug -A2 -B2 -i 'grid dim' -g OUTCAR 
OUTCAR-
OUTCAR-Automatic generation of k-mesh.
OUTCAR: Grid dimensions derived from KSPACING:
OUTCAR- generate k-points for:    7    7    7
OUTCAR-

I do not know how pymatgen sets up there k-point mesh so I cannot comment on the differences to VASP.

Based on the test, in this case, a 0.3 KSPACING will generate the same grid dimension as pymatgen kppa 787.

Please find all input and output files for the VASP dry-run in the attachment.

Regards,
Zhao

[hszhao.cn@gmail.com]

Dear martin.schlipf,

The source code of pymatgen automatic_density is located here, as shown below:

Code: Select all

    def automatic_density(cls, structure: Structure, kppa: float, force_gamma: bool = False) -> Self:
        """Get an automatic Kpoint object based on a structure and a kpoint
        density. Uses Gamma centered meshes for hexagonal cells and face-centered cells,
        Monkhorst-Pack grids otherwise.

        Algorithm:
            Uses a simple approach scaling the number of divisions along each
            reciprocal lattice vector proportional to its length.

        Args:
            structure (Structure): Input structure
            kppa (float): Grid density
            force_gamma (bool): Force a gamma centered mesh (default is to
                use gamma only for hexagonal cells or odd meshes)

        Returns:
            Kpoints
        """
        comment = f"pymatgen with grid density = {kppa:.0f} / number of atoms"
        if math.fabs((math.floor(kppa ** (1 / 3) + 0.5)) ** 3 - kppa) < 1:
            kppa += kppa * 0.01
        lattice = structure.lattice
        lengths: Vector3D = lattice.abc
        ngrid = kppa / len(structure)
        mult: float = (ngrid * lengths[0] * lengths[1] * lengths[2]) ** (1 / 3)

        num_div: Tuple3Ints = cast(Tuple3Ints, [math.floor(max(mult / length, 1)) for length in lengths])

        is_hexagonal: bool = lattice.is_hexagonal()
        is_face_centered: bool = structure.get_space_group_info()[0][0] == "F"
        has_odd: bool = any(idx % 2 == 1 for idx in num_div)
        if has_odd or is_hexagonal or is_face_centered or force_gamma:
            style = cls.supported_modes.Gamma
        else:
            style = cls.supported_modes.Monkhorst

        return cls(
            comment,
            0,
            style,
            [num_div],
            (0, 0, 0),
        )

But it seems that analyzing the relationship between VASP KSPACING and Pymatgen's kppa (k-point per atom) from this is still a very difficult task.

Regards,
Zhao

[martin.schlipf]

I had a look at the pymatgen code and I don't think you will be able to write a code that translates kppa into KSPACING.

The VASP approach is very simple, you take the norm of the reciprocal lattice vector and multiply it by 2π/KSPACING then round it to the next larger integer. This has the problem, that you can get unlucky where with small changes to the cell your rounding jumps by 1 for the k-point grid.

In pymatgen, they try to correct this by dividing the cube root of the volume by the length of a vector. This should be constant for simple scaling of the cell. The expected coarser k-point mesh for larger systems comes from a division by the number of atoms instead. However, this is no way around the quantization problem. It eliminates the volume one but you could still change c/a for a cell at constant volume and run into the jump of 1 for the grid. By dividing by the number of atoms, they introduced a new problem that two structures in the same cell may have a different k-point mesh if one of them has a defect.

I would not recommend use either of these approaches for production calculations where the volume or shape of the cell is allowed to change. In all of these scenarios, I would run this setup once and then copy the generated IBZKPT file to KPOINTS and use that for all further calculations.