diff --git a/toolbox/plotlyjs_wignerseitz.jl b/toolbox/plotlyjs_wignerseitz.jl index 0624519..a066798 100644 --- a/toolbox/plotlyjs_wignerseitz.jl +++ b/toolbox/plotlyjs_wignerseitz.jl @@ -2,8 +2,7 @@ using .PlotlyJS import .PlotlyJS: plot using LinearAlgebra: norm, dot using StaticArrays - -using .WignerSeitz: face_normal +using .WignerSeitz: face_normal, merge_coplanar! # ---------------------------------------------------------------------------------------- # # CONSTANTS @@ -118,8 +117,8 @@ function plot(c::Cell{3}, layout::Layout=Layout(); ts = vcat(tbz, tgs, tgtips, taxs, taxtips) return PlotlyJS.plot(ts, layout; config=config) end - -function plot(clist::Array{Cell{3}}, idx_center::Int, layout::Layout=Layout(); + +function plot(clist::Array{Cell{3}}, idx_center::Int, layout::Layout=Layout(); ibz=nothing, config::PlotConfig=PlotConfig(responsive=true, displaylogo=false)) layout = merge(DEFAULT_PLOTLY_LAYOUT_3D, layout) @@ -134,7 +133,38 @@ function plot(clist::Array{Cell{3}}, idx_center::Int, layout::Layout=Layout(); taxs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 6) taxtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 3) - + if !isnothing(ibz) + facecolor = repeat([ + "rgb(50, 200, 200)", + "rgb(100, 200, 255)", + "rgb(150, 200, 115)", + "rgb(200, 200, 50)", + "rgb(230, 200, 10)", + "rgb(255, 140, 0)" + ], inner=[2]) + setting(ibz) !== CARTESIAN && (ibz = cartesianize(ibz)) + faces = PlotlyJS.mesh3d(x=[ibz.verts[i][1] for i in 1:length(ibz.verts)], + y=[ibz.verts[i][2] for i in 1:length(ibz.verts)], + z=[ibz.verts[i][3] for i in 1:length(ibz.verts)], + i=[ibz.faces[i][1]-1 for i in 1:length(ibz.faces)], + j=[ibz.faces[i][2]-1 for i in 1:length(ibz.faces)], + k=[ibz.faces[i][3]-1 for i in 1:length(ibz.faces)], + facecolor = facecolor, opacity=0.7 + ) + merge_coplanar!(ibz) # Coplanar triangles have to be merged after calling mesh3d. + # BZ + tbz = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(ibz)+1) + for (i, poly) in enumerate(ibz) + tbz[i] = PlotlyJS.scatter3d( + x=push!(getindex.(poly, 1), poly[1][1]), + y=push!(getindex.(poly, 2), poly[1][2]), + z=push!(getindex.(poly, 3), poly[1][3]); + mode="lines", hovertext="Cell", hoverinfo="text+x+y+z", + line=attr(color=BZ_COL[], width=3)) + end + tbz[length(ibz)+1] = faces + push!(tbz_list, tbz) + end for (ci, c) in enumerate(clist) if (!isnothing(c)) setting(c) !== CARTESIAN && (c = cartesianize(c)) @@ -207,6 +237,7 @@ function plot(clist::Array{Cell{3}}, idx_center::Int, layout::Layout=Layout(); end end end + # combine traces and plot ts = tbz_list[1] for i in 2:length(tbz_list) diff --git a/toolbox/test_extend.jl b/toolbox/test_extend.jl index 7ca94a2..bb8ca6d 100644 --- a/toolbox/test_extend.jl +++ b/toolbox/test_extend.jl @@ -10,6 +10,9 @@ using BrillouinZoneMeshes using PyCall const PySpatial = PyNULL() using BrillouinZoneMeshes.LinearAlgebra +using SymmetryReduceBZ.Symmetry:calc_ibz +import SymmetryReduceBZ.Utilities: get_uniquefacets +import QHull include("default_colors.jl") include("plotlyjs_wignerseitz.jl") function wignerseitz_ext(basis::AVec{<:SVector{D,<:Real}}; @@ -41,6 +44,7 @@ function wignerseitz_ext(basis::AVec{<:SVector{D,<:Real}}; vor = PySpatial.Voronoi(lattice) # voronoi tesselation of lattice clist = Cell{D}[] idx_center_final = 0 + # grab all the vertices of the central voronoi region enclosing origo # for idx_cntr in 1:length(vor.point_region) for idx_cntr in 1:length(vor.point_region) @@ -69,6 +73,19 @@ function wignerseitz_ext(basis::AVec{<:SVector{D,<:Real}}; return clist, idx_center_final end +function convert_to_cell(hull::QHull.Chull{Float64}, basis::AVec{<:SVector{D,<:Real}}) where D + # The QHull julia package is used by SymmetryReduceBZ, whereas Brillouin package use Qhull from + # SciPy instead. Therefore we have to reload this function for julia QHull object. + vs′ = hull.points # vertices + simp = hvcat(size(hull.simplices,1),hull.simplices...)' + fs′ = simp # faces + + vs = SVector{D, Float64}.(eachrow(vs′)) + fs = Vector{Int}.(eachrow(fs′)) + return Cell(vs, fs, SVector{D, SVector{D, Float64}}(basis), Ref(CARTESIAN)) +end + + function wignerseitz_ext(basis::AVec{<:AVec{<:Real}}; kwargs...) D = length(first(basis)) all(V -> length(V) == D, @view basis[2:end]) || error(DomainError(basis, "provided `basis` must have identical dimensions")) @@ -77,6 +94,7 @@ function wignerseitz_ext(basis::AVec{<:AVec{<:Real}}; kwargs...) end + # Wigner-Seitz cells visualization # 2D #Rs = [[1.0, 0.0], [-0.5, √3/2]] @@ -89,22 +107,42 @@ end # latvec = [1 0; 0 1]' # 3D -lattice = [[0 1 1.0]; [1 0 1.0]; [1 1 0.0]] -atoms = [1, 1] -positions = [ones(3) / 8, -ones(3) / 8] -br = BZMeshes.Brillouin(lattice=lattice, atoms=atoms, positions=positions) +lattice = [[0 0.5 0.5]; [0.5 0 0.5]; [0.5 0.5 0.0]] +#atoms = [1,1] +#positions = [ones(3) / 8, -ones(3) / 8] + #lattice = [[1.0 0.0 0.0]; [0.0 1.0 0.0]; [0.0 0.0 1.0]] +atoms = [1] +positions = [zeros(3)] + +atom_pos = hvcat(size(positions,1),positions...) +ibzformat = "convex hull" +coordinates = "Cartesian" +makeprim = true +convention = "ordinary" + +br = BZMeshes.Brillouin(lattice=lattice, atoms=atoms, positions=positions) #br = BZMeshes.Brillouin(lattice=lattice) msize = (4, 4, 4) bzmesh = UniformBZMesh(br=br, size=msize) - meshmap = MeshMap(bzmesh) -latvec = mapslices(x -> [x], lattice_vector(bzmesh), dims=1)[:] -clist, idx_center = wignerseitz_ext(latvec, cut=1) +recip_lattice=lattice_vector(bzmesh) +#latvec = mapslices(x -> [x],recip_lattice , dims=1)[:] +latvec = mapslices(x -> [x], recip_lattice , dims=1)[:] + +#generate irreducible brillouin zone +ibz = calc_ibz(lattice./2/π,atoms, atom_pos, coordinates,ibzformat, + makeprim,convention) +#The 1/2π here is due to the convention difference between packages. SymmetryReduceBZ has a \dot a_recip = 1 instead of 2π -P = plot(clist, idx_center) +c = convert_to_cell(ibz, SVector{length(latvec[1]),Float64}.(latvec)) +c = WignerSeitz.reorient_normals!(c) +latticize!(c) # Do not merge coplanar triangles for reduced mesh. This must be done within plot. + +clist, idx_center = wignerseitz_ext(latvec, cut=1) +P = plot(clist, idx_center, ibz=c) fullmesh = [bzmesh[i] for i in 1:length(bzmesh)] reducedmesh = [bzmesh[i] for i in meshmap.irreducible_indices]