high-symmetry k path for arbitrary lattice

[jaemolihm]

Summary

I have implemented a function that gives the high-symmetry k path for an arbitrary 3D lattice, where the k path is equivalent to the one that one would get for the standardized primitive lattice. I would like to know if Brillouin.jl developers are interested to have this functionality. If so, I'll make a PR.

Motivation

Currently Brillouin.jl (and also SeeK-path) computes only "the suggested band paths only of standardized (crystallographic) primitive cells." So, "If you already have done calculations with a non-standardized cell, you will then need to figure out how to remap the labeled k-points in the choice of cell you did." (https://seekpath.readthedocs.io/en/latest/maindoc.html#a-warning-on-how-to-use-and-crystal-structure-standardization) This remapping problem is what I want to address here.

Method

The input to irrfbz_path are the space group number sgnum and the conventional lattice vectors Rs. The returned k-points are given in the basis of the primitive reciprocal basis in the CDML setting. Let's say the CDML lattice vectors Rs_p. So, the problem happens when one has a model whose lattice vectors Rs_model are not equal to Rs_p.

There are two cases.

Case 1: Rs_model describes the same system with Rs_p, but just in a different basis.

In other words, Rs_model[i] = sum_{j=1}^{3} N[i, j] * Rs_p[j] for some 3*3 integer matrix N. In this case, the same k path in Cartesian coordinates can be used for Rs_p and Rs_model. Hence, one can simply transform the k path as follows (could be shorter if #12 is resolved):

kp_cart = cartesianize(kp)
recip_basis = Bravais.ReciprocalBasis(SVector{3}(eachcol(inv(hcat(lattice...))' .* 2π)))
for key in keys(kp.points)
    kp.points[key] = Brillouin.latticize(kp_cart.points[key], recip_basis)
end
# Need to create a new object because kp.basis is not mutable
KPath(kp.points, kp.paths, recip_basis, kp.setting)

Case 2: Rs_model describes a different system from Rs_p.

In this case, I use the following assumption:

Assumption: The system described Rs_model can be obtained by applying a symmetry operation S to the system described by Rs_p. Here, S is a symmetry of the conventional lattice system Rs, but not a symmetry of the primitive lattice Rs_p.

Unfortunately, I do not have a mathematical proof of this assumption. But let me assume this is true.

Then, we can loop over the symmetry operations of Rs and see if it converts Rs_p to Rs_model. I used Spglib.get_symmetry to get the symmetries, and looped over the symmetries.

Example

As an example, I plot the k path for space group 166 (material example: Bi2Se3 https://materialsproject.org/materials/mp-541837/) For the standard cell, the old and new kpath are identical. For the non-standard cell (rotated by 60 degree along the z axis), the L point (marked in red) is at the correct position only in the new kpath.

Code

using StaticArrays
using Brillouin
using PlotlyJS
using LinearAlgebra
import Spglib, Bravais

function irrfbz_path_new(sgnum, lattice, conv_lattice)
    # sgnum: space group number.
    # lattice: primitive lattice vectors. May not be the standard one.
    # conv_lattice: conventional lattice vectors.

    # Convert to 3*3 Matrix to do linear algebra
    lattice_mat = hcat(lattice...)

    # Get the standard primitive lattice
    prim_lattice = Bravais.primitivize(Bravais.DirectBasis(collect(conv_lattice)), Bravais.centering(sgnum, 3))

    # Check whether lattice describes the same physical system with prim_lattice
    lattice_transf = hcat((inv(hcat(prim_lattice...)) * lattice_mat)...)
    if round.(Int, lattice_transf) ≈ lattice_transf
        # In this case, we do not need any additional rotation of the k path.
        sym_cart = Mat3(I(3))
    else
        # The physical systems described by lattice and prim_lattice are different.
        # Need to rotate the k path in Cartesian space using a symmetry that maps prim_lattice to lattice.
        # Look for such a symmetry among the symmetries of the conventional cell.
        conventional_cell = Spglib.Cell(conv_lattice, [0, 0, 0], [0])
        symmetry_conventional, _ = Spglib.get_symmetry(conventional_cell)
        conv_lattice_mat = hcat(conv_lattice...)

        isym_found = 0
        for isym in 1:size(symmetry_conventional, 3)
            s_cart = conv_lattice_mat * symmetry_conventional[:,:,isym]' * inv(conv_lattice_mat)
            # rotate prim_lattice
            rotated_prim_lattice = s_cart * hcat(prim_lattice...)
            # Check wheter lattice = rotated_prim_lattice * lattice_transf where lattice_transf is an integer matrix
            lattice_transf = hcat((inv(rotated_prim_lattice) * lattice_mat)...)
            if round.(Int, lattice_transf) ≈ lattice_transf
                isym_found = isym
                break
            end
        end
        isym_found == 0 && error("Symmetry mapping prim_lattice to lattice not found.")
        sym_cart = Mat3(conv_lattice_mat * symmetry_conventional[:, :, isym_found]' * inv(conv_lattice_mat))
    end
    @show sym_cart

    # Get the k path (in reduced basis for the reciprocal of prim_lattice)
    kp = irrfbz_path(sgnum, conv_lattice)

    # Now, kp is in crystal coordinates of the primitive lattice vector in CDML convention (prim_lattice),
    # not the input lattice. So, we convert kp to crystal coordiantes in lattice by
    # converting as follows:
    # crystal (prim_lattice) -> Cartesian -> rotate using sym_cart -> crystal (input lattice)
    kp_cart = cartesianize(kp)
    recip_basis = Bravais.ReciprocalBasis(SVector{3}(eachcol(inv(hcat(lattice...))' .* 2π)))
    for key in keys(kp.points)
        k_cart_rotated = sym_cart * kp_cart.points[key]
        kp.points[key] = Brillouin.latticize(k_cart_rotated, recip_basis)
    end
    # Need to create a new object because kp.basis is not mutable
    KPath(kp.points, kp.paths, recip_basis, kp.setting)
end

a = 1.0
c = 8.0
# Trigonal lattice, space group R-3m, 166
# Example: Bi2Se3 (https://materialsproject.org/materials/mp-541837/)
lattice_nonstandard = SVector(SVector(a*sqrt(3)/2, -a/2, c/3),
                              SVector(0.0, a, c/3),
                              SVector(-a*sqrt(3)/2, -a/2, c/3))
conv_lattice = SVector(SVector(a*sqrt(3), 0, 0),
                       SVector(-a*sqrt(3)/2, a*3/2, 0),
                       SVector(0, 0, c))
sgnum = 166
lattice_standard = SVector(Bravais.primitivize(Bravais.DirectBasis(collect(conv_lattice)), Bravais.centering(sgnum, 3)))

kp_old = irrfbz_path(sgnum, conv_lattice)
kp_standard_new = irrfbz_path_new(sgnum, lattice_standard, conv_lattice)
kp_nonstandard_new = irrfbz_path_new(sgnum, lattice_nonstandard, conv_lattice)

ws_standard = wignerseitz(Bravais.reciprocalbasis(lattice_standard))
ws_nonstandard = wignerseitz(Bravais.reciprocalbasis(lattice_nonstandard))

p1 = PlotlyJS.plot(ws_standard, kp_old, Layout(title="standard cell, old kpath"))
p2 = PlotlyJS.plot(ws_standard, kp_standard_new, Layout(title="standard cell, new kpath"))
p3 = PlotlyJS.plot(ws_nonstandard, kp_old, Layout(title="non-standard cell, old kpath"))
p4 = PlotlyJS.plot(ws_nonstandard, kp_nonstandard_new, Layout(title="non-standard cell, new kpath"))

[mfherbst]

Cool idea. Would be nice to have something like this.

I wonder if it would make more sense to avoid the dependency of Brillouin on Spglib and rather move the whole symmetry / Brillouin path / band structure functionality to a separate package, that logically sits on top of both. We have quite a bit such code in DFTK and we would would not be too unhappy to move this out of our code base (as it does not really belong to our "core" functionality of DFT). I guess with the emerging AtomsBase.jl package it should not be too difficult to hook everything up in a nice way.

[thchr]

This is very cool and would be really nice to have some form of this in Brillouin: expecting every use-case to necessarily be in the form of a standard lattice is definitely very sub-optimal. Thanks for suggesting and working on this @jaemolihm.

Here's a few thoughts:

• As @mfherbst notes, it would be really nice if we could manage to do this without depending on Spglib. What is the functionality that is needed from Spglib? I didn't completely understand the calls to Spglib in your code example. I have a separate project, Crystalline.jl, that has functionality for space group theory (operations) etc. - but that is also a bit of a heavy dependency because it includes a bunch of other things (irreps, band representations, etc.).

• Regarding Case 1: I'm actually even a bit unsure if this case is trivial. If the desired input lattice basis (let's call this Rs′, assumed to be in a non-primitive, i.e. conventional, but not standard setting) is rotated relative to the standard conventional lattice basis Rs, then what will you input to irrfbz_path to create kp in your example? If we simply give it Rs′, things will probably not work out well (because there's a branch table that assumes a conventional input lattice basis, see e.g. the case of Bravais type oI). Basically, how will you create the initial kp variable in your example if Rs_model (this is what I called Rs′ I believe) is not in a standard setting? The rest of your code makes sense to me.

• Regarding Case 2: The assumption you mention is probably not true, I believe. Case in point: applying any operation from O(3) (i.e., all rotations and inversions/mirrors) will in general produce a new lattice and a new basis - but the described system is of course fundamentally still the same (ignoring chirality-related considerations). So, even if Rs_model and Rs_p describe the same primitive cell, they could differ in general by any orthogonal transformation. I think one general way to map to a standard basis could be to look at the angles and length-ratios between lattice vectors, and then try to obtain a transformation based on needing to map these angles/ratios to the standard setting (by comparing e.g., with this table from Bravais.jl).

• Case 1 versus 2: Having written all this, I am actually not completely sure I understand the distinctions you make between these two cases (sorry!) - could you explain it in a different way?

• Case 3: I wanted to mention another possibility: the user picked a lattice basis which cannot be rotated/reflected into a standard setting. There's nothing we can do then, as far as I can see. (I.e., the basis choice is not unique: any basis matrix B can be mapped to a new, equivalent lattice basis matrix B′ = BT for unimodular integer matrix T). An example could be a 2D hexagonal lattice: the standard basis has an angle of 120° between basis vectors; but 60° would be an equally good choice to describe the same lattice. In general, this seems too hard to deal with.

(A tangentially related issue with the input to irrfbz_path is that it is probably not very ergonomic to need to supply the conventional lattice vectors to get a path referred to the primitive reciprocal lattice: I suspect that most users will usually have the primitive lattice vectors rather than the conventional ones - so we should probably also make it possible to allow irrfbz_path to take primitive lattice vectors instead (e.g., toggled via a keyword argument primitive).)

[thchr]

I think a good interface to solve this would be to have a function, e.g., in Bravais.jl, that takes a conventional (or primitive) direct lattice and a crystal system type (or, if primitive, a Bravais type) and maps it to an equivalent standard lattice basis, returning the standard lattice basis and the associated transformation between the two.

Then we could (1) call that function in irrfbz_path on the input Rs′ to transform to a standard lattice basis Rs; (2) construct the k-path associated with Rs, and (3) apply the inverse transformation to this k-path so that it again refers to the basis Rs (or, really, its primitive reciprocal equivalent).

[jaemolihm]

Thank you both for the comments!

The call to Spglib is to find the point group symmetries of the conventional unit cell. If there is another lightweight package that can do this, one can use it. I like the idea of a new package, but I am not sure about what should be in Brillouin.jl and what in the new package.

@thchr I think I did not understand your comments on cases 1 and 2. So let me explain my ideas in a different way.

First, I think we should consider only the case where two lattices has the same space group number and conventional lattice. So, in your comment about Case 2, a generic O(3) operation (e.g. rotation by 1 degrees) will change the conventional lattice, and this is not my concern. (By the same phsical system, I mean exactly the same lattice, not just up to some rotation.)

EDIT: I was misunderstanding the notion of standard conventional unit cell. You are right that my assumption is false. Thus, my current method does not work in general and we need a more general way like what you mentioned above.

Also, what you said as Case 3 is exactly what I deal with in Case 1.

so we should probably also make it possible to allow irrfbz_path to take primitive lattice vectors instead

Yes, that would be very useful. But the problem is that to use the SeeK-path algorithm, one needs a mapping from the primitive lattice to the standard conventional lattice. In both seek path and DFTK, Spglib is used for this. So one would need the Spglib dependency again.

I think a good interface to solve this would be to have a function, e.g., in Bravais.jl, that takes a conventional (or primitive) direct lattice and a crystal system type (or, if primitive, a Bravais type) and maps it to an equivalent standard lattice basis, returning the standard lattice basis and the associated transformation between the two.

Yes, such a mapping for the primitive case is similar to what I did in the function. It would be nice to make it a generic interface.

[jaemolihm]

It seems there is a similar need (k path for non-standard primitive cell) and some ideas from the seek-path and aiida developers. I'll take a look on them and see if there is something that we can use.

https://github.com/aiidateam/aiida-core/issues/4391

[jaemolihm]

@thchr , I have one question on the relation of Brillouin.jl and seek-path. Are they different (independent) implementations of the same algorithm? Or do you call seek-path here?

[thchr]

@thchr , I have one question on the relation of Brillouin.jl and seek-path. Are they different (independent) implementations of the same algorithm? Or do you call seek-path here?

We use the same paths as in SeeK-path, obtained by parsing their data-files (here; results stored here). This data is then used to build a number of parametrized functions for the possible paths over here. Finally, there's branching to the different possible paths.

So the short answer is that it's the same algorithm but we're not using the SeeK-path library.

[jaemolihm]

I see, thanks! You are right that my function does not work in general. I will spend some time and see how the problem can be solved.

[jaemolihm]

I found a general solution, though it makes a heavy use of spglib.

The essential ingredient is dset.std_rotation_matrix (dset is a spglib dataset), which is the rotation matrix needed to rotate the input system to the standard conventional/primitive lattice. The unimodular integer matrix part (i.e. linear combination of the basis vectors, Case 3 in https://github.com/thchr/Brillouin.jl/issues/13#issuecomment-967433553), can also be calculated, but does not need to be explicitly used. Just cartesianize and latticize does the work.

The only limitation is that one cannot get the ideal k path if the input cell is an undistorted supercell of some smaller cell. I think there is no easy, automatic way to handle this case.

One tricky part remaining is that SeeK-path uses a different "standard" conventional cell than spglib for the triclinic cell. Then, std_lattice might further need to be updated before being passed to Brillouin.irrfbz_path.

function irrfbz_path_for_cell(cell::Spglib.Cell)
    # standardize cell
    dset = Spglib.get_dataset(cell)
    sgnum = dset.spacegroup_number
    std_lattice = Bravais.DirectBasis(collect(eachcol(dset.std_lattice)))

    # If the input cell is a supercell (without any distortion), then the irrfbz algorithm cannot work
    if round(Int, det(cell.lattice) / det(dset.primitive_lattice)) != 1
        @warn "input cell is a supercell. irrfbz Does not give a correct k path."
    end

    # Calculate kpath for standard primitive lattice
    kp = Brillouin.irrfbz_path(sgnum, std_lattice)

    # Now, we convert from the standard primitive lattice to the original lattice.
    # The conversion formula is `cell.lattice = rotation * dset.primitive_lattice * transformation`.
    # `transformation` is not used later.
    rotation = dset.std_rotation_matrix
    transformation = inv(Bravais.primitivebasismatrix(Bravais.centering(sgnum, 3))) * dset.transformation_matrix'

    # Rotate k points in Cartesian space by `rotation`
    recip_basis = Bravais.reciprocalbasis(Bravais.DirectBasis(collect(eachcol(Matrix(cell.lattice)))))
    kp_cart = cartesianize(kp)
    for (lab, kv) in points(kp_cart)
        points(kp_cart)[lab] = rotation * kv
    end
    kp_new = latticize(kp_cart, recip_basis)

    kp_new, kp
end

[thchr]

That implementation is great - we should definitely add some form of this, I think.

I have a few questions:

• @mfherbst: what's the overlap of this implementation and what's currently done in DFTK.jl (mainly with functionality here). From what I understand, the intended work-flow now is that a given structure is first mapped to standard conventional cell using Spglib, and then a path is constructed for that cell using Brillouin? Presumably, @jaemolihm's would make it closer to a "one-call" work-flow?
• @jaemolihm: I don't understand the check for supercell structures. From a quick reading, it seems you are checking whether det(cell.lattice) ≈ det(dset.primitive_lattice), i.e. if the input (cell) and Spglib primitive cell (dset.primitive_lattice) unit cell volumes are identical? I don't understand why this implies a supercell (and wouldn't this break if a conventional cell is supplied? I'd be OK with that though.)
• This could probably live in e.g. /src/spglib.jl. I'm not sure if this should be conditionally available only if Spglib is loaded (i.e., with Requires) or be available by default. As context, loading Spglib.jl takes about ~0.8 s on my laptop, while loading Brillouin.jl is about 3 s currently (most of which comes from PyCall.jl, which I'm keen to drop when there's a stable Julia wrapper for QHull).

[thchr]

As a side note, I mostly wrote Brillouin for research I do in photonic crystals (where we don't have a simple equivalent of "atom positions") so in my ideal world, it would be possible to do some form of this without necessarily specifying an atomic crystal as input - but I realize that Spglib is very nice to leverage here, and that most people will care more about real crystals.

Is there some way to leverage Spglib to make the "plain" irrfbz_path(sgnum::Integer, Rs::DirectLattice) work even when Rs is a non-standard lattice?

[jaemolihm]

@jaemolihm: I don't understand the check for supercell structures. From a quick reading, it seems you are checking whether det(cell.lattice) ≈ det(dset.primitive_lattice), i.e. if the input (cell) and Spglib primitive cell (dset.primitive_lattice) unit cell volumes are identical? I don't understand why this implies a supercell (and wouldn't this break if a conventional cell is supplied? I'd be OK with that though.)

Yes, it checks that the volums are identical. The logic is that if the input cell is a supercell (N identical copies of the primitive lattice), then its volume will be N times that of the primitive lattice.

Also you are correct that the my irrfbz_path_for_cell breaks if a conventional cell is supplied. (For example, it will break if the cubic conventional cell of an fcc system is supplied.)

Is there some way to leverage Spglib to make the "plain" irrfbz_path(sgnum::Integer, Rs::DirectLattice) work even when Rs is a non-standard lattice?

I agree that this is desirable to do. But I could not find a solution. Spglib is designed for solid state systems, so I'm not sure it will be possible.

[thchr]

Okay, makes sense. If you have the cycles, a PR would be very welcome👍!

[mfherbst]

First off, I am very sorry for having taken eternities to react. I had a very heavy past few weeks, but I am glad to see how this discussion has evolved. Thanks for the effort @thchr and @jaemolihm !

@mfherbst: what's the overlap of this implementation and what's currently done in DFTK.jl (mainly with functionality here). From what I understand, the intended work-flow now is that a given structure is first mapped to standard conventional cell using Spglib, and then a path is constructed for that cell using Brillouin? Presumably, @jaemolihm's would make it closer to a "one-call" work-flow?

Essentially it covers lines 33 to 48, but it is more clever and does not bail out so easily if the given lattice is not primitive, so definitely a superior approach and I am all up for using this in DFTK.

Regarding implementation I would vouch for a conditional loading. I would also suggest a light wrapper around spglib in Brillouin.jl by using a named tuple to collect the relevant info of the spglib dataset (and convert it to proper Julia types). That tuple could than be passed onto a irrfbz_path_for_cell function (e.g. as kwargs?). That way Brillouin is fully functional without spglib (just pass on the data determined during the symmetry reduction manually) and does not depend in its interface on C types (which in my opinion is very bad). Also if someone makes a pure-Julia spglib in the future we can easily swap that in. A convenience interface based on spglib.cell could than be optionally supplied with Revise in case spglib is available.

Thoughts?

[jaemolihm]

Hi @mfherbst , thanks a lot for the comments! I made a PR just before I saw your comments 😃

Regarding your comments on passing information on spglib dataset, I think it is a good direction, but it is not included in my PR. My current PR works as is, so the spglib dataset functionality could be added either in the PR or as another PR.