From 833108bae26fd91292577bf3bb7574783df92d4c Mon Sep 17 00:00:00 2001 From: Tao Wang Date: Wed, 26 Oct 2022 16:40:20 -0400 Subject: [PATCH 1/4] add plot --- toolbox/default_colors.jl | 19 ++ toolbox/plotlyjs_wignerseitz.jl | 498 ++++++++++++++++++++++++++++++++ toolbox/test_extend.jl | 97 +++++++ 3 files changed, 614 insertions(+) create mode 100644 toolbox/default_colors.jl create mode 100644 toolbox/plotlyjs_wignerseitz.jl create mode 100644 toolbox/test_extend.jl diff --git a/toolbox/default_colors.jl b/toolbox/default_colors.jl new file mode 100644 index 0000000..e7d17bc --- /dev/null +++ b/toolbox/default_colors.jl @@ -0,0 +1,19 @@ +const TRANSPARENT_COL = Ref("rgba(255, 255, 255, 1)") + +# --- wigner-seitz colors --- +# from https://flatuicolors.com/palette/gb +const BZ_COL = Ref("rgb(47,54,64)") # "electromagnetic" +const BASIS_COL = Ref("rgb(39,60,117)") # "pico void" +const BASIS_LIGHT_COL = Ref("rgb(212,216,227)") # 20% BASIS_COL, 80% white +const AXIS_COL = Ref("rgb(194,54,22)") # "harley davidson orange" +const AXIS_LIGHT_COL = Ref("rgb(242,215,208)") # 20% AXIS_COL, 80% white + +# --- k-path colors --- +# from https://flatuicolors.com/palette/ca +const KPATH_COL = Ref("rgb(95,39,205)") # "nasu purple" + +# --- dispersion colors --- +# from https://flatuicolors.com/palette/gb +const BAND_COL = Ref("rgb(39,60,117)") # "mazarine blue" +const KLINE_COL = Ref("rgb(220,221,225)") # "hint of pensive" (light gray) +const ANNOTATE_COL = Ref("rgb(53, 59, 72)") # "blue nights" (dark gray) \ No newline at end of file diff --git a/toolbox/plotlyjs_wignerseitz.jl b/toolbox/plotlyjs_wignerseitz.jl new file mode 100644 index 0000000..8305f68 --- /dev/null +++ b/toolbox/plotlyjs_wignerseitz.jl @@ -0,0 +1,498 @@ +using .PlotlyJS +import .PlotlyJS: plot +using LinearAlgebra: norm, dot +using StaticArrays + +using .WignerSeitz: face_normal +# ---------------------------------------------------------------------------------------- # +# CONSTANTS + +# default layout +const DEFAULT_PLOTLY_LAYOUT_3D = Layout( + showlegend=false, + scene=attr( + xaxis=attr(tickvals=[], zeroline=false, + showgrid=false, showbackground=false, + title=attr(text=""), + ), + yaxis=attr(tickvals=[], zeroline=false, + showgrid=false, showbackground=false, + title=attr(text=""), + ), + zaxis=attr(tickvals=[], zeroline=false, + showgrid=false, showbackground=false, + title=attr(text=""), + ), + aspectmode = "data", + camera=attr(up = attr(x=0, z=1, y=0), center = attr(x=0, y=0, z=0)), + dragmode = "turntable", + ), + margin=attr(l=0, r=0, b=0, t=0), + autosize=false, + #width=200, height=200, + plot_bgcolor=TRANSPARENT_COL[], paper_bgcolor=TRANSPARENT_COL[], + ) + +# ---------------------------------------------------------------------------------------- # +# 3D + +function plot(c::Cell{3}, layout::Layout = Layout(); + config::PlotConfig = PlotConfig(responsive=true, displaylogo=false)) + + layout = merge(DEFAULT_PLOTLY_LAYOUT_3D, layout) + + setting(c) !== CARTESIAN && (c = cartesianize(c)) + scale = maximum(norm, basis(c)) + + # BZ + tbz = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(c)) + for (i,poly) in enumerate(c) + 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 + + # lattice vectors + tgs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 6) + tgtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 3) + intersects = axis_intersections(c, basis(c)) + for (i,V) in enumerate(basis(c)) + V′ = V./norm(V) + name = "v$(i)" + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : 1.0 + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : 0.0) *scale*V′ + tgs[i+(j-1)*3] = PlotlyJS.scatter3d( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]], z=[V₀[3],V₁[3]]; + mode="lines", hovertext=name, hoverinfo="text", + line=attr(color=ifelse(j==1, BASIS_LIGHT_COL[], BASIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + end + tgtips[i] = PlotlyJS.cone( + x=[V[1]], u=[V′[1]], y=[V[2]], v=[V′[2]], z=[V[3]], w=[V′[3]], + sizeref=.1*scale, showscale=false, anchor="tail", + colorscale=[[0, BASIS_COL[]], [1, BASIS_COL[]]], + hovertext=name, hoverinfo="text+x+y+z") + end + + # Cartesian axes + cart_basis = cartesian_axes(Val(3)) + name_axs = ("x", "y", "z") + taxs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 6) + taxtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 3) + intersects = axis_intersections(c, cart_basis) + for (i, V) in enumerate(cart_basis) + name = ""*string(name_axs[i])*"" + V′ = V./norm(V) + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : intersects[i] + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : -0.2) *scale*V′ + taxs[i+(j-1)*3] = PlotlyJS.scatter3d( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]], z=[V₀[3],V₁[3]]; + mode="lines", hovertext=name, hoverinfo="text", + line=attr(color=ifelse(j==1, AXIS_LIGHT_COL[], AXIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + if j == 2 + taxtips[i] = PlotlyJS.cone( + x=[V₁[1]], u=[V[1]], + y=[V₁[2]], v=[V[2]], + z=[V₁[3]], w=[V[3]], + sizeref=.1*scale, showscale=false, anchor="tail", + colorscale=[[0, AXIS_COL[]], [1, AXIS_COL[]]], + hovertext=name, hoverinfo="text") + end + end + end + + # combine traces and plot + 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(); + config::PlotConfig = PlotConfig(responsive=true, displaylogo=false)) + + layout = merge(DEFAULT_PLOTLY_LAYOUT_3D, layout) + + # scale = maximum(norm, basis(c)) + c_center = clist[idx_center] + setting(c_center) !== CARTESIAN && (c_center = cartesianize(c_center)) + scale = maximum(norm, basis(c_center)) + tbz_list=[] + tgs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 6) + tgtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 3) + taxs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 6) + taxtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 3) + + + for (ci,c) in enumerate(clist) + if(!isnothing(c)) + setting(c) !== CARTESIAN && (c = cartesianize(c)) + # BZ + tbz = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(c)) + for (i,poly) in enumerate(c) + 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 + push!(tbz_list,tbz) + if ci == idx_center + # lattice vectors + intersects = axis_intersections(c, basis(c)) + for (i,V) in enumerate(basis(c)) + V′ = V./norm(V) + name = "v$(i)" + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : 1.0 + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : 0.0) *scale*V′ + tgs[i+(j-1)*3] = PlotlyJS.scatter3d( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]], z=[V₀[3],V₁[3]]; + mode="lines", hovertext=name, hoverinfo="text", + line=attr(color=ifelse(j==1, BASIS_LIGHT_COL[], BASIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + end + tgtips[i] = PlotlyJS.cone( + x=[V[1]], u=[V′[1]], y=[V[2]], v=[V′[2]], z=[V[3]], w=[V′[3]], + sizeref=.1*scale, showscale=false, anchor="tail", + colorscale=[[0, BASIS_COL[]], [1, BASIS_COL[]]], + hovertext=name, hoverinfo="text+x+y+z") + end + + # Cartesian axes + cart_basis = cartesian_axes(Val(3)) + name_axs = ("x", "y", "z") + intersects = axis_intersections(c, cart_basis) + for (i, V) in enumerate(cart_basis) + name = ""*string(name_axs[i])*"" + V′ = V./norm(V) + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : intersects[i] + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : -0.2) *scale*V′ + taxs[i+(j-1)*3] = PlotlyJS.scatter3d( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]], z=[V₀[3],V₁[3]]; + mode="lines", hovertext=name, hoverinfo="text", + line=attr(color=ifelse(j==1, AXIS_LIGHT_COL[], AXIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + if j == 2 + taxtips[i] = PlotlyJS.cone( + x=[V₁[1]], u=[V[1]], + y=[V₁[2]], v=[V[2]], + z=[V₁[3]], w=[V[3]], + sizeref=.1*scale, showscale=false, anchor="tail", + colorscale=[[0, AXIS_COL[]], [1, AXIS_COL[]]], + hovertext=name, hoverinfo="text") + end + end + end + end + end + end + # combine traces and plot + ts = tbz_list[1] + for i in 2:length(tbz_list) + ts = vcat(ts,tbz_list[i]) + end + ts = vcat(ts,tgs,tgtips,taxs,taxtips) + return PlotlyJS.plot(ts, layout; config=config) +end + + + +# ---------------------------------------------------------------------------------------- # +# 2D + +# default layout +const DEFAULT_PLOTLY_LAYOUT_2D = Layout( + showlegend=false, + xaxis=attr(tickvals=[], zeroline=false, + showgrid=false, showbackground=false, + title=attr(text=""), + ), + yaxis=attr(tickvals=[], zeroline=false, + showgrid=false, showbackground=false, + title=attr(text=""), + scaleanchor="x", scaleratio=1 + ), + aspectmode = "data", + hovermode="closest", + margin=attr(l=0, r=0, b=0, t=0), + autosize=false, + plot_bgcolor=TRANSPARENT_COL[], paper_bgcolor=TRANSPARENT_COL[], + annotations=PlotlyBase.PlotlyAttribute[] + ) + +function plot(c::Cell{2}, layout::Layout = Layout(); + config::PlotConfig = PlotConfig(responsive=true, displaylogo=false)) + + layout = merge(DEFAULT_PLOTLY_LAYOUT_2D, layout) + + setting(c) !== CARTESIAN && (c = cartesianize(c)) + + scale = maximum(norm, basis(c)) + max_x, max_y = maximum(v->abs(v[1]), basis(c)), maximum(v->abs(v[2]), basis(c)) + get!(layout[:xaxis], :range, [-max_x-scale/15, max_x+scale/15]) + get!(layout[:yaxis], :range, [-max_y-scale/15, max_y+scale/15]) + + # Cell boundaries + tbz = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(c)) + for (i,poly) in enumerate(c) + tbz[i] = PlotlyJS.scatter( + x=push!(getindex.(poly, 1), poly[1][1]), + y=push!(getindex.(poly, 2), poly[1][2]); + mode="lines", hovertext="Cell", hoverinfo="text+x+y", + line=attr(color=BZ_COL[], width=3) + ) + end + + # lattice vectors + tgs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 4) + tgtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 2) + intersects = axis_intersections(c, basis(c)) + for (i,V) in enumerate(basis(c)) + V′ = V./norm(V) + name = "v$(i)" + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : 1.0 + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : 0.0) *scale*V′ + tgs[i+(j-1)*2] = PlotlyJS.scatter( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]]; + mode="lines", hovertext=name, hoverinfo=ifelse(j==1,"text","text+x+y"), + line=attr(color=ifelse(j==1, BASIS_LIGHT_COL[], BASIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + end + # arrow heads have to be added as annotations to layout in 2D :/ + haskey(layout, :annotations) || (layout[:annotations] = PlotlyBase.PlotlyAttribute[]) + push!(layout[:annotations], + attr(x=V[1]+.05V′[1]*scale, y=V[2]+.05V′[2]*scale, # awful fidgeting; plotly's + ax=V[1]-.05V′[1]*scale, ay=V[2]-.05V′[2]*scale, # arrows are stupid + xref="ax", yref="ay", axref="x", ayref="y", + showarrow=true, arrowhead=2, arrowwidth=6, arrowsize=.5, + arrowcolor=BASIS_COL[])) + end + + # Cartesian axes + cart_basis = cartesian_axes(Val(2)) + name_axs = ("x", "y") + taxs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 4) + taxtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 2) + intersects = axis_intersections(c, cart_basis) + for (i, V) in enumerate(cart_basis) + name = ""*string(name_axs[i])*"" + V′ = V./norm(V) + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : intersects[i] + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : -0.2) *scale*V′ + + taxs[i+2(j-1)] = PlotlyJS.scatter( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]]; + mode="lines", hovertext=name, hoverinfo="text", + line=attr(color=ifelse(j==1, AXIS_LIGHT_COL[], AXIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + if j == 2 + # arrow heads have to be added as annotations to layout in 2D :/ + haskey(layout, :annotations) || (layout[:annotations] = PlotlyBase.PlotlyAttribute[]) + push!(layout[:annotations], + attr(x=V₁[1]+.05V′[1]*scale, y=V₁[2]+.05V′[2]*scale, # awful fidgeting; plotly's + ax=V₁[1]-.05V′[1]*scale, ay=V₁[2]-.05V′[2]*scale, # arrows are stupid + xref="ax", yref="ay", axref="x", ayref="y", + showarrow=true, arrowhead=2, arrowwidth=6, arrowsize=.5, + arrowcolor=AXIS_COL[])) + end + end + end + + # combine traces and plot + ts = vcat(tbz, tgs, taxs) + return PlotlyJS.plot(ts, layout; config=config) +end + +function plot(clist::Array{Cell{2}}, idx_center::Int, layout::Layout = Layout(); + config::PlotConfig = PlotConfig(responsive=true, displaylogo=false)) + + layout = merge(DEFAULT_PLOTLY_LAYOUT_2D, layout) + tbz_list = [] + tgs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 4) + tgtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 2) + taxs = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 4) + taxtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 2) + + c_center = clist[idx_center] + setting(c_center) !== CARTESIAN && (c_center = cartesianize(c_center)) + scale = maximum(norm, basis(c_center)) + max_x, max_y = maximum(v->abs(v[1]), basis(c_center)), maximum(v->abs(v[2]), basis(c_center)) + get!(layout[:xaxis], :range, [-max_x-scale/15, max_x+scale/15]) + get!(layout[:yaxis], :range, [-max_y-scale/15, max_y+scale/15]) + for (ci,c) in enumerate(clist) + if(!isnothing(c)) + # Cell boundaries + setting(c) !== CARTESIAN && (c = cartesianize(c)) + tbz = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(c)) + for (i,poly) in enumerate(c) + tbz[i] = PlotlyJS.scatter( + x=push!(getindex.(poly, 1), poly[1][1]), + y=push!(getindex.(poly, 2), poly[1][2]); + mode="lines", hovertext="Cell", hoverinfo="text+x+y", + line=attr(color=BZ_COL[], width=3) + ) + end + push!(tbz_list,tbz) + # lattice vectors + if ci == idx_center + + intersects = axis_intersections(c, basis(c)) + for (i,V) in enumerate(basis(c)) + V′ = V./norm(V) + name = "v$(i)" + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : 1.0 + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : 0.0) *scale*V′ + tgs[i+(j-1)*2] = PlotlyJS.scatter( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]]; + mode="lines", hovertext=name, hoverinfo=ifelse(j==1,"text","text+x+y"), + line=attr(color=ifelse(j==1, BASIS_LIGHT_COL[], BASIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + end + # arrow heads have to be added as annotations to layout in 2D :/ + haskey(layout, :annotations) || (layout[:annotations] = PlotlyBase.PlotlyAttribute[]) + push!(layout[:annotations], + attr(x=V[1]+.05V′[1]*scale, y=V[2]+.05V′[2]*scale, # awful fidgeting; plotly's + ax=V[1]-.05V′[1]*scale, ay=V[2]-.05V′[2]*scale, # arrows are stupid + xref="ax", yref="ay", axref="x", ayref="y", + showarrow=true, arrowhead=2, arrowwidth=6, arrowsize=.5, + arrowcolor=BASIS_COL[])) + end + + # Cartesian axes + cart_basis = cartesian_axes(Val(2)) + name_axs = ("x", "y") + intersects = axis_intersections(c, cart_basis) + for (i, V) in enumerate(cart_basis) + name = ""*string(name_axs[i])*"" + V′ = V./norm(V) + for j in (1,2) # in/outside BZ + start = j == 1 ? 0.0 : intersects[i] + stop = j == 1 ? intersects[i] : intersects[i] + V₀ = start*V .+ (j == 1 ? 0.0 : 0.025)*scale*V′ + V₁ = stop*V .- (j == 1 ? 0.025 : -0.2) *scale*V′ + + taxs[i+2(j-1)] = PlotlyJS.scatter( + x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]]; + mode="lines", hovertext=name, hoverinfo="text", + line=attr(color=ifelse(j==1, AXIS_LIGHT_COL[], AXIS_COL[]), + width=ifelse(j==1, 5, 6)) + ) + if j == 2 + # arrow heads have to be added as annotations to layout in 2D :/ + haskey(layout, :annotations) || (layout[:annotations] = PlotlyBase.PlotlyAttribute[]) + push!(layout[:annotations], + attr(x=V₁[1]+.05V′[1]*scale, y=V₁[2]+.05V′[2]*scale, # awful fidgeting; plotly's + ax=V₁[1]-.05V′[1]*scale, ay=V₁[2]-.05V′[2]*scale, # arrows are stupid + xref="ax", yref="ay", axref="x", ayref="y", + showarrow=true, arrowhead=2, arrowwidth=6, arrowsize=.5, + arrowcolor=AXIS_COL[])) + end + end + end + # combine traces and plot + #ts = vcat(tbz, tgs, taxs) + end + end + end + #ts=vcat(tb for tb in tbz_list) + ts = tbz_list[1] + for i in 2:length(tbz_list) + ts = vcat(ts,tbz_list[i]) + end + ts = vcat(ts,tgs,taxs) + #ts = tbz_list[1] + return PlotlyJS.plot(ts,layout; config=config) +end + + + + +# ---------------------------------------------------------------------------------------- # +# UTILITIES + +# the "outward" intersections of of lines with direction `Vs` though origo with a cell face +function axis_intersections(c::Cell{3}, + Vs::AbstractVector{<:AbstractVector}=cartesian_axes(Val(3))) + + intersects = MVector{3,Float64}(undef) + fill!(intersects, Inf) + # intersection between plane (r-rₒ) ⋅ n = 0 and line r(t) = l₀ + lt: + # t = (r₀ - l₀) ⋅ n / (l ⋅ n) [https://wikipedia.org/wiki/Line–plane_intersection] + for (i,poly) in enumerate(c) + r₀ = sum(poly)/length(poly) + n = face_normal(c, i) + for (j,V) in enumerate(Vs) + t = dot(r₀, n)/dot(V, n) + if t > 0.0 && t < intersects[j] + intersects[j] = t + end + end + end + return intersects +end + +# the "outward" intersections of of lines with direction `Vs` though origo with a cell segment +function axis_intersections(c::Cell{2}, + Vs::AbstractVector{<:AbstractVector}=cartesian_axes(Val(2))) + + intersects = MVector{2,Float64}(undef) + fill!(intersects, Inf) + # intersection between line rₐ(t) a + nₐt and line rᵦ(t) = β + nᵦt can be gotten by + # linear algebra: [-nₐ|nᵦ][tₐ,tᵦ] = a-β. Here we, pick rᵦ for the axes `Vs` (β=0) and + # let rₐ refer to cell segments. We're then only interested in tᵦ: + fs = faces(c) + vs = vertices(c) + for i in eachindex(vs) # number of vertices and segments are equal in 2D + if length(fs) == 1 # `merge = true` + idxs = i ≠ length(vs) ? (i,i+1) : (i,1) + else + idxs = (fs[i][1], fs[i][2]) + end + a = vs[idxs[1]] + nₐ = vs[idxs[2]] - a + for (j,V) in enumerate(Vs) + nᵦ = V + tᵦ = (nₐ[2]*a[1]-nₐ[1]*a[2])/(-nₐ[1]*nᵦ[2] + nₐ[2]*nᵦ[1]) # ~ inv([-nₐ|nᵦ]) + if tᵦ > 0.0 && tᵦ < intersects[j] + intersects[j] = tᵦ + end + end + end + return intersects +end + +# a set of cartesian basis vectors in D dimensions +cartesian_axes(Dᵛ::Val{D}) where D = SVector(ntuple(i->SVector(ntuple(j->i==j ? 1.0 : 0.0, Dᵛ)), Dᵛ)) diff --git a/toolbox/test_extend.jl b/toolbox/test_extend.jl new file mode 100644 index 0000000..2c1be14 --- /dev/null +++ b/toolbox/test_extend.jl @@ -0,0 +1,97 @@ +using Brillouin, PlotlyJS +using Brillouin: + AVec, + CARTESIAN +import Brillouin: + basis +using Brillouin.KPaths: Bravais +using StaticArrays + +using PyCall +const PySpatial = PyNULL() + +include("default_colors.jl") +include("plotlyjs_wignerseitz.jl") +function wignerseitz_ext(basis::AVec{<:SVector{D,<:Real}}; + merge::Bool = true, + Nmax::Integer = 3, cut = Nmax-1) where D + # bringing in SciPy's Spatial module (for `Voronoi` and `ConvexHull`) + if PyCall.conda + copy!(PySpatial, pyimport_conda("scipy.spatial", "scipy")) + else + copy!(PySpatial, pyimport_e("scipy.spatial")) + end + ispynull(PySpatial) && @warn("scipy python package not found. " * + "WignerSeitz.wignerseitz is nonfunctional.") + + # "supercell" lattice of G-vectors + Ns = -Nmax:Nmax + lattice = Vector{SVector{D,Float64}}(undef, length(Ns)^D) + idx_cntr = 0 + idx_center = 0 + idx_cartesian = CartesianIndices(ntuple(_->Ns, Val(D))) + for (idx, I) in enumerate(idx_cartesian) + V = basis[1]*I[1] + D >= 2 && (V += basis[2]*I[2]) + D == 3 && (V += basis[3]*I[3]) + lattice[idx] = V + iszero(I) && (idx_center = idx) + end + ispynull(PySpatial) && error("You need to install scipy for wignerseitz to work.") + 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) + if sqrt(sum((idx_cartesian[idx_cntr][i]-idx_cartesian[idx_center][i])^2 for i in 1:D))0] + # get convex hull of central vertices + if length(verts_cntr)>D+1 + hull = PySpatial.ConvexHull(verts_cntr) + c = WignerSeitz.convert_to_cell(hull, basis) + c = WignerSeitz.reorient_normals!(c) + + # return either raw simplices or "merged" polygons + if merge + latticize!(merge_coplanar!(c)) + else + latticize!(c) + end + push!(clist,c) + if idx_cntr==idx_center + idx_center_final =length(clist) + end + end + end + end + return clist, idx_center_final +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")) + + return wignerseitz_ext(SVector{D,Float64}.(basis); kwargs...) +end + + +# Wigner-Seitz cells visualization +# 3D +Rs = [[1.0, 0.0], [-0.5, √3/2]] +#Rs = 2π.*[[-1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, -1.0]] +clist,idx_center = wignerseitz_ext(Rs;merge = false) +#P0=plot(clist[idx_center]) +#display(P0) +#readline() +v=[1.0,1.0] +vc = cartesianize(v,Rs) +#println(vc[1:1],vc[2:2]) +linspace = range(0,1,length=30) +lin2 = collect(Iterators.product(linspace,linspace)) +P=plot(clist,idx_center) +addtraces!(P, scatter(x=[vc[1].*lin2[i][1] for i in 1:length(lin2)],y=[vc[2].*lin2[i][2] for i in 1:length(lin2)],mode="markers")) +display(P) +readline() + From b0dbba7c0923b530f261e23b01ba0d967eddc283 Mon Sep 17 00:00:00 2001 From: Tao Wang Date: Sat, 29 Oct 2022 00:05:12 -0400 Subject: [PATCH 2/4] plot using UniformBZMesh --- src/Model.jl | 4 ++-- toolbox/test_extend.jl | 27 +++++++++++++-------------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/src/Model.jl b/src/Model.jl index 980bdbc..a46610c 100644 --- a/src/Model.jl +++ b/src/Model.jl @@ -22,7 +22,7 @@ We use the convention that the reciprocal lattice is the set of G vectors such that G ⋅ R ∈ 2π ℤ for all R in the lattice. """ function _compute_recip_lattice(lattice::Matrix{T}) where {T} - return 2T(π) * _compute_inverse_lattice(lattice) + return 2T(π) * _compute_inverse_lattice(Matrix(lattice')) end """ @@ -303,4 +303,4 @@ function Base.show(io::IO, ::MIME"text/plain", model::Brillouin{T,DIM}) where {T end end -end \ No newline at end of file +end diff --git a/toolbox/test_extend.jl b/toolbox/test_extend.jl index 2c1be14..5fcabf0 100644 --- a/toolbox/test_extend.jl +++ b/toolbox/test_extend.jl @@ -6,10 +6,10 @@ import Brillouin: basis using Brillouin.KPaths: Bravais using StaticArrays - +using BrillouinZoneMeshes using PyCall const PySpatial = PyNULL() - +using BrillouinZoneMeshes.LinearAlgebra include("default_colors.jl") include("plotlyjs_wignerseitz.jl") function wignerseitz_ext(basis::AVec{<:SVector{D,<:Real}}; @@ -78,20 +78,19 @@ end # Wigner-Seitz cells visualization -# 3D +# 2D Rs = [[1.0, 0.0], [-0.5, √3/2]] -#Rs = 2π.*[[-1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, -1.0]] -clist,idx_center = wignerseitz_ext(Rs;merge = false) -#P0=plot(clist[idx_center]) -#display(P0) -#readline() -v=[1.0,1.0] -vc = cartesianize(v,Rs) -#println(vc[1:1],vc[2:2]) -linspace = range(0,1,length=30) -lin2 = collect(Iterators.product(linspace,linspace)) +N, DIM = 4, 2 +origin = [0.0 0.0] +lattice =Matrix([2 0; 1 sqrt(3)]') +msize = (3,3) +# latvec = [1 0; 0 1]' +br = BZMeshes.Brillouin(lattice = lattice) +bzmesh = UniformBZMesh(br=br, size=msize) +latvec = mapslices(x->[x],lattice_vector(bzmesh),dims=1)[:] +clist,idx_center = wignerseitz_ext(latvec;merge = false) P=plot(clist,idx_center) -addtraces!(P, scatter(x=[vc[1].*lin2[i][1] for i in 1:length(lin2)],y=[vc[2].*lin2[i][2] for i in 1:length(lin2)],mode="markers")) +addtraces!(P, scatter(x=[bzmesh[i][1] for i in 1:length(bzmesh)],y=[bzmesh[i][2] for i in 1:length(bzmesh)],mode="markers")) display(P) readline() From 00ab1956c6f95139163c3b5263767617a0770799 Mon Sep 17 00:00:00 2001 From: Tao Wang Date: Sat, 29 Oct 2022 00:30:21 -0400 Subject: [PATCH 3/4] fix test --- test/Model.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/Model.jl b/test/Model.jl index 9819589..a15f529 100644 --- a/test/Model.jl +++ b/test/Model.jl @@ -7,7 +7,7 @@ using BrillouinZoneMeshes.PointSymmetry: spglib_spacegroup_number, spglib_standa DIM = 2 lattice = Matrix([1.0 0; 0 1]') br = Brillouin(lattice=lattice) - @test br.inv_lattice .* 2π ≈ br.recip_lattice + @test br.inv_lattice .* 2π ≈ br.recip_lattice' @test br.unit_cell_volume ≈ abs(det(lattice)) @test br.recip_cell_volume ≈ 1 / abs(det(lattice)) * (2π)^DIM @@ -15,7 +15,7 @@ using BrillouinZoneMeshes.PointSymmetry: spglib_spacegroup_number, spglib_standa DIM = 2 lattice = Matrix([2.0 0; 1 sqrt(3)]') br = Brillouin(lattice=lattice) - @test br.inv_lattice .* 2π ≈ br.recip_lattice + @test br.inv_lattice .* 2π ≈ br.recip_lattice' @test br.unit_cell_volume ≈ abs(det(lattice)) @test br.recip_cell_volume ≈ 1 / abs(det(lattice)) * (2π)^DIM @@ -23,7 +23,7 @@ using BrillouinZoneMeshes.PointSymmetry: spglib_spacegroup_number, spglib_standa DIM = 3 lattice = Matrix([2.0 0 0; 1 sqrt(3) 0; 7 11 19]') br = Brillouin(lattice=lattice) - @test br.inv_lattice .* 2π ≈ br.recip_lattice + @test br.inv_lattice .* 2π ≈ br.recip_lattice' @test br.unit_cell_volume ≈ abs(det(lattice)) @test br.recip_cell_volume ≈ 1 / abs(det(lattice)) * (2π)^DIM end From 0376490c5d0ce11be154ca652b42805b7547365d Mon Sep 17 00:00:00 2001 From: Tao Wang Date: Sat, 29 Oct 2022 11:09:04 -0400 Subject: [PATCH 4/4] minor --- src/Model.jl | 8 ++++---- test/Model.jl | 10 ++++++++-- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/Model.jl b/src/Model.jl index a46610c..5039324 100644 --- a/src/Model.jl +++ b/src/Model.jl @@ -11,7 +11,7 @@ export Brillouin, get_latvec """ Compute the inverse of the lattice. Require lattice to be square matrix """ -function _compute_inverse_lattice(lattice::Matrix{T}) where {T} +function _compute_inverse_lattice(lattice::AbstractMatrix{T}) where {T} @assert size(lattice, 2) == size(lattice, 2) return inv(lattice) end @@ -21,15 +21,15 @@ Compute the reciprocal lattice. We use the convention that the reciprocal lattice is the set of G vectors such that G ⋅ R ∈ 2π ℤ for all R in the lattice. """ -function _compute_recip_lattice(lattice::Matrix{T}) where {T} - return 2T(π) * _compute_inverse_lattice(Matrix(lattice')) +function _compute_recip_lattice(lattice::AbstractMatrix{T}) where {T} + return 2T(π) * _compute_inverse_lattice(lattice') end """ Return I-th lattice vector of lattice. Lattice vectors are specified column-wise in lattice::Matrix. """ -function get_latvec(lattice::Matrix{T}, I::Int) where {T} +function get_latvec(lattice::AbstractMatrix{T}, I::Int) where {T} return view(lattice, :, I) end diff --git a/test/Model.jl b/test/Model.jl index a15f529..72feb44 100644 --- a/test/Model.jl +++ b/test/Model.jl @@ -2,7 +2,7 @@ using BrillouinZoneMeshes.PointSymmetry: spglib_spacegroup_number, spglib_standa @testset "Brillouin" begin Brillouin = BrillouinZoneMeshes.Model.Brillouin - + get_latvec = BrillouinZoneMeshes.Model.get_latvec # square lattice DIM = 2 lattice = Matrix([1.0 0; 0 1]') @@ -18,7 +18,9 @@ using BrillouinZoneMeshes.PointSymmetry: spglib_spacegroup_number, spglib_standa @test br.inv_lattice .* 2π ≈ br.recip_lattice' @test br.unit_cell_volume ≈ abs(det(lattice)) @test br.recip_cell_volume ≈ 1 / abs(det(lattice)) * (2π)^DIM - + for i in 1:DIM + @test dot(get_latvec(br.recip_lattice,i),get_latvec(br.lattice,i)) ≈ 2π + end # 3d testing lattice DIM = 3 lattice = Matrix([2.0 0 0; 1 sqrt(3) 0; 7 11 19]') @@ -26,6 +28,10 @@ using BrillouinZoneMeshes.PointSymmetry: spglib_spacegroup_number, spglib_standa @test br.inv_lattice .* 2π ≈ br.recip_lattice' @test br.unit_cell_volume ≈ abs(det(lattice)) @test br.recip_cell_volume ≈ 1 / abs(det(lattice)) * (2π)^DIM + for i in 1:DIM + @test dot(get_latvec(br.recip_lattice,i),get_latvec(br.lattice,i)) ≈ 2π + end + end @testset "Standard Brillouin" begin