BUG: Dependance on the k mesh of VASP

[Lattay]

Describe the bug Despite being able to load data from other codes through the Pymatgen BandStructure object, IFermi produces absurd results when provided with a K mesh different from VASP's.

In particular, I tried to do the same computation in VASP and Quantum-espresso, computing the Fermi surface of metal copper.

With VASP no problem, but with QE I get completely absurd result. I tracked down the problem to the fact that interpolation produce very high and very low energies in all bands, causing each band to intersect the Fermi level.

I am trying to understand what is going wrong in the interpolation, but I think it is related to the fact that QE and VASP uses different IBZ (see following picture).

To Reproduce

1. Band data here (don't extract them, they are npz files, but I renamed them for github to take them) qe.zip vasp.zip
2. Fermi surface script:

   
   #!/usr/bin/env python3
   import numpy as np
   import sys

from pymatgen.core import Lattice, Structure from pymatgen.electronic_structure.bandstructure import BandStructure, Spin

from ifermi.surface import FermiSurface from ifermi.interpolate import FourierInterpolator from ifermi.plot import FermiSurfacePlotter, show_plot from ifermi.kpoints import kpoints_from_bandstructure

def load_banstructure(path): data = dict(np.load(path, allow_pickle=True)) kpts = data["kpts"] lattice = Lattice(data["lattice"]) struct = Structure.from_dict(data["struct"][()])

eigs = data["eigs"]
efermi = data["efermi"]

return BandStructure(
    kpts,
    {Spin.up: eigs},
    lattice,
    efermi,
    structure=struct,
)

def plot_fermi_surface(bs: BandStructure): for i, band in enumerate(bs.bands[Spin.up]): print(np.min(band - bs.efermi), np.max(band - bs.efermi))

# Complete the BZ (source is just in the IBZ)
interpolator = FourierInterpolator(bs)
dense_bs, velocities = interpolator.interpolate_bands(
    return_velocities=True,
    interpolation_factor=2,
)

# generate the Fermi surface and calculate the group velocity at the
# center of each triangular face
dense_kpoints = kpoints_from_bandstructure(dense_bs)

fs = FermiSurface.from_band_structure(
    dense_bs,
    mu=0.0,
    wigner_seitz=True,
    calculate_dimensionality=True,
    property_data=velocities,
    property_kpoints=dense_kpoints,
)

# number of isosurfaces in the Fermi surface
print("n surf =", fs.n_surfaces)

# number of isosurfaces for each Spin channel
print("n surf/spin =", fs.n_surfaces_per_spin)

# the total area of the Fermi surface
print("area =", fs.area)

# plot the Fermi surface
fs_plotter = FermiSurfacePlotter(fs)
plot = fs_plotter.get_plot()
show_plot(plot)  # displays an interactive plot

if name == "main": if "-h" in sys.argv or "--help" in sys.argv or len(sys.argv) < 2: print("usage: fermi_surface.py ") sys.exit(0) bs = load_banstructure(sys.argv[1]) plot_fermi_surface(bs)

3. `python fermi_surface.py vasp.zip` behave as expected
4. `python fermi_surface.py qe.zip` uses huge amounts of ram, takes a really long time and produces hundreds of bogus surfaces

**Expected behavior**
There should be no dependance on the exact IBZ chosen.

[utf]

Hi @Lattay, thanks for this bug report. Can you confirm that your reciprocal k-points for QE are correct. I see that they span into the neighbouring Brillouin zone (0 to 0.67343506 in fractional coordinates). Is this definitely right?

[utf]

Actually, I just figured it out. The QE k-points are off by a factor of 1/sqrt(2).

If you update it to:

kpts = data["kpts"] / np.sqrt(2) 

That will fix your issue. Also, FYI, the lattice used in the BandStructure object is supposed to be the reciprocal lattice. It doesn't matter in this case since we're already in fractional coordinates, but a more robust function would look like:

def load_banstructure(path, cart_coords=True):
    data = dict(np.load(path, allow_pickle=True))
    kpts = data["kpts"] / np.sqrt(2)
    lattice = Lattice(data["lattice"])
    struct = Structure.from_dict(data["struct"][()])

    eigs = data["eigs"]
    efermi = data["efermi"]

    return BandStructure(
        kpts,
        {Spin.up: data["eigs"]},
        lattice.reciprocal_lattice,
        efermi,
        structure=struct,
        coords_are_cartesian=False
    )

[Lattay]

Actually, I just figured it out. The QE k-points are off by a factor of 1/sqrt(2).

Oh, that's anoying... thank you for catching that!

Edit: turns out the bug was in the QE reader from ASE: https://gitlab.com/ase/ase/-/merge_requests/3527

Also, FYI, the lattice used in the BandStructure object is supposed to be the reciprocal lattice.

That is what is stored in the (arguably poorly named) data["lattice"] field.