thchr / brillouin.jl Goto Github PK

Currently the __init__ code in

From skimming the code I noted that PySpatial is only used inside the wignerseitz function, which is not used a lot inside Brillouin. Would it be ok to move the python import into that function (i.e. remove it from module scope) and make it optional such that only this function fails in case scipy is not installed? What I have in mind, is a pattern like

PySpatial = pyimport_e("scipy.spatial")
if ispynull(PySpatial)
   PyCall.conda || error("Install scipy")
   PySpatial = pyimport_conda("scipy.spatial", "scipy")
end

That way, for Conda.jl users everything is still automatic, otherwise you get an error, but only if the function is actually needed. Would that be ok with you @thchr ?

Hi,
the docstring mentions latticize!(kp::KPath{D}, Gs::BasisLike{D}), but the corresponding function was deleted in f32a826.

If I understood correctly, latticize!(kp, Gs) allows one to convert the Cartesian k path to reduced coordinates using any reciprocal basis Gs, not just using the default kp.basis. This functionality is very useful, so I hope the function is implemented again.

We don't need to create a separate .jll or anything like that, but if we want to keep parsing SeeK as a separate step (maybe this is just a terrible idea) we should not pollute the package directory during build but use the artifact system instead.

The paths for plane group 2 (and, similarly, other space groups with mP Bravais type; e.g., space group 3), currently only go from the origin to [0, ½], [½, 0], and [½, ½]. This is inherited from seekpath's choices but doesn't make a ton of sense: e.g., the lines are not actually symmetry lines in plane group 2 (since there aren't any) and, practically, one would likely like to include parts of the BZ boundary as well.

There are extra points in seekpath that aren't included in the path - they don't seem to generally get the BZ boundary right, but frequently (there probably needs to be a branch and more points in general). The code below illustrates the extra points in space group 3 (the "parent" space group of plane group 2):

using LinearAlgebra, Brillouin, Crystalline, PlotlyJS

Rs = directbasis(3)
Gs = reciprocalbasis(Rs)
cell = wignerseitz(Gs)
kp = irrfbz_path(3, Rs)

a, b, c, = norm(Rs[1]), norm(Rs[2]), norm(Rs[3])
cosβ = dot(Rs[3], Rs[1])/(c*a)
sinβ = sin(acos(cosβ))
Y = ((((1 + (a / c) * cosβ) / 2) / sinβ) / sinβ)
N = (1 / 2 + (Y * c * cosβ) / a)
pts = Dict(:Γ => [0, 0, 0], 
       :Z => [0, 1/2, 0], 
       :B => [0, 0, 1/2], 
       :Y₂ => [-1/2, 0, 0], 
       :C₂ => [-1/2, 1/2, 0], 
       :D => [0, 1/2, 1/2], 
       :A => [-1/2, 0, 1/2],
       :E => [-1/2, 1/2, 1/2],
       # extra points not included in the path (but in seekpath's point lists)
       :H => [-Y, 0, 1-N],
       :H₂ => [-1+Y, 0, N],
       :H₄ => [-Y, 0, -N],
       :M => [-Y, 1/2, 1-N],
       :M₂ => [-1+Y, 1/2, N],
       :M₄ => [-Y, 1/2, -N])

kp′ = deepcopy(kp)
merge!(kp′.points, pts)
plot(cell, kp′)

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Hi, thank you for using Spglib.jl in your code!
I've noticed that Brillouin.jl is currently using an older version of Spglib.jl. As the latest version of Spglib.jl includes several new features and bug fixes, it would be beneficial to update this dependency.

I've made several API changes in the latest updates. I understand these changes might require some adjustments in Brillouin.jl to ensure compatibility. Please let me know if you need any assistance with this update, I would be more than happy to help resolve the API changes and ensure Brillouin.jl works smoothly with the latest version of Spglib.jl.

Might be fun to implement https://doi.org/10.1038/s41524-020-00383-7 as an alternative to SeeK-path paths. The key appeal is that it is a general method that takes a lattice basis and a set of symmetry operations as input (as opposed to trying to enumerate/tabulate the relevant k-points explicitly).

A significant down-side, though, is that this wouldn't give us meaningful labels for the high-symmetry points.

I am currently experimenting with moving to PythonCall instead of using PyCall for python dependencies. It seems to have a few advantages (especially with respect to managing python dependencies from Julia), so I'm giving it a try. One unfortunate issue is, however, that using both PythonCall and PyCall in the same project can lead to issues (in the sense of a dying Julia session). For that reason I am wondering if there is an easy way to make Brillouin independent of requiring PyCall as a dependency. It seems to only be used in one submodule. How much effort would this be @thchr ?

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"))

Everything referenced from assets should probably be using https://github.com/JuliaPackaging/RelocatableFolders.jl.

Hi! We are using Brillouin.jl in our code DFTK.jl and since this morning update we are encountering Illegal instruction signals. I can't really manage to reproduce a minimal working example, sorry, but if we downgrade Brillouin.jl to 0.5.8, it works fine.

As far as I know, this should not happen so this might be linked e.g. to the usage of the spglib external library to generate the lattice we use as input to Brillouin.jl (here) or something like that, but maybe you already have encountered such errors and have some ideas on where it comes from ?

Thanks in advance :)

It seems like this package requires Julia v1.5 and above, is it necessary? Can we relax the version in Project.toml to v1.0?

It will be nice to have a function to find a periodic image of some point r inside the WS cell.
(https://discourse.julialang.org/t/ann-brillouin-jl-k-space-utilities-for-crystalline-eigenproblems/64597/5)

It seems that SymmetryReduceBZ.jl has a similar function, but I haven't looked closely.
https://jerjorg.github.io/SymmetryReduceBZ.jl/Documentation/#SymmetryReduceBZ.Symmetry.mapto_bz-Tuple{AbstractVector{var%22#s130%22}%20where%20var%22#s130%22%3C:Real,%20AbstractMatrix{var%22#s129%22}%20where%20var%22#s129%22%3C:Real,%20AbstractMatrix{var%22#s128%22}%20where%20var%22#s128%22%3C:Real,%20String}

Right now, wignerseitz will iterate to give points in a cartesian basis (but embeds the underlying basis) but irrfbz_path iterates to return points in the lattice basis (and doesn't embed the underlying basis - but can be converted via cartesianize!).

We should make these two consistent. Since the core utility for other applications will be irrfbz_path, and since they - I guess - would be more likely to want elements in a lattice basis, I guess it makes more sense to have everything in the lattice basis by default, and make it more ergonomic to convert to a cartesian basis (by referencing an embedded basis).

The downside to embedding the basis e.g. in a KPath is that we then need to be able to go between conventional and primitive basis vectors - which means either depending on Crystalline or vendoring parts of it.