diff --git a/src/Crystalline.jl b/src/Crystalline.jl index fe613b6e..a47a3a90 100644 --- a/src/Crystalline.jl +++ b/src/Crystalline.jl @@ -32,15 +32,17 @@ using .SquareStaticMatrices # exports `SSqMatrix{D,T}` # include vendored SmithNormalForm.jl package from ../.vendor/ include("../.vendor/SmithNormalForm/src/SmithNormalForm.jl") -using .SmithNormalForm -import .SmithNormalForm: smith, Smith # TODO: remove explicit import when we update SmithNormalForm -export smith, Smith # export, so that loading Crystalline effectively also provides SmithNormalForm +using .SmithNormalForm: smith, Smith +export smith, Smith # export, so that loading Crystalline also namespaces these @reexport using Bravais import Bravais: primitivize, conventionalize, cartesianize, transform, centering using Bravais: stack, all_centeringtranslations, centeringtranslation, centering_volume_fraction +include("surface_unitcell.jl") +export surface_basis # TODO: move to Bravais (but tricky cf. SmithNormalForm dependency) + # included files and exports include("constants.jl") export MAX_SGNUM, MAX_SUBGNUM, MAX_MSGNUM, MAX_MSUBGNUM, ENANTIOMORPHIC_PAIRS diff --git a/src/surface_unitcell.jl b/src/surface_unitcell.jl new file mode 100644 index 00000000..2255f90f --- /dev/null +++ b/src/surface_unitcell.jl @@ -0,0 +1,356 @@ +using Bravais +using StaticArrays +using LinearAlgebra: dot, cross, norm, normalize, I +using .SmithNormalForm: smith + +# what are the constraints we impose on surface-lattice points {r}? they are: +# (1) n ⋅ r = 0 for surface lattice vectors r +# (2) r = ∑ⱼ aⱼRⱼ +# with integer coefficients aⱼ (j = 1,2,3); combining (1) & (2), we have +# n ⋅ ∑ⱼ aⱼRⱼ = 0 ⇔ ∑ⱼ aⱼ(n ⋅ Rⱼ) = 0 +function surface_basis_from_miller_indices_via_snf( + n :: AbstractVector{<:Integer}; + # surface's Miller indices (normal vector ∝ n[1]*Gs[1] + n[2]*Gs[2] + n[3]*Gs[3]) + ) + length(Rs) == 3 || error("only 3D implemented; number of basis vectors is not 3") + all(R->length(R) == 3, Rs) || error("only 3D implemented; basis vectors are not 3D") + length(n) == 3 || error("only 3D implemented; normal vector is not 3D") + + # build a matrix to represent the relation n ⋅ r = 0; since we specify n in the a basis + # of the reciprocal lattice vectors, and r in the direct basis, the dot product is + # simply n ⋅ r = nᵢrⱼGᵢRⱼ = 2πnᵢrᵢ; the factors of 2π do not matter for the 0-relation, + # so we can represent the action of n ⋅ r simply by a "matrix" nᵀ; the corresponding + # surface basis is simply the (integer) null space of this matrix, i.e., of nᵀ. In turn, + # this integer null-space can be obtained as the last columns of T⁻¹ from the Smith + # normal form of nᵀ = SΛT (the columns associated with the with zero-columns of the + # diagonal matrix Λ) + nᵀ = Matrix(transpose(n)) + F = smith(nᵀ) + if F.SNF ≠ [1] + error(lazy"failed to determine surface basis; n may have non-unity factors (n=$n)") + end + rs = F.Tinv[:, 2:3] # these two columns are a basis for the surface lattice + @assert size(rs,1) == 3 && size(rs,2) == 2 + + return rs +end + +signed_volume(R₁, R₂, R₃) = dot(R₁, cross(R₂, R₃)) # R₁⋅(R₂×R₃) = R₃⋅(R₁×R₂) = R₂⋅(R₃×R₁) +signed_volume(r₁, r₂) = r₁[1] * r₂[2] - r₁[2] * r₂[1] # z⋅(r₁×r₂) w/ r₁,r₂ lifted to ℝ³ + +function rotation_matrix_n_to_z( + n :: AbstractVector{<:Real}, + z :: AbstractVector{<:Real} = (@SVector[0.0,0.0,1.0]) + ) + # purpose: determine a rotation matrix P such that P⁻¹n = z + # strategy: rotate around an axis perpendicular to both `n` and ̂z by an angle θ (the + # angle that `n` makes to ̂z). The rotation matrix for a rotation around an axis `u` at + # an angle is given by Rodrigues' rotation formula + + length(n) == length(z) == 3 || error(lazy"input must be 3D vectors; got n=$n & z=$z") + + n = normalize(SVector{3,Float64}(n)) # ̂n + z = normalize(SVector{3,Float64}(z)) # ̂z + I₃ = SMatrix{3,3,Float64}(I) + n ≈ z && return I₃ # already aligned with ̂z; no rotation needed + + u = normalize(cross(n, z)) # = ̂n × ̂z (rotation axis) + θ = acos(dot(n, z)) + + # Rodrigues' rotation formula (rotate around `u` by θ): P = I + sin(θ)K + (1 - cos(θ))K² + K = SMatrix{3,3,Float64}([0 -u[3] u[2]; u[3] 0 -u[1]; -u[2] u[1] 0]) + P = I₃ + sin(θ) * K + (1 - cos(θ)) * K^2 + + # computed `P` acts on `n` according to `P*n = z`; however, we want `P` to act on + # as `P⁻¹n = z` (i.e., `P` rotates `n` to `z`); so return P⁻¹ = Pᵀ (cf. P ∈ O(3)) + return P' +end + +## --------------------------------------------------------------------------------------- # + +""" + surface_basis(Rs, n; cartesian=true) + +Compute a basis for the surface-cell obtained from terminating a 3D lattice `Rs` over a +surface surface specified by its normal vector `n` (or Miller indices if `cartesian=false`). + +## Output + +The function returns a tuple `(rs³ᴰ, rs′²ᴰ, P)`, whose elements are described below: + +- `rs³ᴰ`: a `DirectBasis{3}`, whose first two basis vectors lie in the plane of the surface, + and whose third vector is (positively) aligned with the surface normal. + All basis vectors correspond to points in the original direct lattice. + +- `rs′²ᴰ`: a `DirectBasis{2}`, whose basis vectors are given in the local coordinate system + of the surface unit cell; effectively, this is ``(x,y)``-components of `rs³ᴰ` after a + rotation that takes `n` to ``\\hat{\\mathbf{z}}``. The first basis vector is aligned with + the ``\\hat{\\mathbf{x}}``-direction of the local coordinate system. + +- `P`: a rotation matrix that takes `rs³ᴰ` to the local `n`-to-``\\hat{\\mathbf{z}}`` + rotated coordinate system of `rs′²ᴰ`. + In particular, defining `rs′³ᴰ = transform.(DirectPoint.(rs³ᴰ), Ref(P))`, the following + holds: + - `rs′³ᴰ[i] ≈ [rs′²ᴰ[i]..., 0]` for `i ∈ 1:2`, + - `rs′³ᴰ[3]` is (positively) aligned with `[0,0,1]`. + To transform the other way, i.e., from surface-local to lattice-global coordinates, simply + use `P⁻¹ = transpose(P)` instead. + +The returned basis is right-handed. + +## Keyword arguments +- `cartesian :: Bool = true`: whether the surface normal `n` is specified in Cartesian + coordinates (`true`) or in the basis of the reciprocal lattice vectors (`false`), i.e., + corresponding to the Cartesian vector `n[1]*Gs[1] + n[2]*Gs[2] + n[3]*Gs` with + `Gs` denoting the Cartesian representation of the reciprocal lattice vectors, i.e., + `Gs = reciprocalbasis(Rs)`. + The latter case (`false`) is a specification of the surface in terms of its Miller + indices: the coordinates of `n` can then equivalently be interpreted as the inverse of the + surface's intercept with each of the axes spanned by `Rs`. + +## Example +```jldoctest +julia> sgnum = 195; # P23, a cubic lattice + +julia> Rs = directbasis(sgnum, 3) +DirectBasis{3} (cubic): + [1.0, 0.0, 0.0] + [0.0, 1.0, 0.0] + [0.0, 0.0, 1.0] + +julia> n = [1,1,1]; # project along 111 direction + +julia> rs³ᴰ, rs′²ᴰ, P = surface_basis(Rs, n; cartesian=true); + +julia> rs³ᴰ # the surface basis in the original coordinate system +DirectBasis{3} (hexagonal): + [1.0, -1.5700924586837747e-16, -0.9999999999999998] + [-0.9999999999999998, 0.9999999999999997, 0.0] + [1.0, 1.0, 1.0 + +julia> rs′²ᴰ # the in-plane surface vectors in a local "surface coordinate system" +DirectBasis{2} (hexagonal): + [1.414213562373095, 0.0] + [-0.7071067811865475, 1.2247448713915887] + +julia> DirectBasis(transform.(DirectPoint.(rs³ᴰ), Ref(P))) # from rs³ᴰ to rs′²ᴰ coordinates +DirectBasis{3} (hexagonal): + [1.414213562373095, -1.1102230246251563e-16, 0.0] + [-0.7071067811865476, 1.2247448713915885, 0.0] + [0.0, -1.1102230246251563e-16, 1.7320508075688772] + +!!! warning + This function will likely be moved from Crystalline to Bravais.jl at some point in the + future. +``` +""" +function surface_basis( + Rs :: DirectBasis{3}, + n :: AbstractVector{<:Real}; + cartesian :: Bool = true + ) + # TODO: Add tests + # TODO: Move this from Crystalline to Bravais (but this is annoying, because we need + # SmithNormalForm.jl, which is not a dependency of Bravais - and not even + # registered, but vendored by Crystalline). Messy :( + + length(n) == 3 || throw(DimensionMismatch(lazy"normal vector must be 3D, got $n")) + + # if `cartesian=true` (default), we assume `n` has been given in Cartesian coordinates; + # otherwise, if `cartesian=false`, `n` is assumed to be given in the basis of + # `reciprocalbasis(Rs)`, i.e., to be a set of Miller indices for the surface + local n_cartesian + if cartesian + n_cartesian = n + n = latticize(n, reciprocalbasis(Rs)) + else + n_cartesian = cartesianize(float.(n), reciprocalbasis(Rs)) + end + + # cast `n` as the smallest integer vector proportional to `n`, so we can do integer math + # from here on out + local n_int + if !(eltype(n) <: Integer) + d = float_gcd(n...; atol=1e-9) + _n = n ./ d + n_int = round.(Int, _n) + if norm(_n - n_int) > 1e-9 + error(LazyString( + "failed to convert `n` to an integer vector with precision 1e-9; ", + " got n=", _n, " vs. n_int=", n_int, + "; if possible, try to give `n` as integer Miller indices")) + end + else + d = gcd(n) + n_int = isone(d) ? convert.(Int, n) : convert.(Int, div.(n, d)) + end + + # find a basis for the surface lattice, using the constraint n ⋅ rⱼ = 0 for each surface + # basis vector rⱼ; the below returns two basis vectors, in matrix form, specified in the + # basis of the direct lattice vectors of the bulk + rs_snf = surface_basis_from_miller_indices_via_snf(n_int) + + _rs = stack(Rs) * rs_snf # equiv. to `cartesianize(r₃_snf, Rs)`, but need different type + r₁ = DirectPoint{3}(_rs[1,1], _rs[2,1], _rs[3,1]) + r₂ = DirectPoint{3}(_rs[1,2], _rs[2,2], _rs[3,2]) + + # rotate the surface unit cell to have a normal along z, so we can embed it in 2D by + # ignoring the third dimension + P_to_z = rotation_matrix_n_to_z(n_cartesian) + r₁′ = transform(r₁, P_to_z) + r₂′ = transform(r₂, P_to_z) + atol = 1e-9 * maximum(norm, Rs) + if abs(r₁′[3]) > atol && abs(r₂′[3]) > atol + error(lazy"rotation to 2D plane failed; projection to z is appreciably nonzero") + end + + # if the basis is already right-handed (with out-of-plane direction =z), we can proceed + # directly but otherwise we swap r₁ and r₂ before proceeding (specifically, we require + # that a positive signed volume, i.e. ̃n ⋅ (r₁ × r₂) > 0, for the final output) + if signed_volume(r₁′, r₂′) < 0 + r₁′, r₂′ = r₂′, r₁′ + end + tmp_rs′_2D = DirectBasis{2}(r₁′[SOneTo(2)], r₂′[SOneTo(2)]) + + # compute the associated 2D Niggli cell of the rotated basis vectors + rs′²ᴰ, _ = nigglibasis(tmp_rs′_2D) + r₁′²ᴰ = DirectPoint{2}(rs′²ᴰ[1]) + r₂′²ᴰ = DirectPoint{2}(rs′²ᴰ[2]) + + # in case we to use the 2D projection later (i.e., `rs′²ᴰ`), it would be nice if + # its first vector was aligned with its local x-axis; one more rotation to ensure this + θ = atan(r₁′²ᴰ[2], r₁′²ᴰ[1]) + P_to_x = @SMatrix [cos(θ) -sin(θ) 0.0; sin(θ) cos(θ) 0.0; 0.0 0.0 1.0] + P_to_x_2D = P_to_x[SOneTo(2), SOneTo(2)] + r₁′²ᴰ, r₂′²ᴰ = transform.((r₁′²ᴰ, r₂′²ᴰ), Ref(P_to_x_2D)) + rs′²ᴰ = DirectBasis{2}(r₁′²ᴰ, r₂′²ᴰ) + + # overall rotation matrix going back from 2D to 3D + P = P_to_z * P_to_x # correct order since P acts inversely + P⁻¹ = transpose(P) # P⁻¹ = Pᵀ + + # rotate back from the 2D projection to 3D + r₁′³ᴰ = DirectPoint{3}(r₁′²ᴰ..., 0.0) + r₂′³ᴰ = DirectPoint{3}(r₂′²ᴰ..., 0.0) + r₁³ᴰ = transform(r₁′³ᴰ, P⁻¹) + r₂³ᴰ = transform(r₂′³ᴰ, P⁻¹) + + # now: find the "last" or "3rd" basis vector r₃, which we want to be aligned with `n`, + # i.e., is normal to the surface; this requires that n×r₃ = 0; to ensure this, we must + # construct a matrix that represents the action of n× on r₃ (accounting for the fact + # that n specified in the basis of the reciprocal lattice vectors, and r₃ is specified + # in the direct basis); we denote this matrix C, s.t. C*r₃ = n×r₃ = 0; we seek a basis + # for the null-space of C and we will use the smith normal form to get it again + _Gs = reciprocalbasis(Rs) ./ (2π) # factors of 2π are irrelevant & troublesome; remove + C = reduce(hcat, (sum(cross(n_int[k]*_Gs[k], Rs[j]) for k in 1:3) for j in 1:3)) # (n×r₃)ᵢ = Cᵢⱼ(r₃)ⱼ + d = float_gcd(C...; atol=1e-9) # convert to nearest integer matrix + norm(d) ≈ 1 || (C = C ./ d) + C_int = round.(Int, C) + isapprox(C_int, C, atol=1e-8) || error("failed to convert to integer matrix") + F = smith(Matrix(C_int)) + if !(length(F.SNF) == 3 && F.SNF[1] ≠ 0 && F.SNF[2] ≠ 0 && F.SNF[3] == 0) + error("unexpected Smith normal form for 3rd surface-cell basis vector calculation") + end + r₃_snf = DirectPoint{3}(F.Tinv[:,3]) # 3rd basis vector for surface-cell, in `Rs`-basis + r₃³ᴰ = cartesianize(r₃_snf, Rs) + if dot(r₃³ᴰ, n_cartesian) < 0 + r₃³ᴰ = -r₃³ᴰ # turn "negative" parallel to "positive" parallel + end + + # check that everything is still good: r₃ is aligned with `n_cartesian` and the basis is + # right-handed + rs³ᴰ = DirectBasis{3}(r₁³ᴰ, r₂³ᴰ, r₃³ᴰ) + if (signed_volume(r₁′²ᴰ, r₂′²ᴰ) ≤ 0 || + signed_volume(r₁³ᴰ, r₂³ᴰ, r₃³ᴰ) ≤ 0 || + norm(cross(r₃³ᴰ, n_cartesian)) > 1e-9 || + dot(r₃³ᴰ, n_cartesian) < 0) + error("something went wrong during the surface basis computation; signed volume \ + is negative or the out-of-plane (third) surface-unit-cell basis vector is \ + not (positively) oriented towards the surface normal vector") + end + + return rs³ᴰ, rs′²ᴰ, P +end + +float_gcd(x,y; atol=1e-9) = abs(y) < atol ? x : float_gcd(y, x % y; atol) +float_gcd(x, y, z...; kws...) = float_gcd(x, float_gcd(y, z...); kws...) + +## +# Example for visualizing the surface Brillouin zone, and a surface BZ path, together +# with the bulk BZ and the bulk path + +# TODO: Add this either to a tutorial/example/docstring thing, or port parts of it to +# Brillouin.jl +#= +function plot_bz_projection(sgnum, Rs, rs³ᴰ, rs′²ᴰ; axis=NamedTuple()) + Gs = reciprocalbasis(Rs) + kp = irrfbz_path(sgnum, conventionalize(Rs, centering(sgnum))) + + gs³ᴰ = reciprocalbasis(rs³ᴰ) + # pick least-symmetric space group number in 2D + csys = crystalsystem(rs′²ᴰ) + sgnum²ᴰ = csys == "hexagonal" ? 13 #= p3 =# : + csys == "square" ? 10 #= p4 =# : + csys == "rectangular" ? 3 #= p1m1 =# : + csys == "oblique" ? 1 #= p1 =# : error() + + kp′²ᴰ = irrfbz_path(sgnum²ᴰ, rs′²ᴰ) + gs′²ᴰ = basis(kp′²ᴰ) + c′²ᴰ = wignerseitz(gs′²ᴰ) + + vol_Gs = Crystalline.volume(Gs) + vol_gs²ᴰ = Crystalline.volume(gs′²ᴰ) + α = vol_Gs/vol_gs²ᴰ / norm(gs³ᴰ[3]) + + kp³ᴰ = KPath{3}( + Dict(klab => SVector{3,Float64}(kv..., 0) for (klab, kv) in kp′²ᴰ.points), # points + kp′²ᴰ.paths, # paths + gs³ᴰ, # basis + Ref(Brillouin.LATTICE) # setting + ) + kp³ᴰ_shift = KPath{3}( + Dict(klab => SVector{3,Float64}(kv..., α/2) for (klab, kv) in kp′²ᴰ.points), # points + kp′²ᴰ.paths, # paths + gs³ᴰ, # basis + Ref(Brillouin.LATTICE) # setting + ) + #kp³ᴰ = latticize(cartesianize(kp³ᴰ), Gs) + #kp³ᴰ_shift = latticize(cartesianize(kp³ᴰ_shift), Gs) + c = wignerseitz(Gs) + + faxp = plot(c; axis) + plot!(kp; linecolor=:red, markercolor=:red) + proj_col = RGBf(.75, .75, .75) + for f = (-1, 1) + _c = Cell{3}([SVector(v²ᴰ..., f*α/2) for v²ᴰ in c′²ᴰ.verts], + c′²ᴰ.faces, gs³ᴰ, Ref(Brillouin.LATTICE)) + plot!(_c; color= f == 1 ? :black : proj_col) + end + plot!(kp³ᴰ_shift; linecolor=:green, markercolor=:green) + for v²ᴰ in c′²ᴰ.verts + lines!(cartesianize.([SVector(v²ᴰ..., -α/2), SVector(v²ᴰ..., α/2)], Ref(gs³ᴰ)); + color=proj_col) + end + + return faxp, kp, kp³ᴰ, α +end + +## +using Crystalline +using Brillouin +#using CairoMakie +#CairoMakie.activate!(type = "svg") +using GLMakie +GLMakie.activate!() + +sgnum = 216 +Rs = primitivize(directbasis(sgnum), centering(sgnum)) +n = [1,1,1] +rs³ᴰ, rs′²ᴰ, P = surface_basis(Rs, n; cartesian=true) +rs³ᴰ + +faxp, kp, kp³ᴰ, α=plot_bz_projection(sgnum, Rs, rs³ᴰ, rs′²ᴰ; + axis=(;elevation=0.4, azimuth=-4.94, perspectiveness=0.15)) +faxp +#save("bz_projection_$(join(n)).pdf", faxp) +=# \ No newline at end of file