Hello everybody,
I am trying to do surface energy calculations and naturally want to set the kpoint meshes equivalent for the slab calculation and the bulk calculation, other than for the non-periodic slab axis, where the number of k-points should be set to 1.
Before I converge the kpoint density with respect to the total energy for a simple unit cell of the material using gamma-centered meshes, giving me a converged kpoint density kappa.
Now I construct a slab and an oriented unit cell (OUC) and use the automatic_gamma_density method of the pymatgen.io.vasp.inputs.Kpoints class again to construct k-point meshes with the converged density value. However, although the first two lattice vectors are equivalent, for the slab cell and the OUC, the resulting meshes are not.
After some investigation I found that this is rooted in the scaling of kappa with the number of sites from the input structure. This makes no sense to me and I wanted to ask if this is indeed intended behaviour and why this is set up like this. Of course I should have noticed that the scaling with respect to the number of atoms is mentioned in the comment line of the Kpoint object, but nevertheless I am confused by this choice.
In the docstring of the method the algorithm is described without that scaling:
Algorithm:
Uses a simple approach scaling the number of divisions along each
reciprocal lattice vector proportional to its length.
Here is a small example for a 001 slab of gold (disregarding that for the third axis of the slab the number of points should be set to 1 manually):
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 12 17:16:57 2021
@author: Michael Wolloch
"""
from pymatgen.core.structure import Structure
from pymatgen.core.surface import Slab
from pymatgen.core.surface import SlabGenerator
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.io.vasp import Kpoints
prim_bulk_dict = {'@module': 'pymatgen.core.structure',
'@class': 'Structure',
'charge': None,
'lattice': {'matrix': [[0.0, 2.079028483254921, 2.079028483254921],
[2.079028483254921, 0.0, 2.079028483254921],
[2.079028483254921, 2.079028483254921, 0.0]],
'a': 2.9401902775790747,
'b': 2.9401902775790747,
'c': 2.9401902775790747,
'alpha': 60.00000000000001,
'beta': 60.00000000000001,
'gamma': 60.00000000000001,
'volume': 17.972616757073553},
'sites': [{'species': [{'element': 'Au', 'occu': 1}],
'abc': [0.0, 0.0, 0.0],
'xyz': [0.0, 0.0, 0.0],
'label': 'Au',
'properties': {'magmom': -0.0}}]}
kappa = 1250
prim_bulk = Structure.from_dict(prim_bulk_dict)
bulk_conv = SpacegroupAnalyzer(prim_bulk).get_conventional_standard_structure()
SG = SlabGenerator(initial_structure=bulk_conv,
miller_index=[0,0,1],
center_slab=True,
primitive=True,
lll_reduce=True,
max_normal_search=1,
min_slab_size=20,
min_vacuum_size=20)
slab = SG.get_slabs(bonds=None, ftol=0.1, tol=0.1, max_broken_bonds=0,
symmetrize=False, repair=False)[0]
ouc = slab.oriented_unit_cell
kpoints_slab = Kpoints.automatic_gamma_density(slab, kappa)
kpoints_ouc = Kpoints.automatic_gamma_density(ouc, kappa)
print('slab lattice abc:\n{}'.format(slab.lattice.abc))
print('ouc lattice abc:\n{}'.format(ouc.lattice.abc))
print()
print('kpoint mesh for the slab:\n{}'.format(kpoints_slab))
print('kpoint mesh for the ouc:\n{}'.format(kpoints_ouc))
Which results in:
slab lattice abc:
(2.940190277579075, 2.9401902775790743, 41.58056966509842)
ouc lattice abc:
(2.940190277579075, 2.9401902775790743, 4.158056966509842)
kpoint mesh for the slab:
pymatgen v2020.12.31 with grid density = 1250 / number of atoms
0
Gamma
13 13 2
kpoint mesh for the ouc:
pymatgen v2020.12.31 with grid density = 1250 / number of atoms
0
Gamma
11 11 8
I thought about just scaling kappa with the number of sites as well:
kpoints_slab = Kpoints.automatic_gamma_density(slab, kappa*slab.num_sites)
kpoints_ouc = Kpoints.automatic_gamma_density(ouc, kappa*slab.num_sites)
print()
print('adjusted kpoint mesh for the slab:\n{}'.format(kpoints_slab))
print('adjusted kpoint mesh for the ouc:\n{}'.format(kpoints_ouc))
but this leads to even more different meshes:
adjusted kpoint mesh for the slab:
pymatgen v2020.12.31 with grid density = 12500 / number of atoms
0
Gamma
27 27 2
adjusted kpoint mesh for the ouc:
pymatgen v2020.12.31 with grid density = 2500 / number of atoms
0
Gamma
13 13 9
although now the set density should be 1250 for real for both systems since the slab has 10 sites and the OUC has 2.
I would appreciate some help or explanation for this issue, since it does not seem to be very practical.
Thanks and all the best, Michael