From 5c78497323a31507dadf2d713c126431299fc09a Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Tue, 30 Jun 2020 20:39:04 -0600 Subject: [PATCH 01/21] Lattice provided to BZ functions is now the real-space lattice. --- src/IBZ.jl | 139 ++++++++++++++++++++++++------------ test/make_bz.jl | 24 ++++--- test/make_ibz.jl | 17 ++--- test/make_ibz_rand.jl | 8 ++- test/testHelperFunctions.jl | 4 +- 5 files changed, 122 insertions(+), 70 deletions(-) diff --git a/src/IBZ.jl b/src/IBZ.jl index e2b76ad..fff288b 100644 --- a/src/IBZ.jl +++ b/src/IBZ.jl @@ -7,6 +7,7 @@ export make_bz, make_bz_2d export reduce_bz_old export rand_sc,rand_fcc,rand_bcc,rand_hex,rand_rhom_a,rand_rhom_b,rand_st, rand_bct_a,rand_bct_b,rand_so export lattices,lattice_labels,expectedOrder +export make_recip_latvecs @doc """ function projectToPlane(plane,point) @@ -137,7 +138,7 @@ function reflectionReduce(rPlane,ch) end @doc """ -reduce_bz(lat, at, atom_pos) +reduce_bz(lat, at, atom_pos, convention) returns the reduced Brillouin zone @@ -145,11 +146,14 @@ returns the reduced Brillouin zone - `lat:Array{Float64,3,3}`: The lattice vectors - `at:Array{Int64,n}`: The atom types - `atom_pos:Array{Int64,n,3}`: position of the atoms - +- `convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. """ -function reduce_bz(lat, at, atom_pos) +function reduce_bz(lat, at, atom_pos, convention) - bz = make_bz(lat,false) + bz = make_bz(lat,convention,false) verts = collect(points(polyhedron(bz,CDDLib.Library()))) @@ -200,9 +204,9 @@ function reduce_bz(lat, at, atom_pos) end #function -function reduce_bz_2d(lat, spaceGroup) +function reduce_bz_2d(lat, spaceGroup, convention) - bz = make_bz_2d(lat,false) + bz = make_bz_2d(lat,convention,false) verts = collect(points(polyhedron(bz,CDDLib.Library()))) @@ -259,18 +263,20 @@ returns the brillouin zone for the given lattice, by default returns the verticies the of the bz, can also return the Half space representation """ -function make_bz(lat,vertsOrHrep = true) +function make_bz(lat,convention,vertsOrHrep = true) + + rlat = make_recip_latvecs(lat, convention) # Minkowski reduce the lattice. vector_utils = pyimport("phenum.vector_utils") eps=10^-9; - lat = vector_utils._minkowski_reduce_basis(lat',eps)' + rlat = vector_utils._minkowski_reduce_basis(rlat',eps)' #@show vertsOrHrep #enumerate all lattice points 2 on each side lattice_points = collect(Iterators.product(-2:2,-2:2,-2:2)) #make a array of vectors lattice_points = [collect(i) for i in vec(lattice_points)] #calculate the Cartesian position of each lattice point - cart_point = [lat*i for i in lattice_points] + cart_point = [rlat*i for i in lattice_points] #sort based on distance from origin distance = [norm(i) for i in cart_point] cart_point = cart_point[sortperm(distance)] @@ -293,13 +299,14 @@ function make_bz(lat,vertsOrHrep = true) end function make_bz_2d(lat,vertsOrHrep = true) + rlat = make_recip_latvecs(lat, convention) #@show vertsOrHrep #enumarate all lattice points 2 on each side lattice_points = collect(Iterators.product(-2:2,-2:2)) #make a array of vectors lattice_points = [collect(i) for i in vec(lattice_points)] #calculate the Cartesian position of each lattice point - cart_point = [lat*i for i in lattice_points] + cart_point = [rlat*i for i in lattice_points] #sort based on distance from origin distance = [norm(i) for i in cart_point] cart_point = cart_point[sortperm(distance)] @@ -327,7 +334,7 @@ end reduce_bz_old(lat, at, atom_pos) returns the reduced Brillouin zone, uses a more a cutting plane to reduce the bz instead of -halfspaces, used only for refrence +halfspaces, used only for reference @@ -340,9 +347,9 @@ reflection reduce !!!! - `atom_pos:Array{Int64,n,3}`: position of the atoms """ -function reduce_bz_old(lat, at, atom_pos) +function reduce_bz_old(lat, at, atom_pos,convention) - bz = make_bz(lat) + bz = make_bz(lat,convention) bz = chull(bz) obz = bz @@ -411,7 +418,6 @@ function simplify_faces(chull) points = chull.points faces = chull.simplices polygons = get_polygons(faces,points) - @show polygons for i in 1:size(polygons)[1] polygons[i] = detriangulate(points,polygons[i]) end @@ -432,24 +438,20 @@ function detriangulate(points,polygon) #distance total += euclidean(points[index,:],points[next_point,:]) end - #sometimes total is 0? TEMP FIX - if total < min_distance #&& total !=0 - @show min_distance - @show total + if total < min_distance + #@show min_distance + #@show total min_distance = total min_order = order_num if total == 0 - @show points - @show order - @show possible_orderings - @show polygon + #@show points + #@show order end end end return possible_orderings[min_order] end -#TODO sometimes returns an empty array function get_polygons(faces,points) #each element contains the indices of the vertices that make up a polygon polygons = [] @@ -486,10 +488,7 @@ function get_polygons(faces,points) push!(processed_faces,findex2) end end - if size(polygon)[1] != 0 - @show polygon - push!(polygons,polygon) - end + push!(polygons,polygon) end end return polygons @@ -600,7 +599,11 @@ end #1/a^2 > 1/b^2 + 1/c^2 #TODO figure out why we need certain ordering? (a < b < c) -function rand_fco_a(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) +function rand_fco_a() + a = rand(Uniform(.1,5)) + b = rand(Uniform(a,5)) + c = rand(Uniform(b,5)) + if a > b a,b = b,a end @@ -662,31 +665,73 @@ function rand_basecm_a(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand( end +lattice_labels = ["sc" "fcc" "bcc" "hex" "rhoma" "rhomb" "st" "bcta" "bctb" "so" "baseco" "bco" "fcoa" "fcob" "fcoc" "sm" "basecm" "tric"] #labels = ["sc"] sym = pyimport("bzip.symmetry") sc = rand_sc(1) fcc = rand_fcc(1) bcc = rand_bcc(1) -hex = rand_hex(2,3) -rhoma = rand_rhom_a(.5,pi/7) -rhomb = rand_rhom_b(.5,4*pi/7) -st = rand_st(1,2) -bcta = rand_bct_a(2,3) -bctb = rand_bct_b(2,3) -so = rand_so(1,2,3) -baseco = rand_baseco(1,2,3) -bco = rand_bco(1,2,3) -fcoa = rand_fco_a(1,2,3) -fcob = rand_fco_b(1,.5,.25) -fcoc = rand_fco_c(1,2,3) -sm = rand_sm(1,2,3,pi/4) -#not suer if this is correct -basecm = rand_basecm_a(1,2,3,pi/4) +hex = rand_hex() +rhoma = rand_rhom_a() +rhomb = rand_rhom_b() +st = rand_st() +bcta = rand_bct_a() +bctb = rand_bct_b() +so = rand_so() +baseco = rand_baseco() +bco = rand_bco() +fcoa = rand_fco_a() +fcob = rand_fco_b() +fcoc = rand_fco_c() +sm = rand_sm() +basecm = rand_basecm_a() tric = sym.make_lattice_vectors("triclinic",[1, 2, 3],[pi/13, pi/7, pi/5]) lattices = [sc, fcc, bcc, hex, rhoma,rhomb, st, bcta, bctb, so, baseco, bco, fcoa, fcob, fcoc, sm, basecm, tric] expectedOrder = [48, 48, 48, 24, 12, 12, 16, 16, 16, 8, 8, 8, 8, 8, 8, 4, 4, 2] -#lattices = [sc, fcc, bcc, hex, rhoma,rhomb, st, bcta, bctb, so, baseco, bco, fcoa, sm, tric] -#expectedOrder = [48, 48, 48, 24, 12, 12, 16, 16, 16, 8, 8, 8, 8, 4, 2] -lattice_labels = ["sc" "fcc" "bcc" "hex" "rhoma" "rhomb" "st" "bcta" "bctb" "so" "baseco" "bco" "fcoa" "fcob" "fcoc" "sm" "tric"] - end #module + + +@doc """ + make_recip_latvecs(real_latvecs, convention) + +Calculate the reciprocal lattice vectors. + +# Arguments +- `real_latvecs::AbstractArray{<:Real,2}`: the real space lattice vectors or + primitive translation vectors as columns of a 2x2 or 3x3 array. +- `convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. + +# Returns +- `recip_latvecs::Array{<:Real,2}` the reciprocal lattice vectors (reciprocal + primitive translation vectors) as columns of a 2x2 or 3x3 array. + +# Examples +```jldoctest +using IBZ +real_latvecs=[1 0 0; 0 1 0; 0 0 1] +convention="angular" +make_recip_latvecs(real_latvecs,convention) +# output +3×3 Array{Float64,2}: + 6.28319 0.0 0.0 + 0.0 6.28319 0.0 + 0.0 0.0 6.28319 +``` +""" +function make_recip_latvecs(real_latvecs::AbstractArray{<:Real,2}, + convention::String="ordinary")::Array{Float64,2} + if convention == "ordinary" + recip_latvecs = Array(inv(real_latvecs)') + elseif convention == "angular" + recip_latvecs = Array(2π*inv(real_latvecs)') + else + throw(ArgumentError("The allowed conventions are \"ordinary\" and + \"angular\".")) + end + recip_latvecs +end + +end #module diff --git a/test/make_bz.jl b/test/make_bz.jl index d516faa..e59b27b 100644 --- a/test/make_bz.jl +++ b/test/make_bz.jl @@ -4,19 +4,23 @@ lats = lattices labels =lattice_labels @testset "make_bz_test" begin + symmetry = pyimport("phenum.symmetry") i=1 + at = [1] + atom_pos = [0 0 0] + convention = "ordinary" @testset "$l" for l in labels lat = lats[i] - at = [1] - atom_pos = [0 0 0] - make_IBZ = pyimport("bzip.make_IBZ") - ch = make_IBZ.find_bz(lat) - #convert to julia - ch = chull(ch.points) - chjulia = chull(make_bz(lat)) - @test ch.volume ≈ chjulia.volume - @test comparePolygon(ch, chjulia) + rlat = make_recip_latvecs(lat,convention) + bz = make_bz(lat,convention) + ch = chull(bz) + @test ch.volume ≈ det(rlat) + + spaceGroup = symmetry.get_spaceGroup(lats[i]',at,atom_pos)[1] + ops = [spaceGroup[j,:,:] for j=1:size(spaceGroup)[1]] + unfoldedpts=unique([op*bz[j,:] for op in ops, j in 1:size(bz,1)]) + @test all([any([unfoldedpts[j] ≈ bz[k,:] for j in + 1:length(unfoldedpts)]) for k in 1:size(bz,1)]) i+=1 end end - diff --git a/test/make_ibz.jl b/test/make_ibz.jl index 47fe7f6..ffbe7bf 100644 --- a/test/make_ibz.jl +++ b/test/make_ibz.jl @@ -6,15 +6,16 @@ labels =lattice_labels @testset "bz reduce fixed" begin i = 0 @testset "$l" for l in labels + symmetry = pyimport("phenum.symmetry") + at = [1] + atom_pos = [0 0 0] + convention = "ordinary" i+=1 println(labels[i]) lat = lats[i] println(lat) - at = [1] - atom_pos = [0 0 0] - symmetry = pyimport("phenum.symmetry") println("Getting BZ...") - ch = make_bz(lat) + ch = make_bz(lat,convention) #convert to format ch = chull(ch) # origionalPoints = ch.points @@ -22,14 +23,14 @@ labels =lattice_labels println("Getting symmetry...") spaceGroupTranspose = symmetry.get_spaceGroup(transpose(lat),at,atom_pos)[1] #spaceGroup = symmetry.get_spaceGroup(lat,at,atom_pos)[1] - + println("Getting IBZ...") - rch = reduce_bz(lat,at,atom_pos) + rch = reduce_bz(lat,at,atom_pos,convention) #rch_old = reduce_bz_old(lat,at,atom_pos) - @test expectedOrder[i] ≈ size(spaceGroupTranspose)[1] + @test expectedOrder[i] ≈ size(spaceGroupTranspose,1) - @test rch.volume/ch.volume ≈ 1/size(spaceGroupTranspose)[1] + @test ch.volume/rch.volume ≈ size(spaceGroupTranspose,1) println("Testing unflold...") @test comparePolygon(ch, unfold(rch,spaceGroupTranspose)) println("Testing bz mapping...") diff --git a/test/make_ibz_rand.jl b/test/make_ibz_rand.jl index 7bd457f..7672462 100644 --- a/test/make_ibz_rand.jl +++ b/test/make_ibz_rand.jl @@ -4,8 +4,9 @@ lats = lattices labels =lattice_labels @testset "bz reduce rand" begin + convention = "ordinary" lats = Array[] - test_size =20 + test_size =20 for i in 1:test_size push!(lats,rand_sc()) end for i in 1:test_size push!(lats,rand_fcc()) end for i in 1:test_size push!(lats,rand_bcc()) end @@ -21,7 +22,7 @@ labels =lattice_labels atom_pos = [0 0 0] symmetry = pyimport("phenum.symmetry") - ch = make_bz(lat) + ch = make_bz(lat,convention) #convert to format #ch = copy(hcat(ch...)') ch = chull(ch) @@ -31,7 +32,8 @@ labels =lattice_labels spaceGroupTranspose = symmetry.get_spaceGroup(transpose(lat),at,atom_pos)[1] #println(size(spaceGroupTranspose)[1]) #println(size(spaceGroupTranspose)[1]) - rch = reduce_bz(lat,at,atom_pos) + + rch = reduce_bz(lat,at,atom_pos,convention) #rch_old = reduce_bz_old(lat,at,atom_pos) @test rch.volume/ch.volume ≈ 1/size(spaceGroupTranspose)[1] @test comparePolygon(ch, unfold(rch,spaceGroupTranspose)) diff --git a/test/testHelperFunctions.jl b/test/testHelperFunctions.jl index 5b422ab..b35bf45 100644 --- a/test/testHelperFunctions.jl +++ b/test/testHelperFunctions.jl @@ -26,7 +26,7 @@ function BZMaps(BZ,spaceGroup) for m in ops for i in 1:size(BZ.points)[1] #newVert = transpose(transpose(BZ.points[i,:]) * m) - newVert = m*BZ.points[i,:] + newVert = m*BZ.points[i,:] #check if the newVert maps to another point point_maps = false for j in 1:size(BZ.points)[1] @@ -43,7 +43,7 @@ function BZMaps(BZ,spaceGroup) end function comparePolygon(bz1,bz2) num_points_in_2 = 0 - for i in bz1.vertices + for i in bz1.vertices point = bz1.points[i,:] point_in_2 = false for j in bz2.vertices From a903c70e2549ef003ce9ff27e379a752468ce56d Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 1 Jul 2020 10:14:02 -0600 Subject: [PATCH 02/21] Removed function --- src/IBZ.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/IBZ.jl b/src/IBZ.jl index fff288b..d8af138 100644 --- a/src/IBZ.jl +++ b/src/IBZ.jl @@ -1,6 +1,6 @@ module IBZ -using LinearAlgebra,PyCall,QHull,Polyhedra,CDDLib, Combinatorics,Distances, Distributions +using LinearAlgebra,PyCall,QHull,Polyhedra,CDDLib, Combinatorics,Distances,Distributions export reduce_bz,reducePoints export make_bz, make_bz_2d @@ -298,7 +298,7 @@ function make_bz(lat,convention,vertsOrHrep = true) return bz end -function make_bz_2d(lat,vertsOrHrep = true) +function make_bz_2d(lat,convention,vertsOrHrep = true) rlat = make_recip_latvecs(lat, convention) #@show vertsOrHrep #enumarate all lattice points 2 on each side From 7203a844ef7c0cd8ace44f0540ca2ebd70c37bcf Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Fri, 17 Jul 2020 23:01:41 -0600 Subject: [PATCH 03/21] Overhaul of code --- src/IBZ.jl | 753 ++---------------------------------------- src/lattices.jl | 849 ++++++++++++++++++++++++++++++++++++++++++++++++ src/plotting.jl | 406 +++++++++++++++++++++++ src/symmetry.jl | 590 +++++++++++++++++++++++++++++++++ 4 files changed, 1865 insertions(+), 733 deletions(-) create mode 100644 src/lattices.jl create mode 100644 src/plotting.jl create mode 100644 src/symmetry.jl diff --git a/src/IBZ.jl b/src/IBZ.jl index d8af138..4abdd3f 100644 --- a/src/IBZ.jl +++ b/src/IBZ.jl @@ -1,737 +1,24 @@ module IBZ -using LinearAlgebra,PyCall,QHull,Polyhedra,CDDLib, Combinatorics,Distances,Distributions - -export reduce_bz,reducePoints -export make_bz, make_bz_2d -export reduce_bz_old -export rand_sc,rand_fcc,rand_bcc,rand_hex,rand_rhom_a,rand_rhom_b,rand_st, rand_bct_a,rand_bct_b,rand_so -export lattices,lattice_labels,expectedOrder -export make_recip_latvecs - -@doc """ -function projectToPlane(plane,point) - -returns the projection of point onto plane -""" -function projectToPlane(plane,point) - plane = normalize(plane) - t = (-plane[1]*point[1] - plane[2]*point[2] -plane[3]*point[3])/ - (plane[1]^2+plane[2]^2+plane[3]^2) - return [point[1] + t*plane[1] - point[2] + t*plane[2] - point[3] + t*plane[3]] -end - - -@doc """ -getIntercept(p1,p2,n) -returns a point between p1 and p2 if the plane normal to n intersects -return Nothing if it does not -""" -function getIntercept(p1,p2,n,pPlane=[0.0,0.0,0.0]) - l = p1-p2 - if dot(l,n) ≈ 0 - #lines are parallel do nothing - return Nothing - end - d = dot((pPlane-p1),n)/dot(l,n) - pIntersect = d*l + p1 - #check if pIntersect is inbetween p1 and p2 - dP1P2 = norm(p1-p2) - dP1PIntersect = norm(p1-pIntersect) - dP2PIntersect = norm(p2-pIntersect) - if dP1P2 ≈ dP1PIntersect + dP2PIntersect - return pIntersect - end - # doesn't intersect the line segment - return Nothing -end - -@doc """ -returns the ch points plus the points of intersection from the plane -and edges of the convex hull -""" -function getIntersectPoints(plane,ch) - - oldPoints = ch.points - newPoints = Array{Float64,2}(undef,0,3) - - #create vertices at plane intersect - for s ∈ ch.simplices - #test the 3 lines making up each simplex - i1 = getIntercept(oldPoints[s[1],:],oldPoints[s[2],:],plane) - #println(i1) - if i1 != Nothing - newPoints = vcat(newPoints,i1') - end - i2 = getIntercept(oldPoints[s[2],:],oldPoints[s[3],:],plane) - if i2 != Nothing - newPoints = vcat(newPoints,i2') - end - i3 = getIntercept(oldPoints[s[3],:],oldPoints[s[1],:],plane) - if i3 != Nothing - newPoints = vcat(newPoints,i3') - end - end - return newPoints - -end - -@doc """ -function reducePoints(points) -gets rid of duplicate points and points that lie along an already -defined edge. also rounds points near zero to zero -""" -function reducePoints(points) - ch = chull(points) - reducedPoints = Array{Float64,2}(undef,0,3) - for i ∈ ch.vertices - if isapprox(ch.points[i,:],[0;0;0],atol=eps(Float64)) - reducedPoints = vcat(reducedPoints,[0 0 0]) - else - for j ∈ 1:3 - if isapprox(ch.points[i,j],0,atol=eps(Float64)) - ch.points[i,j] = 0 - end - end - reducedPoints = vcat(reducedPoints,ch.points[i,:]') - end - end - return reducedPoints -end - - - -@doc """ -reflectionReduce(rPlane,ch) - -removes all points from the convex hull ch that are on the opposite side of -the plane defined by the vector rPlane. Creates additional verticies at -plane intersections of the simplicies of the convex hull - -""" -function reflectionReduce(rPlane,ch) - #println("rPlane") - #println(rPlane) - - oldPoints = ch.points - - newPoints = getIntersectPoints(rPlane,ch) - #remove all points that have a negative dot product - for i ∈ ch.vertices - p = oldPoints[i,:] - if isapprox(dot(p,rPlane), 0,atol=eps(Float64)) || dot(p,rPlane) > 0 - #println("is above plane") - #println(p) - newPoints = vcat(newPoints,p') - end - end - #println(newPoints) - try - rPoints = reducePoints(newPoints) - return chull(rPoints) - catch - println("attempted to reduced already reduced volume") - return ch - end -end - -@doc """ -reduce_bz(lat, at, atom_pos, convention) - -returns the reduced Brillouin zone - -# Arguments -- `lat:Array{Float64,3,3}`: The lattice vectors -- `at:Array{Int64,n}`: The atom types -- `atom_pos:Array{Int64,n,3}`: position of the atoms -- `convention::String="ordinary"`: the convention used to go between real and - reciprocal space. The two conventions are ordinary (temporal) frequency and - angular frequency. The transformation from real to reciprocal space is - unitary if the convention is ordinary. -""" -function reduce_bz(lat, at, atom_pos, convention) - - bz = make_bz(lat,convention,false) - - verts = collect(points(polyhedron(bz,CDDLib.Library()))) - - obz = copy(hcat(verts...)') - obz = chull(obz) - - symmetry = pyimport("phenum.symmetry") - #transpose to use phenumlib - spaceGroup = symmetry.get_spaceGroup(transpose(lat),at,atom_pos)[1] - ops = [spaceGroup[i,:,:] for i=1:size(spaceGroup)[1]] - #keep track of the indices of the stabilzer operators ( the operators that did not move a given point) - stabilizer_index = [] - - index = 1 - while size(ops)[1] != 1 - for i=1:(size(ops)[1]) - m = ops[i] - - vert = obz.points[obz.vertices[index],:] - - rotated_vert = m*vert - - # check if m is a stabilizer, if it is, mark it to be kept - if (vert ≈ rotated_vert) - append!(stabilizer_index,i) - else - cut_plane = (vert - rotated_vert) - - bz = bz ∩ HalfSpace(cut_plane,0) - end - end - - #if we only have stabilizers move to a new point - if size(stabilizer_index)[1] == size(ops)[1] - index = (index + 1)%size(obz.vertices)[1] - else - index = 1 - end - #remove all but the stabilizers - ops = [ops[i,:,:][1] for i in stabilizer_index] - stabilizer_index= [] - end - - #convert to verts - verts = collect(points(polyhedron(bz,CDDLib.Library()))) - # convert to 2d array - return chull(copy(hcat(verts...)')) - -end #function - -function reduce_bz_2d(lat, spaceGroup, convention) - - bz = make_bz_2d(lat,convention,false) - - verts = collect(points(polyhedron(bz,CDDLib.Library()))) - - obz = copy(hcat(verts...)') - obz = chull(obz) - - #symmetry = pyimport("phenum.symmetry") - #transpose to use phenumlib - #spaceGroup = symmetry.get_spaceGroup(transpose(lat),at,atom_pos)[1] - #ops = [spaceGroup[i,:,:] for i=1:size(spaceGroup)[1]] - #keep track of the indices of the stabilzer operators ( the operators that did not move a given point) - stabilizer_index = [] - - index = 1 - while size(ops)[1] != 1 - for i=1:(size(ops)[1]) - m = ops[i] - - vert = obz.points[obz.vertices[index],:] - - rotated_vert = m*vert - - # check if m is a stabilizer, if it is, mark it to be kept - if (vert ≈ rotated_vert) - append!(stabilizer_index,i) - else - cut_plane = (vert - rotated_vert) - - bz = bz ∩ HalfSpace(cut_plane,0) - end - end - - #if we only have stabilizers move to a new point - if size(stabilizer_index)[1] == size(ops)[1] - index = (index + 1)%size(obz.vertices)[1] - else - index = 1 - end - #remove all but the stabilizers - ops = [ops[i,:,:][1] for i in stabilizer_index] - stabilizer_index= [] - end - - #convert to verts - verts = collect(points(polyhedron(bz,CDDLib.Library()))) - # convert to 2d array - return chull(copy(hcat(verts...)')) - -end #function - -@doc """ -returns the brillouin zone for the given lattice, - -by default returns the verticies the of the bz, -can also return the Half space representation -""" -function make_bz(lat,convention,vertsOrHrep = true) - - rlat = make_recip_latvecs(lat, convention) - # Minkowski reduce the lattice. - vector_utils = pyimport("phenum.vector_utils") - eps=10^-9; - rlat = vector_utils._minkowski_reduce_basis(rlat',eps)' - #@show vertsOrHrep - #enumerate all lattice points 2 on each side - lattice_points = collect(Iterators.product(-2:2,-2:2,-2:2)) - #make a array of vectors - lattice_points = [collect(i) for i in vec(lattice_points)] - #calculate the Cartesian position of each lattice point - cart_point = [rlat*i for i in lattice_points] - #sort based on distance from origin - distance = [norm(i) for i in cart_point] - cart_point = cart_point[sortperm(distance)] - #create the first half space - p1 =cart_point[1]/2 - d = p1[1]^2+p1[2]^2+p1[3]^2 - bz = HalfSpace(p1,d) - for lp in cart_point[2:end] - p = lp/2 - d = p[1]^2+p[2]^2+p[3]^2 - bz = bz ∩ HalfSpace(p,d) - end - if vertsOrHrep - #return polyhedron(bz,CDDLib.Library()) - verts = collect(points(polyhedron(bz,CDDLib.Library()))) - #2d array - bz = copy(hcat(verts...)') - end - return bz -end - -function make_bz_2d(lat,convention,vertsOrHrep = true) - rlat = make_recip_latvecs(lat, convention) - #@show vertsOrHrep - #enumarate all lattice points 2 on each side - lattice_points = collect(Iterators.product(-2:2,-2:2)) - #make a array of vectors - lattice_points = [collect(i) for i in vec(lattice_points)] - #calculate the Cartesian position of each lattice point - cart_point = [rlat*i for i in lattice_points] - #sort based on distance from origin - distance = [norm(i) for i in cart_point] - cart_point = cart_point[sortperm(distance)] - #create the first half space - p1 =cart_point[1]/2 - d = p1[1]^2+p1[2]^2 - bz = HalfSpace(p1,d) - for lp in cart_point[2:end] - p = lp/2 - d = p[1]^2+p[2]^2 - bz = bz ∩ HalfSpace(p,d) - end - if vertsOrHrep - #return polyhedron(bz,CDDLib.Library()) - verts = collect(points(polyhedron(bz,CDDLib.Library()))) - #2d array - bz = copy(hcat(verts...)') - order = detriangulate(bz,collect(1:size(bz)[1])) - print(order,:) - end - return bz[order,:] -end - -@doc """ -reduce_bz_old(lat, at, atom_pos) - -returns the reduced Brillouin zone, uses a more a cutting plane to reduce the bz instead of -halfspaces, used only for reference - - - -!!! this method fails some test cases, I think it is some sort of finite preciesion error in -reflection reduce !!!! - -# Arguments -- `lat:Array{Float64,3,3}`: The lattice vectors -- `at:Array{Int64,n}`: The atom types -- `atom_pos:Array{Int64,n,3}`: position of the atoms - -""" -function reduce_bz_old(lat, at, atom_pos,convention) - - bz = make_bz(lat,convention) - bz = chull(bz) - obz = bz - - symmetry = pyimport("phenum.symmetry") - #transpose to use phenumlib - spaceGroup = symmetry.get_spaceGroup(transpose(lat),at,atom_pos)[1] - ops = [spaceGroup[i,:,:] for i=1:size(spaceGroup)[1]] - stabilizer_index = [] - - index = 1 - while size(ops)[1] != 1 - for i=1:(size(ops)[1]) - - m = ops[i] - - #create point, rotated point pair to determine cutting plane - vert = obz.points[obz.vertices[index],:] - rotated_vert = m*vert - - # check if m is a stabilizer, if it is, mark it to be kept - if (vert ≈ rotated_vert) - append!(stabilizer_index,i) - else - cut_plane = (vert - rotated_vert) - bz = reflectionReduce(cut_plane,bz) - end - end - #if we only have stabilzers, move to a new point - if size(stabilizer_index)[1] == size(ops)[1] - index = (index + 1)%size(obz.vertices)[1] - else - index = 1 - end - - #remove all but the stabilizers - ops = [ops[i,:,:][1] for i in stabilizer_index] - stabilizer_index= [] - end - - return bz - -end #function - -function write_obj(verts,faces,path) - io = open(path,"w") - for v in eachrow(verts) - write(io, "v ") - for p in v - write(io,string(p, " " )) - end - write(io,"\n") - end - for f in faces - write(io, "f ") - for i in f - write(io,string(i, " " )) - end - write(io,"\n") - end - - close(io) -end - -function simplify_faces(chull) - vertices = chull.vertices - points = chull.points - faces = chull.simplices - polygons = get_polygons(faces,points) - for i in 1:size(polygons)[1] - polygons[i] = detriangulate(points,polygons[i]) - end - return polygons -end -function detriangulate(points,polygon) - indices = unique!(vcat(polygon...)) - possible_orderings = collect(permutations(indices)) - siz = size(possible_orderings[1])[1] - min_distance = Inf - min_order = 0 - for (order_num,order) in enumerate(possible_orderings) - #@show order - total = 0 - for (i,index) in enumerate(order) - #wrap to the beginning - next_point = i == siz ? order[1] : order[i+1] - #distance - total += euclidean(points[index,:],points[next_point,:]) - end - if total < min_distance - #@show min_distance - #@show total - min_distance = total - min_order = order_num - if total == 0 - #@show points - #@show order - end - end - - end - return possible_orderings[min_order] -end -function get_polygons(faces,points) - #each element contains the indices of the vertices that make up a polygon - polygons = [] - processed_faces = [] - for (findex, face) in enumerate(faces) - #@show findex - #@show processed_faces - #check if the face has already been processed - if indexin(findex,processed_faces)[1] === nothing - push!(processed_faces,findex) - polygon = [] - p0 = points[face[1],:] - v1 = points[face[1],:] - points[face[2],:] - v2 = points[face[1],:] - points[face[3],:] - normal = normalize(cross(v1,v2)) - #equation of a plane ax+by+cz = d - a = normal[1] - b = normal[2] - c = normal[3] - d = a *p0[1] + b*p0[2] +c*p0[3] - #find all faces whose points are all in this plane - for (findex2,face2) in enumerate(faces) - in_plane = true - #@show face2 - for i in face2 - p = points[i,:] - if !(a*p[1]+b*p[2]+c*p[3] ≈ d) - in_plane=false - break - end - end - if in_plane - push!(polygon,face2) - push!(processed_faces,findex2) - end - end - push!(polygons,polygon) - end - end - return polygons -end -#lattices used for testing -function order_params(a,b,c) - if a > b - a,b = b,a - end - if b > c - b,c= c,b - end - if a > b - a,b = b,a - end - return a,b,c -end - -function rand_sc(a=rand(Uniform(.1,5))) - return [a 0 0; - 0 a 0; - 0 0 a] -end -function rand_fcc(a=rand(Uniform(.1,5))) - return [0 a/2 a/2; - a/2 0 a/2; - a/2 a/2 0] -end -function rand_bcc(a=rand(Uniform(.1,5))) - return [-a/2 a/2 a/2; - a/2 -a/2 a/2; - a/2 a/2 -a/2] -end -function rand_hex(a=rand(Uniform(.1,5)),c=rand(Uniform(.1,5))) - #I don't know why this is necessary, but phenum lib hangs without this. - #TODO fix - if a < c - temp = a - a = c - c = temp - end - return transpose([a/2 -(a*sqrt(3))/2 0; - a/2 (a*sqrt(3))/2 0; - 0 0 c]) -end -function rand_rhom_a(a = rand(Uniform(.1,3)),α = rand(Uniform(π/10,π/2))) - return transpose([a*cos(α/2) -a*sin(α/2) 0; - a*cos(α/2) a*sin(α/2) 0; - a*cos(α)/cos(α/2) 0 a*sqrt(1-(cos(α)^2)/(cos(α/2)^2))]) -end -function rand_rhom_b(a = rand(Uniform(.1,5)),α = rand(Uniform(π/2,2))) - # goes to 2 so that I don't have a negative root - return transpose([a*cos(α/2) -a*sin(α/2) 0; - a*cos(α/2) a*sin(α/2) 0; - a*cos(α)/cos(α/2) 0 a*sqrt(1-(cos(α)^2)/(cos(α/2)^2))]) -end -function rand_st(a = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - return transpose([a 0 0; - 0 a 0; - 0 0 c]) -end - -function rand_bct_a(a = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - if c > a - #swap - a,c = c,a - end - return transpose([-a/2 a/2 c/2; - a/2 -a/2 c/2; - a/2 a/2 -c/2]) -end - -function rand_bct_b(a = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - if c < a - #swap - a,c = c,a - end - return transpose([-a/2 a/2 c/2; - a/2 -a/2 c/2; - a/2 a/2 -c/2]) -end - -function rand_so(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - a,b,c = order_params(a,b,c) - if a > b || b > c - print("error in so conditions") - end - return transpose([a 0 0; - 0 b 0; - 0 0 c]) -end - -function rand_baseco(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - a,b,c = order_params(a,b,c) - if a > b || b > c - print("error in caseco so conditions") - end - return transpose([a/2 -b/2 0; - a/2 b/2 0; - 0 0 c]) -end - -function rand_bco(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - return transpose([-a/2 b/2 c/2; - a/2 -b/2 c/2; - a/2 b/2 -c/2]) -end - -#1/a^2 > 1/b^2 + 1/c^2 -#TODO figure out why we need certain ordering? (a < b < c) -function rand_fco_a() - a = rand(Uniform(.1,5)) - b = rand(Uniform(a,5)) - c = rand(Uniform(b,5)) - - if a > b - a,b = b,a - end - if 1/a^2 <= 1/b^2 + 1/c^2 - range = sqrt((1/a^2 - 1/b^2)^(-1)) - c = rand(Uniform(range,2*range)) - if 1/a^2 <= 1/b^2 + 1/c^2 - print("incorrect conditions for rand_fco_a") - end - end - return transpose([0 b/2 c/2; - a/2 0 c/2; - a/2 b/2 0]) -end - -#1/a^2 < 1/b^2 + 1/c^2 -function rand_fco_b(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - if a > b - a,b = b,a - end - if 1/a^2 >= 1/b^2 + 1/c^2 - range = sqrt((1/a^2 - 1/b^2)^(-1)) - c = rand(Uniform(0,range)) - if 1/a^2 <= 1/b^2 + 1/c^2 - print("incorrect conditions for rand_fco_a") - end - end - return transpose([0 b/2 c/2; - a/2 0 c/2; - a/2 b/2 0]) -end - -#1/a^2 = 1/b^2 + 1/c^2 -function rand_fco_c(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5))) - if a > b - a,b = b,a - end - if 1/a^2 != 1/b^2 + 1/c^2 - c = sqrt((1/a^2 - 1/b^2)^(-1)) - if 1/a^2 != 1/b^2 + 1/c^2 - print("incorrect conditions for rand_fco_c") - end - end - return transpose([0 b/2 c/2; - a/2 0 c/2; - a/2 b/2 0]) -end - -function rand_sm(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5)),α = rand(Uniform(0,π/2))) - return transpose([a 0 0; - 0 b 0; - 0 c*cos(α) c*sin(α)]) -end -#need help here, what is k_gamma? -function rand_basecm_a(a = rand(Uniform(.1,5)),b = rand(Uniform(.1,5)),c = rand(Uniform(.1,5)),α = rand(Uniform(0,π/2))) - return transpose([a/2 b/2 0; - -a/2 b/2 0; - 0 c*cos(α) c*sin(α)]) -end - - -lattice_labels = ["sc" "fcc" "bcc" "hex" "rhoma" "rhomb" "st" "bcta" "bctb" "so" "baseco" "bco" "fcoa" "fcob" "fcoc" "sm" "basecm" "tric"] -#labels = ["sc"] -sym = pyimport("bzip.symmetry") -sc = rand_sc(1) -fcc = rand_fcc(1) -bcc = rand_bcc(1) -hex = rand_hex() -rhoma = rand_rhom_a() -rhomb = rand_rhom_b() -st = rand_st() -bcta = rand_bct_a() -bctb = rand_bct_b() -so = rand_so() -baseco = rand_baseco() -bco = rand_bco() -fcoa = rand_fco_a() -fcob = rand_fco_b() -fcoc = rand_fco_c() -sm = rand_sm() -basecm = rand_basecm_a() -tric = sym.make_lattice_vectors("triclinic",[1, 2, 3],[pi/13, pi/7, pi/5]) - -lattices = [sc, fcc, bcc, hex, rhoma,rhomb, st, bcta, bctb, so, baseco, bco, fcoa, fcob, fcoc, sm, basecm, tric] -expectedOrder = [48, 48, 48, 24, 12, 12, 16, 16, 16, 8, 8, 8, 8, 8, 8, 4, 4, 2] - - -@doc """ - make_recip_latvecs(real_latvecs, convention) - -Calculate the reciprocal lattice vectors. - -# Arguments -- `real_latvecs::AbstractArray{<:Real,2}`: the real space lattice vectors or - primitive translation vectors as columns of a 2x2 or 3x3 array. -- `convention::String="ordinary"`: the convention used to go between real and - reciprocal space. The two conventions are ordinary (temporal) frequency and - angular frequency. The transformation from real to reciprocal space is - unitary if the convention is ordinary. - -# Returns -- `recip_latvecs::Array{<:Real,2}` the reciprocal lattice vectors (reciprocal - primitive translation vectors) as columns of a 2x2 or 3x3 array. - -# Examples -```jldoctest -using IBZ -real_latvecs=[1 0 0; 0 1 0; 0 0 1] -convention="angular" -make_recip_latvecs(real_latvecs,convention) -# output -3×3 Array{Float64,2}: - 6.28319 0.0 0.0 - 0.0 6.28319 0.0 - 0.0 0.0 6.28319 -``` -""" -function make_recip_latvecs(real_latvecs::AbstractArray{<:Real,2}, - convention::String="ordinary")::Array{Float64,2} - if convention == "ordinary" - recip_latvecs = Array(inv(real_latvecs)') - elseif convention == "angular" - recip_latvecs = Array(2π*inv(real_latvecs)') - else - throw(ArgumentError("The allowed conventions are \"ordinary\" and - \"angular\".")) - end - recip_latvecs -end +include("lattices.jl") +include("plotting.jl") +include("symmetry.jl") +using .lattices, .plotting, .symmetry + +# lattices +export get_recip_latvecs, minkowski_reduce, check_reduced +# 2D Bravais lattices +export genlat_SQR, genlat_HXG, genlat_REC, genlat_RECI, genlat_OBL +# 3D Bravais lattices +export genlat_CUB, genlat_FCC, genlat_BCC, genlat_TET, genlat_BCT, genlat_ORC, + genlat_ORCF, genlat_ORCI, genlat_ORCC, genlat_HEX, genlat_RHL, genlat_MCL, + genlat_MCLC, genlat_TRI + +# plotting +export plot_cells + +# symmetry +export mapto_unitcell, calc_pointgroup, calc_spacegroup, sampleCircle, + sampleSphere, calc_bz, calc_ibz end #module diff --git a/src/lattices.jl b/src/lattices.jl new file mode 100644 index 0000000..a89df67 --- /dev/null +++ b/src/lattices.jl @@ -0,0 +1,849 @@ +module lattices + +using LinearAlgebra, Combinatorics + +export get_recip_latvecs, minkowski_reduce, check_reduced +# 2D Bravais lattices +export genlat_SQR, genlat_HXG, genlat_REC, genlat_RECI, genlat_OBL +# 3D Bravais lattices +export genlat_CUB, genlat_FCC, genlat_BCC, genlat_TET, genlat_BCT, genlat_ORC, + genlat_ORCF, genlat_ORCI, genlat_ORCC, genlat_HEX, genlat_RHL, genlat_MCL, + genlat_MCLC, genlat_TRI + +@doc """ + get_recip_latvecs(real_latvecs, convention) + +Calculate the reciprocal lattice vectors. + +# Arguments +- `real_latvecs::AbstractArray{<:Real,2}`: the real-space lattice vectors or + primitive translation vectors as columns of a 2x2 or 3x3 array. +- `convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. + +# Returns +- `recip_latvecs::Array{<:Real,2}` the reciprocal lattice vectors (reciprocal + primitive translation vectors) as columns of a 2x2 or 3x3 array. + +# Examples +```jldoctest +using IBZ +real_latvecs=[1 0 0; 0 1 0; 0 0 1] +convention="angular" +get_recip_latvecs(real_latvecs,convention) +# output +3×3 Array{Float64,2}: + 6.28319 0.0 0.0 + 0.0 6.28319 0.0 + 0.0 0.0 6.28319 +``` +""" +function get_recip_latvecs(real_latvecs::AbstractArray{<:Real,2}, + convention::String="ordinary")::Array{Float64,2} + if convention == "ordinary" + recip_latvecs = Array(inv(real_latvecs)') + elseif convention == "angular" + recip_latvecs = Array(2π*inv(real_latvecs)') + else + throw(ArgumentError("The allowed conventions are \"ordinary\" and + \"angular\".")) + end + recip_latvecs +end + + +""" + reduce_basis!(basis,k) + +Reduces the `k`th lattice vector. This is accomplished by locating the +lattice point closest to the projection of the `k`th lattice vector onto +the line or plane given by the other lattice vectors, subtracting the +closest lattice point from the `k`th lattice vector, and reordering the +lattice vectors by increasing Euclidean norms. + +# Arguments +-`basis::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. +-`k::Int`: Keeps track of which lattice vector needs to be reduced. + +# Returns +-`basis::AbstractArray{<:Real,2}`: the partially reduced lattice basis as + columns of an array. +-`k::Int`: The index of the lattice vector that needs to be reduced next. + +# Examples +```jldoctest +using IBZ +basis = Array([1 2 0; 0 1 0; 3 2 1]') +k=2 +k=reduce_basis!(basis,k) +basis +# output +3×3 Array{Int64,2}: + 0 1 3 + 1 2 2 + 0 0 1 +""" +function reduce_basis!(basis::AbstractArray{<:Real,2},k::Int)::Int + if k == 2 + v1,v2=[basis[:,i] for i=1:k] + i = round(Int64,dot(v1,v2)/dot(v1,v1)) + vecs = [v2-j*v1 for j=i-1:i+1] + v2 = vecs[sortperm(norm.(vecs))[1]] + if norm(v1) <= norm(v2) + basis[:,k] = v2 + k=3 + else + k=2 + basis[:,2] = v1 + basis[:,1] = v2 + end + + elseif k==3 + v1,v2,v3=[basis[:,i] for i=1:k] + + i = round(Int64,dot(v3,v1)/dot(v1,v1)) + j = round(Int64,dot(v3,v2)/dot(v2,v2)) + + vecs = [[v3-m*v2-l*v1 for m=j-1:j+1,l=i-1:i+1]...] + v3 = vecs[sortperm(norm.(vecs))[1]] + if norm(v2) <= norm(v3) + basis[:,k] = v3 + k=4 + else + if norm(v1) <= norm(v3) + k=3 + basis[:,2] = v3 + basis[:,3] = v2 + else + k=2 + basis[:,1] = v3 + basis[:,2] = v1 + basis[:,3] = v2 + end + end + else + throw(ArgumentError("Argument `k` can be 2 or 3.")) + end + k +end + +@doc """ + minkowski_reduce(basis) + +Minkowski reduce a lattice basis. Follows the logic of Fig. 4 in +\"Low-Dimensional Lattice Basis Reduction Revisited\" by Nguyen, 2009. + +# Arguments +- `basis::AbstractArray{<:Real,2}`: the lattice basis given by the columns + of a 2x2 or 3x3 array. + +# Returns +- `rbasis:: the Minkowski reduced lattice basis as columns of an array. + +# Examples +```jldoctest +using IBZ +basis = [1 2 0; 0 1 0; 0 0 1] +minkowski_reduce(basis) +# output +3×3 Array{Int64,2}: + 0 1 0 + 0 0 1 + 1 0 0 +""" +function minkowski_reduce( + basis::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} + + # Sort the lattice vectors by increasing norm. + order = sortperm([basis[:,i] for i=1:size(basis,1)]) + rbasis = basis[:,order] + + k=2 + while k <= size(rbasis,1) + k=reduce_basis!(rbasis,k) + end + rbasis +end + +@doc """ + check_reduced(basis) + +Verify a lattice basis is Minkowski reduced + +# Arguments +- `basis::AbstractArray{<:Real,2}`: the lattice basis given by the columns + of a 2x2 or 3x3 matrix. + +# Returns +- `Bool`: a boolean that indicates if the lattice basis is reduced. + +# Examples +``` jldoctest +using IBZ +basis = [1 0; 0 1] +check_reduced(basis) +# output +true +``` +""" +function check_reduced(basis::AbstractArray{<:Real,2})::Bool + if size(basis,2) == 2 + (a,b) = [basis[:,i] for i=1:2] + all([norm(a) <= norm(b), + norm(b) <= norm(b+a), + norm(b) <= norm(b-a)]) + elseif size(basis,1) == 3 + (a,b,c) = [basis[:,i] for i=1:3] + all([norm(a) <= norm(b), + norm(b) <= norm(b+a), + norm(b) <= norm(b-a), + norm(b) <= norm(c), + norm(c) <= norm(c+a), + norm(c) <= norm(c-a), + norm(c) <= norm(c+b), + norm(c) <= norm(c-b), + norm(c) <= norm(c+b-a), + norm(c) <= norm(c-b-a), + norm(c) <= norm(c+b+a), + norm(c) <= norm(c-b+a)]) + else + throw(ArgumentError("Can only verify Minkowski reduction in 2D and 3D.")) + end +end + + +@doc """ + genlat_SQR(a) + +Generate a square lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +genlat_SQR(a) +# output +2×2 Array{Int64,2}: + 1 0 + 0 1 +``` +""" +function genlat_SQR(a::Real)::AbstractArray{<:Real,2} + [a 0; 0 a] +end + +@doc """ + genlat_HXG(a) + +Generate a 2D hexagonal lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +genlat_HXG(a) +# output +2×2 Array{Float64,2}: + 1.0 -0.5 + 0.0 0.866025 +``` +""" +function genlat_HXG(a::Real)::AbstractArray{<:Real,2} + Array([a 0; -a/2 a*√3/2]') +end + +@doc """ + genlat_REC(a,b) + +Generate a rectangular lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +genlat_REC(a,b) +# output +2×2 Array{Float64,2}: + 1.0 0.0 + 0.0 1.2 +``` +""" +function genlat_REC(a::Real,b::Real)::AbstractArray{<:Real,2} + Array([a 0; 0 b]') +end + +@doc """ + genlat_RECI(a,b) + +Generate a body-centered rectangular lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +genlat_RECI(a,b) +# output +2×2 Array{Float64,2}: + 0.5 0.5 + -0.6 0.6 +``` +""" +function genlat_RECI(a::Real,b::Real)::AbstractArray{<:Real,2} + if a ≈ b + throw(ArgumentError("The lattice constant must be different sizes for a + rectangular lattice.")) + end + Array([a/2 -b/2; a/2 b/2]') +end + +@doc """ + genlat_OBL(a,b) + +Generate a body-centered rectangular lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +θ=π/3 +test=genlat_OBL(a,b,θ) +# output +2×2 Array{Float64,2}: + 1.0 0.6 + 0.0 1.03923 +``` +""" +function genlat_OBL(a::Real,b::Real,θ::Real)::AbstractArray{<:Real,2} + if a ≈ b + throw(ArgumentError("The lattice constant must be different sizes for an + oblique lattice.")) + end + if θ ≈ π/2 + throw(ArgumentError("The lattice angle must not equal π/2 for an oblique + lattice.")) + end + Array([a 0; b*cos(θ) b*sin(θ)]') +end + +@doc """ + genlat_CUB(a) + +Generate a simple cubic lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +genlat_CUB(a) +# output +3×3 Array{Int64,2}: + 1 0 0 + 0 1 0 + 0 0 1 +``` +""" +function genlat_CUB(a::Real)::AbstractArray{<:Real,2} + [a 0 0; 0 a 0; 0 0 a] +end + +@doc """ + genlat_FCC(a) + +Generate a face-centered cubic lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +genlat_CUB(a) +# output +3×3 Array{Float64,2}: + 0.0 0.5 0.5 + 0.5 0.0 0.5 + 0.5 0.5 0.0 +""" +function genlat_FCC(a::Real)::AbstractArray{<:Real,2} + Array([0 a/2 a/2; a/2 0 a/2; a/2 a/2 0]') +end + +@doc """ + genlat_BCC(a) + +Generate a body-centered cubic lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +genlat_BCC(a) +# output +3×3 Array{Float64,2}: + 0.0 0.5 0.5 + 0.5 0.0 0.5 + 0.5 0.5 0.0 +""" +function genlat_BCC(a::Real)::AbstractArray{<:Real,2} + a/2*Array([-1 1 1; 1 -1 1; 1 1 -1]') +end + +@doc """ + genlat_TET(a) + +Generate a simple tetragonal lattice. + +# Arguments +- `a::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +c=1.2; +genlat_TET(a,c) +# output +3×3 Array{Float64,2}: + 1.0 0.0 0.0 + 0.0 1.0 0.0 + 0.0 0.0 1.2 +""" +function genlat_TET(a::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ c + throw(ArgumentError("The lattices constants must be different for a + tetragonal lattice.")) + end + Array([a 0 0; 0 a 0; 0 0 c]') +end + +@doc """ + genlat_BCT(a) + +Generate a body-centered tetragonal lattice. + +# Arguments +- `a::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +c=1.2; +genlat_BCT(a,c) +# output +3×3 Array{Float64,2}: + -0.5 0.5 0.5 + 0.5 -0.5 0.5 + 0.6 0.6 -0.6 +""" +function genlat_BCT(a::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ c + throw(ArgumentError("The lattices constants must be different for a + tetragonal lattice.")) + end + Array([-a/2 a/2 c/2; a/2 -a/2 c/2; a/2 a/2 -c/2]') +end + +@doc """ + genlat_ORC(a,b,c) + +Generate an orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.4; +c=1.2; +genlat_ORC(a,b,c) +# output +3×3 Array{Float64,2}: + 1.0 0.0 0.0 + 0.0 1.2 0.0 + 0.0 0.0 1.4 +""" +function genlat_ORC(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([a 0 0; 0 b 0; 0 0 c]') +end + +@doc """ + genlat_ORCF(a,b,c) + +Generate a face-centered orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.4; +c=1.2; +genlat_ORCF(a,b,c) +# output +3×3 Array{Float64,2}: + 0.0 0.5 0.5 + 0.6 0.0 0.6 + 0.7 0.7 0.0 +""" +function genlat_ORCF(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([0 b/2 c/2; a/2 0 c/2; a/2 b/2 0]') +end + +@doc """ + genlat_ORCI(a,b,c) + +Generate a body-centered orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.4; +c=1.2; +genlat_ORCI(a,b,c) +# output +3×3 Array{Float64,2}: + -0.5 0.5 0.5 + 0.6 -0.6 0.6 + 0.7 0.7 -0.7 +""" +function genlat_ORCI(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([-a/2 b/2 c/2; a/2 -b/2 c/2; a/2 b/2 -c/2]') +end + +@doc """ + genlat_ORCC(a,b,c) + +Generate a base-centered orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2; +c=1.4; +genlat_ORCC(a,b,c) +# output +3×3 Array{Float64,2}: + 0.5 0.5 0.0 + -0.6 0.6 0.0 + 0.0 0.0 1.4 +""" +function genlat_ORCC(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([a/2 -b/2 0; a/2 b/2 0; 0 0 c]') +end + +@doc """ + genlat_HEX(a,c) + +Generate a hexagonal lattice. + +# Arguments +- `a::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +c=1.2; +genlat_HEX(a,c) +# output +3×3 Array{Float64,2}: + 0.5 0.5 0.0 + -0.866025 0.866025 0.0 + 0.0 0.0 1.2 +""" +function genlat_HEX(a::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ c + throw(ArgumentError("The lattices constants must be different for a + hexagonal lattice.")) + end + Array([a/2 -a*√3/2 0; a/2 a*√3/2 0; 0 0 c]') +end + +@doc """ + genlat_RHL(a,α) + +Generate a rhombohedral lattice. + +# Arguments +- `a::Real`: a lattice constant +- `α::Real`: a lattice angle in radians + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +α=π/6; +genlat_RHL(a,α) +# output +3×3 Array{Float64,2}: + 0.965926 0.965926 0.896575 + -0.258819 0.258819 0.0 + 0.0 0.0 0.442891 +""" +function genlat_RHL(a::Real,α::Real)::AbstractArray{<:Real,2} + if α ≈ π/2 + throw(ArgumentError("The lattice angle cannot be π/2 for a rhombohedral + lattice.")) + end + Array([a*cos(α/2) -a*sin(α/2) 0; a*cos(α/2) a*sin(α/2) 0; + a*cos(α)/cos(α/2) 0 a*√(1-cos(α)^2/cos(α/2)^2)]') +end + +@doc """ + genlat_MCL(a,b,c,α) + +Generate a monoclinic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant +- `α::Real`: a lattice angle in radians less than π/2 + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +c=1.4 +α=π/6; +genlat_MCL(a,b,c,α) +# output +3×3 Array{Float64,2}: + 1.0 0.0 0.0 + 0.0 1.2 1.21244 + 0.0 0.0 0.7 +""" +function genlat_MCL(a::Real,b::Real,c::Real,α::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for a + monoclinic lattice.")) + end + if α ≈ π/2 || α > π/2 + throw(ArgumentError("The lattice angle must be less than π/2.")) + end + Array([a 0 0; 0 b 0; 0 c*cos(α) c*sin(α)]') +end + +@doc """ + genlat_MCLC(a,b,c,α) + +Generate a base-centered monoclinic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant +- `α::Real`: a lattice angle in radians less than π/2 + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +c=1.4 +α=π/6; +genlat_MCLC(a,b,c,α) +# output +3×3 Array{Float64,2}: + 0.5 -0.5 0.0 + 0.6 0.6 1.21244 + 0.0 0.0 0.7 +""" +function genlat_MCLC(a::Real,b::Real,c::Real,α::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for a + monoclinic lattice.")) + end + if α ≈ π/2 || α > π/2 + throw(ArgumentError("The lattice angle must be less than π/2.")) + end + Array([a/2 b/2 0; -a/2 b/2 0; 0 c*cos(α) c*sin(α)]') +end + +@doc """ + genlat_TRI(a,b,c,α,β,γ) + +Generate a triclinic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant +- `α::Real`: a lattice angle in radians less than π/2 + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +c=1.4 +α=π/6; +β=π/3; +γ=π/4; +genlat_TRI(a,b,c,α,β,γ) +# output +3×3 Array{Float64,2}: + 1.0 0.848528 0.7 + 0.0 0.848528 1.01464 + 0.0 0.0 0.663702 +""" +function genlat_TRI(a::Real,b::Real,c::Real,α::Real,β::Real, + γ::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for a + triclinic lattice.")) + end + if α ≈ β || α ≈ γ || β ≈ γ || α ≈ π/2 || β ≈ π/2 || γ ≈ π/2 + throw(ArgumentError("No lattice angles can be the same or equal to π/2 + for a triclinic lattice.")) + end + Array([a 0 0; b*cos(γ) b*sin(γ) 0; + c*cos(β) c/sin(γ)*(cos(α)-cos(β)*cos(γ)) c/sin(γ)* + √(sin(γ)^2-cos(α)^2-cos(β)^2+2*cos(α)*cos(β)*cos(γ))]') +end + +end diff --git a/src/plotting.jl b/src/plotting.jl new file mode 100644 index 0000000..30f5874 --- /dev/null +++ b/src/plotting.jl @@ -0,0 +1,406 @@ +module plotting + +include("symmetry.jl") +include("lattices.jl") + +using .symmetry +using .lattices + +ENV["MPLBACKEND"]="qt5agg" +using PyCall, PyPlot, QHull, LinearAlgebra + +export plot_cells + +@doc """ + function sort_points_comparison(a,b,c) + +A "less than" function for sorting Cartesian points in 2D. + +# Arguments +- `a::AbstractArray{<:Real,1}`: a point in Cartesian coordinates. +- `b::AbstractArray{<:Real,1}`: a point in Cartesian coordinates. +- `c::AbstractArray{<:Real,1}`: the position of the origin in Cartesian + coordinates. + +# Returns +- `Bool`: a boolean that indicates if `a` is less than `b`. + + +# Examples +```jldoctest +using IBZ +a=[0,0] +b=[1,0] +c=[0.5,0.5] +sort_points_comparison(a,b,c) +# output +false +``` +""" +function sort_points_comparison(a::AbstractArray{<:Real,1}, + b::AbstractArray{<:Real,1},c::AbstractArray{<:Real,1})::Bool + (ax,ay)=a + (bx,by)=b + (cx,cy)=c + + if ((ay - cy) > 0) & ((by - cy) < 0) + return true + elseif ((ay - cy) < 0) & ((by - cy) > 0) + return false + elseif ((ay - cy) >= 0) & ((by - cy) >= 0) + return bx > ax + else #((ay - cy) < 0) & ((by - cy) < 0) + return ax > bx + end +end + +@doc """ + affine_trans(pts) + +Calculate the affine transformation that maps the points to the xy-plane. + +# Arguments +- `pts::AbstractArray{<:Real,2}`: an array of Cartesian points as the columns + of a 2D array. The points must all lie on a plane in 3D. + +# Returns +- `M::AbstractArray{<:Real,2}`: the affine transformation matrix that operates + on points in homogeneous coordinates from the left. + +# Examples +```jldoctest +using IBZ +pts = [0.5 0.5 0.5; 0.5 -0.5 0.5; -0.5 0.5 0.5; -0.5 -0.5 0.5]' +affine_trans(pts) +# output +4×4 Array{Float64,2}: + 0.0 -1.0 0.0 0.5 + -1.0 0.0 0.0 0.5 + 0.0 0.0 -1.0 0.5 + 0.0 0.0 0.0 1.0 +""" +function affine_trans(pts::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} + a,b,c = [pts[:,i] for i=1:3] + + # Create a coordinate system with two vectors lying on the plane the points + # lie on. + u = b-a + v = c-a + u = u/norm(u) + v = v - (u⋅v)*u/(u⋅u) + v = v/norm(v) + w = u×v + + # Augmented matrix of affine transform + inv(vcat(hcat([u v w],a),[0 0 0 1])) +end + +@doc """ + function mapto_xyplane(pts) + +Map Cartesian points embedded in 3D to the xy-plane embedded in 2D. + +# Arguments +- `pts::AbstractArray{<:Real,2}`: Cartesian points in 3D as columns of a 2D + array. + +# Returns +- `AbstractArray{<:Real,2}`: Cartesian points in2D as columns of a 2D array. + +# Examples +```jldoctest +using IBZ +pts = [0.5 -0.5 0.5; 0.5 -0.5 -0.5; 0.5 0.5 -0.5; 0.5 0.5 0.5]' +mapto_xyplane(pts) +# output +2×4 Array{Float64,2}: + 0.0 1.0 1.0 0.0 + 0.0 0.0 1.0 1.0 +``` +""" +function mapto_xyplane(pts::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} + + M = affine_trans(pts) + reduce(hcat,[(M*[pts[:,i]..., 1])[1:2] for i=1:size(pts,2)]) +end + +@doc """ + function sortslice_perm(xypts,sxypts) + +Return the permutation vector that maps Cartesian 2D points `xypts` to `sxypts`. + +# Arguments +- `xypts::AbstractArray{<:Real,2}`: 2D Cartesian points as columns of a 2D + array. +- `sxypts::AbstractArray{<:Real,2}`: sorted 2D Cartesian points as columns of a + 2D array. + +# Returns +- `perm::AbstractArray{<:Real,1}`: a permutation vector that sorts an array of + 2D Cartesian coordinates. + +# Examples +```jldoctest +usig IBZ +xypts = [0 0; 0 1; 1 0; 1 1]' +sxypts = [0 1; 1 1; 1 0; 0 0]' +perm=sortslice_perm(xypts,sxypts) +xypts[:,perm] +# output +2×4 Array{Int64,2}: + 0 1 1 0 + 1 1 0 0 +``` +""" +function sortslice_perm(xypts::AbstractArray{<:Real,2}, + sxypts::AbstractArray{<:Real,2})::AbstractArray{<:Real,1} + perm = zeros(Int64,size(xypts,2)) + for i=1:size(xypts,2) + for j=1:size(xypts,2) + if isapprox(xypts[:,i],sxypts[:,j]) + perm[j]=i + continue + end + end + end + perm +end + +@doc """ + function sortpts_perm(pts) + +Calculate the permutation vector that sorts Cartesian points embedded in 3D that + lie on a plane (counter)clockwise with respect to the average of all points. + +# Arguments +- `pts::AbstractArray{<:Real,2}`: Cartesian points embedded in 3D that all lie + on a plane. The points are columns of a 2D array. + +# Returns +- `perm::AbstractArray{<:Real,1}`: the permutation vector that orders the points + clockwise or counterclockwise. + +# Examples +```jldoctest +using IBZ +pts = [0.5 -0.5 0.5; 0.5 -0.5 -0.5; 0.5 0.5 -0.5; 0.5 0.5 0.5]' +perm=sortpts_perm(pts) +pts[:,perm] +# output +3×4 Array{Float64,2}: + 0.5 0.5 0.5 0.5 + 0.5 0.5 -0.5 -0.5 + 0.5 -0.5 -0.5 0.5 +``` +""" +function sortpts_perm(pts::AbstractArray{<:Real,2}) + xypts=mapto_xyplane(pts) + middlept=[sum(xypts[i,:])/size(pts,2) for i=1:2] + sxypts=sortslices(xypts, lt=(x,y)->sort_points_comparison(x,y,middlept), + dims=2) + perm=sortslice_perm(xypts,sxypts) + return perm +end + +@doc """ + get_uniquefacets(ch) + +Calculate the unique facets of a convex hull object. +""" +function get_uniquefacets(ch) + facets = ch.facets + unique_facets = [] + removed=zeros(Int64,size(facets,1)) + for i=1:size(facets,1) + if removed[i] == 0 + removed[i]=1 + face=ch.simplices[i] + for j=i+1:size(facets,1) + if isapprox(facets[i,:],facets[j,:]) + removed[j]=1 + append!(face,ch.simplices[j]) + end + end + face = unique(reduce(hcat,face)[:]) + # Order the corners of the face either clockwise or + # counterclockwise. + face = face[sortpts_perm(Array(ch.points[face,:]'))] + append!(unique_facets,[face]) + end + end + unique_facets +end + +@doc """ + plot_2Dcell(convexhull, axes, color) + +Plot a 2D convex hull + +# Arguments +-`convexhull::Chull{<:Real}`: a convex hull object. +-`axes::PyObject`: an axes object from matplotlib. +-`color::String`: the face color of the convex hull. + +# Returns +-`axes::PyObject`: updated `axes` that includes a plot of the convex hull. + +# Examples +```jldoctest +using IBZ +real_latvecs = [1 0; 0 1] +convention = "ordinary" +bzformat = "convex hull" +bz = calc_bz(real_latvecs,convention,bzformat) + +fig,axes=subplots(figsize=figaspect(1)*1.5) +axes = plot_2Dcell(bz,axes,"deepskyblue"); +color="deepskyblue" +plot_2Dcell(bz,axes,color) +# output +PyObject +``` +""" +function plot_2Dcell(convexhull::Chull{<:Real}, axes::PyObject, + color::String)::PyObject + + patch=pyimport("matplotlib.patches") + collections=pyimport("matplotlib.collections") + + bzpts=Array(convexhull.points') + middlept=[sum(bzpts[i,:])/size(bzpts,2) for i=1:2] + bzpts=sortslices(bzpts, lt=(x,y)->sort_points_comparison(x,y,middlept), + dims=2) + + (x,y)=[bzpts[i,:] for i=1:2] + axes.fill(x,y, facecolor=color,edgecolor="black",linewidth=3) + axes +end + +@doc """ + plot_3Dcell(convexhull, axes, color) + +Plot a 3D convex hull + +# Arguments +-`convexhull::Chull{<:Real}`: a convex hull object. +-`axes::PyObject`: an axes object from matplotlib. +-`color::String`: the face color of the convex hull. + +# Returns +-`axes::PyObject`: updated `axes` that includes a plot of the convex hull. + +# Examples +```jldoctest +using IBZ +real_latvecs = [1 0 0; 0 1 0; 0 0 1] +convention = "ordinary" +bzformat = "convex hull" +bz = calc_bz(real_latvecs,convention,bzformat) +fig = figure(figsize=figaspect(1)*1.5) +axes = fig.add_subplot(111, projection="3d") +axes = plot_3Dcell(bz,axes,"deepskyblue") +# output +PyObject +``` +""" +function plot_3Dcell(convexhull, axes, color, plotrange=false) + + art3d=pyimport("mpl_toolkits.mplot3d.art3d") + + if plotrange ==false + ϵ=0.1*convexhull.volume + plotrange=[[minimum(convexhull.points[:,i])-ϵ, + maximum(convexhull.points[:,i])+ϵ] for i=1:3] + end + + facesᵢ=get_uniquefacets(convexhull) + edgesᵢ=deepcopy(facesᵢ) + for i=1:length(edgesᵢ) + append!(edgesᵢ[i],edgesᵢ[i][1]) + end + + faces=[convexhull.points[i,:] for i in facesᵢ] + edges=[convexhull.points[i,:] for i in edgesᵢ] + + # Reshape faces and edges + faces=[[faces[j][i,:] for i=1:size(faces[j],1)] for j=1:length(faces)] + edges=[[edges[j][i,:] for i=1:size(edges[j],1)] for j=1:length(edges)] + + p=art3d.Poly3DCollection(faces, alpha=0.2, facecolors=color) + l=art3d.Line3DCollection(edges, colors="black", linewidths=1) + + axes.add_collection3d(p) + axes.add_collection3d(l) + + #ax.view_init(-45,-45) + axes.set_xlim3d(plotrange[1]...) + axes.set_ylim3d(plotrange[2]...) + axes.set_zlim3d(plotrange[3]...) + return axes +end + +@doc """ + plot_cells(real_latvecs,atom_types,atom_pos,coords,convention,rtol,atol) + +Plot the Brillouin and Irreducible Brillouin zone in 2D or 3D. + +# Arguments +-`real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as + columns of an array. +-`atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. +-`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal + structure as columns of an array. +-`coords::String`: indicates the positions of the atoms are in \"lattice\" or + \"Cartesian\" coordinates. +-`convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. +-`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for + floating point comparisons. +-`atol::Real=0.0`: an absolute tolerance for floating point comparisons. + +# Returns +-`(fig,axes)`:the figure and axes Python objects. + +# Examples +```jldoctest +real_latvecs = [1 0; .5 1] +atom_types=[0] +coords = "Cartesian" +atom_pos = Array([0 0]') +convention = "ordinary" +(fig,axes)=plot_cells(real_latvecs,atom_types,atom_pos,coords,convention) +``` +""" +function plot_cells(real_latvecs,atom_types,atom_pos,coords,convention, + rtol::Real=sqrt(eps(float(maximum(real_latvecs)))),atol::Real=0.0) + + bzformat = "convex hull" + bz = calc_bz(real_latvecs,convention,bzformat) + ibzformat = "convex hull" + ibz = calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention, + rtol,atol) + + dim = size(real_latvecs,1) + + if dim == 2 + fig,axes=subplots(figsize=figaspect(1)*1.5) + axes = plot_2Dcell(bz,axes,"deepskyblue") + axes = plot_2Dcell(ibz,axes,"coral") + elseif dim == 3 + fig = figure(figsize=figaspect(1)*1.5) + axes = fig.add_subplot(111, projection="3d") + ϵ=0.1*bz.volume + plotrange=[[minimum(bz.points[:,i])-ϵ, + maximum(bz.points[:,i])+ϵ] for i=1:3] + axes = plot_3Dcell(bz,axes,"blue",plotrange) + axes = plot_3Dcell(ibz,axes,"red",plotrange) + else + throw(ArgumentError("The lattice vectors must be in a 2x2 or 3x3 + array.")) + end + (fig,axes) +end + +end #module diff --git a/src/symmetry.jl b/src/symmetry.jl new file mode 100644 index 0000000..ddc39ec --- /dev/null +++ b/src/symmetry.jl @@ -0,0 +1,590 @@ +module symmetry + +using Distances, LinearAlgebra, CDDLib, QHull, Polyhedra, Combinatorics + +include("lattices.jl") +using .lattices + +export mapto_unitcell, calc_pointgroup, calc_spacegroup, sampleCircle, + sampleSphere, calc_bz, calc_ibz + +@doc """ + shoelace(vertices) + +Calculate the area of a polygon with the shoelace algorithm. + +# Arguments +- `vertices::AbstractArray{<:Real,2}`: the xy-coordinates of the vertices + of the polygon as the columns of a 2D array. + +# Returns +- `<:Real`: the area of the polygon. +""" +function shoelace(vertices) + xs = vertices[1,:] + ys = vertices[2,:] + abs(xs[end]*ys[1] - xs[1]*ys[end] + + sum([xs[i]*ys[i+1]-xs[i+1]*ys[i] for i=1:(size(vertices,2)-1)]))/2 +end + +@doc """ + edgelengths(basis,radius) + +Calculate the edge lengths of a parallelepiped circumscribed by a sphere. + +# Arguments +- `basis::Array{<:Real,2}`: a 2x2 or 3x3 matrix whose columns give the + parallelogram or parallelepiped directions, respectively. +- `radius::Real`: the radius of the sphere. +- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for + floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point + comparisons. + +# Returns +- `[la,lb,lc]::Array{Float64,1}`: a list of parallelepiped lengths. + +# Examples +```jldoctest +using IBZ +basis=Array([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]) +radius=3.0 +edgelengths(basis,radius) +# output +3-element Array{Float64,1}: + 3.0 + 3.0 + 3.0 +""" +function edgelengths(basis::Array{<:Real,2}, radius::Real, + rtol::Real=sqrt(eps(float(radius))), atol::Real=0.0)::Array{Float64,1} + + if radius < 0 | isapprox(radius, 0, rtol=rtol, atol=atol) + error("The radius has to be a positive number.") + end + + if size(basis,1) == 2 + (a,b)=[basis[:,i] for i=1:2] + ax,ay=a + bx,by=b + la=2*abs(radius*sqrt(bx^2+by^2)/(ay*bx-ax*by)) + lb=2*abs(radius*sqrt(ax^2+ay^2)/(ay*bx-ax*by)) + return [la,lb] + + elseif size(basis,1) == 3 + (a,b,c)=[basis[:,i] for i=1:3] + ax,ay,az=a + bx,by,bz=b + cx,cy,cz=c + + la=abs(radius*sqrt((by*cx-bx*cy)^2+(bz*cx-bx*cz)^2+(bz*cy-by*cz)^2)/ + (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) + lb=abs((radius*sqrt((ay*cx-ax*cy)^2+(az*cx-ax*cz)^2+(az*cy-ay*cz)^2))/ + (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) + lc=abs(radius*sqrt((ay*bx-ax*by)^2+(az*bx-ax*bz)^2+(az*by-ay*bz)^2)/ + (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) + return [la,lb,lc] + else + throw(ArgumentError("Basis has to be a 2x2 or 3x3 matrix.")) + end +end + + +@doc """ + sampleCircle(basis,radius,offset,rtol,atol) + +Sample uniformly within a circle centered about a point. + +## Arguments +- `basis::Array{<:Real,2}`: a 2x2 matrix whose columns are the grid generating + vectors. +- `radius::Real`: the radius of the circle. +- `offset::Array{<:Real,1}=[0.,0.]`: the xy-coordinates of the center of the + circle. +- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for + floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point + comparisons. + +## Returns +- `pts::Array{Float64,2}` a matrix whose columns are sample points in Cartesian + coordinates. + +## Examples +```jldoctest +using IBZ +basis=Array([1. 0.; 0. 1.]') +radius=1.0 +offset=[0.,0.] +sampleCircle(basis,radius,offset) +# output +2×5 Array{Float64,2}: + 0.0 -1.0 0.0 1.0 0.0 + -1.0 0.0 0.0 0.0 1.0 +``` +""" +function sampleCircle(basis::AbstractArray{<:Real,2}, radius::Real, + offset::AbstractArray{<:Real,1}=[0.,0.], + rtol::Real=sqrt(eps(float(radius))), atol::Real=0.0)::Array{Float64,2} + + # Put the offset in lattice coordinates and round. + (o1,o2)=round.(inv(basis)*offset) + lens=edgelengths(basis,radius) + n1,n2=round.(lens) .+ 1 + + l=0; + pt=Array{Float64,1}(undef,2) + pts=Array{Float64,2}(undef,2,Int((2*n1+1)*(2*n2+1))); + distances=Array{Float64,1}(undef,size(pts,2)) + for (i,j) in Iterators.product((-n1+o1):(n1+o1),(-n2+o2):(n2+o2)) + l+=1 + pt=BLAS.gemv('N',float(basis),[i,j]) + pts[:,l]=pt + distances[l]=euclidean(pt,offset) + end + + return pts[:,distances.<=radius] + +end + +@doc """ + sampleSphere(basis,radius,offset,rtol,atol) + +Sample uniformly within a circle centered about a point. + +# Arguments +- `basis::Array{<:Real,2}`: a 3x3 matrix whose columns are the grid generating + vectors. +- `radius::Real`: the radius of the sphere. +- `offset::Array{<:Real,1}=[0.,0.]`: the xy-coordinates of the center of the + circle. +- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for + floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point + comparisons. + +# Returns +- `pts::Array{Float64,2}` a matrix whose columns are sample points in Cartesian + coordinates. + +# Examples +```jldoctest +using IBZ +basis=Array([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]) +radius=1.0 +offset=[0.,0.,0.] +sampleSphere(basis,radius,offset) +# output +3×7 Array{Float64,2}: + 0.0 0.0 -1.0 0.0 1.0 0.0 0.0 + 0.0 -1.0 0.0 0.0 0.0 1.0 0.0 + -1.0 0.0 0.0 0.0 0.0 0.0 1.0 +``` +""" +function sampleSphere(basis::AbstractArray{<:Real,2}, radius::Real, + offset::Array{<:Real,1}=[0.,0.,0.], rtol::Real=sqrt(eps(float(radius))), + atol::Real=0.0)::Array{Float64,2} + + # Put the offset in lattice coordinates and round. + (o1,o2,o3)=round.(inv(basis)*offset) + lens=edgelengths(basis,radius) + n1,n2,n3=round.(lens) .+ 1 + + l=0; + pt=Array{Float64,1}(undef,3) + pts=Array{Float64,2}(undef,3,Int((2*n1+1)*(2*n2+1)*(2*n3+1))); + distances=Array{Float64,1}(undef,size(pts,2)) + for (i,j,k) in Iterators.product((-n1+o1):(n1+o1),(-n2+o2):(n2+o2), + (-n3+o3):(n3+o3)) + l+=1 + pt=BLAS.gemv('N',float(basis),[i,j,k]) + pts[:,l]=pt + distances[l]=euclidean(pt,offset) + end + + return pts[:,distances.<=radius] +end + +@doc """ + calc_pointgroup(real_latvecs,rtol,atol) + +Calculate the point group of lattice in 2D or 3D. + +# Arguments +-`real_latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns + of an array. +-`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))`: a relative tolerance for + floating point comparisons. It is used to compare lengths of vectors and the + volumes of primitive cells. +-`atol::Real=0.0`: an absolute tolerance for floating point comparisons. It is + used to compare lengths of vectors and the volumes of primitive cells. + +# Returns +-`pointgroup::Array{Array{Float64,2},1}`: the point group of the lattice. The + operators operate on points in Cartesian coordinates. + +# Examples +```jldoctest +using IBZ +basis = [1 0; 0 1] +calc_pointgroup(basis) +# output +8-element Array{Array{Float64,2},1}: + [0.0 -1.0; -1.0 0.0] + [0.0 1.0; -1.0 0.0] + [-1.0 0.0; 0.0 -1.0] + [-1.0 0.0; 0.0 1.0] + [1.0 0.0; 0.0 -1.0] + [1.0 0.0; 0.0 1.0] + [0.0 -1.0; 1.0 0.0] + [0.0 1.0; 1.0 0.0] +``` +""" +function calc_pointgroup(real_latvecs::AbstractArray{<:Real,2}, + rtol::Real=sqrt(eps(float(maximum(real_latvecs)))), + atol::Real=0.0)::Array{Array{Float64,2},1} + + dim=size(real_latvecs,1) + real_latvecs = minkowski_reduce(real_latvecs) + radius = maximum([norm(real_latvecs[:,i]) for i=1:dim]) + if dim==2 + pts=sampleCircle(real_latvecs,radius,[0,0],rtol,atol) + elseif dim==3 + pts=sampleSphere(real_latvecs,radius,[0,0,0],rtol,atol) + else + throw(ArgumentError("The lattice basis must be a 2x2 or 3x3 array.")) + end + + normsᵢ=[norm(real_latvecs[:,i]) for i=1:dim] + sizeᵢ=abs(det(real_latvecs)) + pointgroup=Array{Array{Float64,2}}([]) + inv_latvecs=inv(real_latvecs) + for perm=permutations(1:size(pts,2),dim) + real_latvecsᵣ=pts[:,perm] + normsᵣ=[norm(real_latvecsᵣ[:,i]) for i=1:dim] + sizeᵣ=abs(det(real_latvecsᵣ)) + # Point operation preserves length of lattice vectors and + # volume on primitive cell. + if (isapprox(normsᵢ,normsᵣ,rtol=rtol,atol=atol) & + isapprox(sizeᵢ,sizeᵣ,rtol=rtol,atol=atol)) + op=inv_latvecs*real_latvecsᵣ + # Point operators are orthogonal. + if isapprox(op*op',I,rtol=rtol,atol=atol) + append!(pointgroup,[op]) + end + end + end + pointgroup +end + + +@doc """ + mapto_unitcell(pt,real_latvecs,inv_latvecs,coords,rtol) + +Map a point to the first unit (primitive) cell. + +# Arguments +-`pt::AbstractArray{<:Real,1}`: a point in lattice or Cartesian coordinates. +-`latvecs::AbstractArray{<:Real,2}`: the basis vectors of the lattice as columns + of an array. +-`inv_latvecs::AbstractArray{<:Real,2}`: the inverse of the matrix of that + contains the lattice vectors. +-`coords::String`: indicates whether `pt` is in \"Cartesian\" or \"lattice\" + coordinates. +-`rtol::Real=sqrt(eps(float(maximum(inv_latvecs))))`: a relative tolerance for + floating point comparisons. Finite precision errors creep up when `pt` is + transformed to lattice coordinates because the transformation requires + calculating a matrix inverse. The components of the point in lattice + coordinates are checked to ensure that values close to 1 are equal to 1. + +# Returns +-`AbstractArray{<:Real,1}`: a translationally equivalent point to `pt` in the + first unit cell in the same coordinates. + +#Examples +```jldoctest +real_latvecs = [0 1 2; 0 -1 1; 1 0 0] +inv_latvecs=inv(real_latvecs) +pt=[1,2,3.2] +coords = "Cartesian" +mapto_unitcell(pt,real_latvecs,inv_latvecs,coords) +# output +3-element Array{Float64,1}: + 0.0 + 0.0 + 0.20000000000000018 +``` +""" +function mapto_unitcell(pt::AbstractArray{<:Real,1}, + latvecs::AbstractArray{<:Real,2},inv_latvecs::AbstractArray{<:Real,2}, + coords::String,rtol::Real=sqrt(eps(float(maximum(inv_latvecs)))) + )::Array{Float64,1} + if coords == "Cartesian" + Array{Float64,1}(latvecs*[isapprox(mod(c,1),1,rtol=rtol) ? 0 : mod(c,1) + for c=inv_latvecs*pt]) + elseif coords == "lattice" + Array{Float64,1}([isapprox(mod(c,1),1,rtol=rtol) ? 0 : mod(c,1) for + c=pt]) + else + throw(ArgumentError("Allowed coordinates of the point are \"Cartesian\" + and \"lattice\".")) + end +end + +@doc """ + calc_spacegroup(real_latvecs,atom_types,atom_pos,coords,rtol,atol) + +Calculate the space group of a crystal structure. + +# Arguments +-`real_latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns + of an array. +-`atom_types::AbstractArray{<:Int,1}`: a list of atom types as integers. +-`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal + structure as columns of an array. +-`coords::String`: indicates the positions of the atoms are in \"lattice\" or + \"Cartesian\" coordinates. +-`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for + floating point comparisons. +-`atol::Real=0.0`: an absolute tolerance for floating point comparisons. + +# Returns +-`spacegroup`: the space group of the crystal structure. The first element of + `spacegroup` is a list of fractional translations, and the second element is + a list of point operators. The translations are in Cartesian coordinates, + and the operators operate on points in Cartesian coordinates. + +# Examples +```jldoctest +real_latvecs = Array([1 0; 2 1]') +atom_types = [0, 1] +atom_pos = Array([0 0; 0.5 0.5]') +coords = "Cartesian" +spacegroup=calc_spacegroup(real_latvecs,atom_types,atom_pos,coords) +# output +(Any[[0.0, 0.0], [0.0, 0.0]], Any[[1.0 0.0; 0.0 1.0], [0.0 1.0; 1.0 0.0]]) +``` +""" +function calc_spacegroup(real_latvecs::AbstractArray{<:Real,2}, + atom_types::AbstractArray{<:Int,1},atom_pos::AbstractArray{<:Real,2}, + coords::String,rtol::Real=sqrt(eps(float(maximum(real_latvecs)))), + atol::Real=0.0) + + if length(atom_types) != size(atom_pos,2) + throw(ArgumentError("The number of atom types and positions must be the + same.")) + end + + real_latvecs = minkowski_reduce(real_latvecs) + inv_latvecs = inv(real_latvecs) + pointgroup = calc_pointgroup(real_latvecs,rtol,atol) + numatoms = length(atom_types) + + # Map points to the primitive cell. + atom_pos = reduce(hcat,[mapto_unitcell(atom_pos[:,i],real_latvecs, + inv_latvecs,coords,rtol) for i=1:numatoms]) + + # Place atom positions in Cartesian coordinates. + if coords == "lattice" + atom_pos = real_latvecs*atom_pos + end + + ops_spacegroup=[] + trans_spacegroup=[] + kindᵣ=atom_types[1] + + # Atoms of the same type as the first. + same_atoms=findall(x->x==kindᵣ,atom_types) + for op in pointgroup + # Use the first atom to calculate allowed fractional translations. + posᵣ=op*atom_pos[:,1] + + for atomᵢ=same_atoms + # Calculate a candidate fractional translation. + ftrans = mapto_unitcell(atom_pos[:,atomᵢ]-posᵣ,real_latvecs, + inv_latvecs,"Cartesian") + + # Check that all atoms land on the lattice when rotated and + # translated. + mapped = false + for atomⱼ = 1:numatoms + kindⱼ = atom_types[atomⱼ] + posⱼ = op*atom_pos[:,atomⱼ] + ftrans + mapped = any([kindⱼ == atom_types[i] && + isapprox(posⱼ,atom_pos[:,i],rtol=rtol,atol=atol) ? + true : false for i=1:numatoms]) + if !mapped continue end + end + if mapped + append!(ops_spacegroup,[op]) + append!(trans_spacegroup,[ftrans]) + end + end + end + (trans_spacegroup,ops_spacegroup) +end + +@doc """ + calc_bz(real_latvecs, convention, vertsOrHrep=true, eps=1e-9) + +Calculate the Brillouin zone for the given real-spcae lattice. + +# Arguments +- `real_latvecs::AbstractArray{<:Real,2}`: the real-space lattice vectors or + primitive translation vectors as columns of a 3x3 array. +- `convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. +- `bzformat::String`: the format of the Brillouin zone. Options include + \"convex-hull\" and \"half-space\". + +# Returns +- `bz`: the vertices or half-space representation of the Brillouin zone + depending on the value of `vertsOrHrep`. + +# Examples +```jldoctest +using IBZ +real_latvecs=[1 0 0; 0 1 0; 0 0 1] +convention="ordinary" +vertsOrHrep=true +make_bz(real_latvecs,convention,vertsOrHrep) +# output +8×3 Array{Float64,2}: + 0.5 -0.5 0.5 + 0.5 -0.5 -0.5 + 0.5 0.5 -0.5 + 0.5 0.5 0.5 + -0.5 0.5 -0.5 + -0.5 0.5 0.5 + -0.5 -0.5 -0.5 + -0.5 -0.5 0.5 +``` +""" +function calc_bz(real_latvecs::AbstractArray{<:Real,2},convention::String, + bzformat::String) + + recip_latvecs = get_recip_latvecs(real_latvecs,convention) + + if size(real_latvecs) == (2,2) + latpts = reduce(hcat,[recip_latvecs*[i,j] for + (i,j)=Iterators.product(-2:2,-2:2)]) + else + latpts = reduce(hcat,[recip_latvecs*[i,j,k] for + (i,j,k)=Iterators.product(-2:2,-2:2,-2:2)]) + end + + distances = [norm(latpts[:,i]) for i=1:size(latpts,2)] .^2 ./2 + bz = HalfSpace(latpts[:,2],distances[2]) + for i=3:size(latpts,2) + bz = bz ∩ HalfSpace(latpts[:,i],distances[i]) + end + verts = reduce(hcat,points(polyhedron(bz,CDDLib.Library()))) + bzvolume = chull(Array(verts')).volume + + if !(bzvolume ≈ abs(det(recip_latvecs))) + error("The area or volume of the Brillouin zone is incorrect.") + end + + if bzformat == "half-space" + bz + elseif bzformat == "convex hull" + chull(Array(verts')) + else + throw(ArgumentError("Formats for the BZ include \"half-space\" and + \"convex hull\".")) + end +end + +@doc """ + calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention,rtol, + atol) + +Calculate the irreducible Brillouin zone of a crystal structure in 2D or 3D. + +# Arguments +-`real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as + columns of an array. +-`atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. +-`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal + structure as columns of an array. +-`coords::String`: indicates the positions of the atoms are in \"lattice\" or + \"Cartesian\" coordinates. +-`ibzformat::String`: the format of the irreducible Brillouin zone. Options + include \"convex-hull\" and \"half-space\". +-`convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. +-`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for + floating point comparisons. +-`atol::Real=0.0`: an absolute tolerance for floating point comparisons. + +# Returns +-`ibz`: the irreducible Brillouin zone as a convex hull or intersection of + half-spaces. + +# Examples +```jldoctest +real_latvecs = [1 0; 0 1] +convention="ordinary" +atom_types=[0] +atom_pos = Array([0 0]') +coords = "Cartesian" +ibzformat = "convex hull" +ibz = calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention) +# output +Convex Hull of 3 points in 2 dimensions +Hull segment vertex indices: +[1, 2, 3] +Points on convex hull in original order: + +[0.0 0.0; 0.5 0.0; 0.5 0.5] +``` +""" +function calc_ibz(real_latvecs::AbstractArray{<:Real,2}, + atom_types::AbstractArray{<:Int,1},atom_pos::AbstractArray{<:Real,2}, + coords::String,ibzformat::String,convention::String="ordinary", + rtol::Real=sqrt(eps(float(maximum(real_latvecs)))), + atol::Real=0.0)::Chull{<:Real} + + pointgroup = calc_spacegroup(real_latvecs,atom_types,atom_pos,coords)[2] + copy_pg = deepcopy(pointgroup) + + bzformat = "half-space" + bz = calc_bz(real_latvecs,convention,bzformat) + bz_vertices = collect(points(polyhedron(bz,CDDLib.Library()))) + ibz = bz + while length(copy_pg) > 0 + op = pop!(copy_pg) + for v=bz_vertices + vʳ=op*v + if vʳ != v + a = vʳ-v + ibz = ibz ∩ HalfSpace(a,0) + break + end + end + end + + ibz_vertices = reduce(hcat,points(polyhedron(ibz,CDDLib.Library()))) + bz_vertices = reduce(hcat,bz_vertices) + bzvolume = chull(Array(bz_vertices')).volume + ibzvolume = chull(Array(ibz_vertices')).volume + + if !(ibzvolume ≈ bzvolume/size(pointgroup,1)) + error("The area or volume of the irreducible Brillouin zone is + incorrect.") + end + if ibzformat == "half-space" + ibz + elseif ibzformat == "convex hull" + chull(Array(ibz_vertices')) + else + throw(ArgumentError("Formats for the IBZ include \"half-space\" and + \"convex hull\".")) + end +end + +end #module From 98a3f23d21b6410276e9964567a757a1e7444668 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 13:15:19 -0600 Subject: [PATCH 04/21] Added unit tests --- Project.toml | 5 + src/IBZ.jl | 10 +- src/Utilities.jl | 43 +++++++++ test/lattices.jl | 76 +++++++++++++++ test/runtests.jl | 7 +- test/symmetry.jl | 245 +++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 378 insertions(+), 8 deletions(-) create mode 100644 src/Utilities.jl create mode 100644 test/lattices.jl create mode 100644 test/symmetry.jl diff --git a/Project.toml b/Project.toml index f1a975e..4540a83 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,12 @@ CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" QHull = "a8468747-bd6f-53ef-9e5c-744dbc5c59e7" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/src/IBZ.jl b/src/IBZ.jl index 4abdd3f..eacf5dc 100644 --- a/src/IBZ.jl +++ b/src/IBZ.jl @@ -1,9 +1,11 @@ module IBZ -include("lattices.jl") -include("plotting.jl") -include("symmetry.jl") -using .lattices, .plotting, .symmetry +include("Lattices.jl") +include("Plotting.jl") +include("Symmetry.jl") +include("Utilities.jl") + +using .Lattices, .Plotting, .Symmetry, .Utilities # lattices export get_recip_latvecs, minkowski_reduce, check_reduced diff --git a/src/Utilities.jl b/src/Utilities.jl new file mode 100644 index 0000000..83563ba --- /dev/null +++ b/src/Utilities.jl @@ -0,0 +1,43 @@ +module Utilities + +export unique + +@doc """ + unique(points,rtol,atol) + +Remove duplicate points from an array. + +# Arguments +-`points::AbstractArray{<:Real,2}`: the points are columns of a 2D array. +-`rtol::Real=sqrt(eps(float(maximum(points))))`: a relative tolerance for floating + point comparisons. +-`atol::Real=0.0`: an absolume tolerance for floating point comparisons. + +# Returns +-`uniquepts::AbstractArray{<:Real,2}`: a 2D array of unique points as columns. + +# Examples +```jldoctest +using IBZ +points=Array([1 2; 2 3; 3 4; 1 2]') +IBZ.Utilities.unique(points) +# output +2×3 Array{Int64,2}: + 1 2 3 + 2 3 4 +``` +""" +function unique(points::AbstractArray{<:Real,2}, + rtol::Real=sqrt(eps(float(maximum(points)))), + atol::Real=0.0)::AbstractArray{<:Real,2} + uniquepts=[] + for i=1:size(points,2) + pt=points[:,i] + if !any([isapprox(pt,uniquepts[i],rtol=rtol) for i=1:length(uniquepts)]) + append!(uniquepts,[pt]) + end + end + reduce(hcat,uniquepts) +end + +end # module diff --git a/test/lattices.jl b/test/lattices.jl new file mode 100644 index 0000000..253a04b --- /dev/null +++ b/test/lattices.jl @@ -0,0 +1,76 @@ +using IBZ, Test + +@testset "lattices" begin + @testset "minkowski_reduce" begin + # Test 1 + basis = [0.997559590093 0.327083693383 0.933708574257; + 0.246898693832 0.853865606330 0.534685117471; + 0.470554239317 0.503012665094 0.122191325919] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 2 + basis = [0.566824707611 0.230580598095 0.644795207463; + 0.062164441073 0.069592500822 0.444446612167; + 0.400346380548 0.334526702097 0.184730262940] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 3 + basis = [0.965915375496 0.913559885031 0.715882418258; + 0.907651021009 0.074522909849 0.841837872220; + 0.273128578139 0.826489413457 0.040093668849] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 4 + basis = [0.747927685471 0.288681880791 0.388321675735; + 0.519739433266 0.818689570053 0.616453247920; + 0.886262688190 0.399796817611 0.970347244023] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 5 + basis = [0.367390897375 0.660775660474 0.211656825515; + 0.084605591875 0.481734488916 0.429855314823; + 0.585439077728 0.797785266051 0.636905287099] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 6 + basis = [0.480636898755 0.799208677212 0.981910531646; + 0.395947716585 0.894180214969 0.229684979406; + 0.466195536761 0.321402105618 0.558040043342] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 7 + basis = [0.073754263634 0.919171496982 0.084716875229; + 0.101486421912 0.827727284628 0.362098268238; + 0.199341503544 0.463651019283 0.902481585522] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 8 + basis = [0.005355229278 0.995622271352 0.278599636291; + 0.047892267806 0.312276051592 0.364509680531; + 0.921399596836 0.679640515605 0.122704568964] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 9 + basis = [0.347081288518 0.164663213375 0.877540564643; + 0.239622809292 0.463104682188 0.998022806449; + 0.154273976230 0.627982335854 0.966889709417] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + # Test 10 + basis = [0.342243816073 0.811478344011 0.567158427035; + 0.643460104906 0.993952140553 0.710611242854; + 0.941907678956 0.015738060731 0.598940997350] + rbasis=minkowski_reduce(basis) + @test check_reduced(rbasis) + + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 849af3d..7bf520e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,4 @@ using IBZ,Test,LinearAlgebra,PyCall,QHull -include("testHelperFunctions.jl") -include("make_bz.jl") -include("make_ibz.jl") -include("make_ibz_rand.jl") + +include("lattices.jl") +include("symmetry.jl") diff --git a/test/symmetry.jl b/test/symmetry.jl new file mode 100644 index 0000000..e2c6ea4 --- /dev/null +++ b/test/symmetry.jl @@ -0,0 +1,245 @@ +using Test,IBZ + +# Lattice vectors +# 2D +a=ℯ +b=π +θ=π/3.2 +sqr_latvecs=genlat_SQR(a) +rec_latvecs=genlat_REC(a,b) +reci_latvecs=genlat_RECI(a,b) +hxg_latvecs=genlat_HXG(a) +obl_latvecs=genlat_OBL(a,b,θ) + +# 3D +a=π +cub_latvecs=genlat_CUB(a) +fcc_latvecs=genlat_FCC(a) +bcc_latvecs=genlat_BCC(a) + +c=ℯ +#1 c < a +#2 c > a +tet_latvecs=genlat_TET(a,c) +bct_latvecs1=genlat_BCT(a,c) +bct_latvecs2=genlat_BCT(c,a) + +b=5/7+(π+ℯ)/3 +orc_latvecs=genlat_ORC(a,b,c) + +# Three face-centered orthorhombic cases +# 1/a² <,=,> 1/b²+1/c² +a=π +b=5/7+(π+ℯ)/3 +c=ℯ +orcf_latvecs1=genlat_ORCF(a,b,c) +a=π/2 +b=5/7+(π+ℯ)/3 +c=ℯ +orcf_latvecs2=genlat_ORCF(a,b,c) +b=5/7 +c=ℯ +a=b*c/√(b^2 + c^2) +orcf_latvecs3=genlat_ORCF(a,b,c) + +# One body-centered orthorhombic case +a=ℯ +b=π +c=4.0+1/3 +orci_latvecs=genlat_ORCI(a,b,c) +# One base-centered orthorhombic case +orcc_latvecs=genlat_ORCC(a,b,c) +# One hexagonal case +hex_latvecs=genlat_HEX(a,c) + +# Two rhombohedral case α < π/2 and α > π/2 +α=π/7+π/2 +rhl_latvecs1=genlat_RHL(a,α) +α=π/7 +rhl_latvecs2=genlat_RHL(a,α) + +# One monoclinic case +mcl_latvecs=genlat_MCL(a,b,c,α) + +# Five base-centered monoclinic cases +# kᵧ > π/2 +a=1+1/3 +b=a+1/3 +c=b+1/3 +α=7*π/20 +mclc_latvecs1=genlat_MCLC(a,b,c,α) +# rlatvecs=get_recip_latvecs(mclc_latvecs1) +# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show γ>π/2 + +# kᵧ == π/2 +a=ℯ +b=π +c=4+1/3 +α=1.0456607472366295 +mclc_latvecs2=genlat_MCLC(a,b,c,α) +# rlatvecs=get_recip_latvecs(mclc_latvecs2) +# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show γ==π/2 + +# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 <1 +a=ℯ +b=π +c=4+1/3 +α=.3 +mclc_latvecs3=genlat_MCLC(a,b,c,α) +#@show b*cos(α)/c + b^2*sin(α)^2/a^2 <1 +# rlatvecs=get_recip_latvecs(mclc_latvecs3) +# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show γ<π/2 + +# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 +a=1.2 +b=1.3 +c=1.4 +α=0.6 +mclc_latvecs4=genlat_MCLC(a,b,c,α) +# @show b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 +# rlatvecs=get_recip_latvecs(mclc_latvecs4) +# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show γ<π/2 + +# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 +b=1.1 +c=1.2 +α=1.1 +a=(b*√(c)*sin(α))/√(c - b*cos(α)); +mclc_latvecs5=genlat_MCLC(a,b,c,α) +# @show b*cos(α)/c + b^2*sin(α)^2/a^2 == 1 +# rlatvecs=get_recip_latvecs(mclc_latvecs5) +# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show γ==π/2 + +# Four triclinic cases +# kₐ>π/2,kᵦ>π/2,kᵧ>π/2 and kₐ = min(kₐ,kᵦ,kᵧ) +a=1.1 +b=1.2 +c=1.3 +α=3*π/20 +β=3.5*π/20 +γ=4*π/20 +tri_latvecs1=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs1) +((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show (α,β,γ) +# @show α>π/2 +# @show β>π/2 +# @show γ>π/2 +# @show γ==minimum([α,β,γ]) + +# kₐ>π/2,kᵦ>π/2,kᵧ>π/2 and kₐ = max(kₐ,kᵦ,kᵧ) +a=1.1 +b=1.2 +c=1.3 +α=3.5*π/20 + π/3 +β=5*π/20 + π/3 +γ=3.3*π/20 + π/3 +tri_latvecs2=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs2) +((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show (α,β,γ) +# @show α<π/2 +# @show β<π/2 +# @show γ<π/2 +# @show γ==maximum([α,β,γ]) + +# kₐ>π/2,kᵦ>π/2,kᵧ==π/2 +a=1.1 +b=1.2 +c=1.3 +α=3.5*π/20 +β=5*π/20 +γ=6*π/20-0.0188221060748125 +tri_latvecs3=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs3) +((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show γ-π/2 +# @show (α,β,γ) +# @show α>π/2 +# @show β>π/ +# @show γ≈π/2 + +# kₐ<π/2,kᵦ<π/2,kᵧ==π/2 +a=1.1 +b=1.2 +c=1.3 +α=2.5*π/20 +β=2*π/20 +γ=3.1*π/20 + 0.01079768761229742 +tri_latvecs4=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs4) +((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) +# @show γ-π/2 +# @show (α,β,γ) +# @show α>π/2 +# @show β>π/2 +# @show γ≈π/2 + +listreal_latvecs = [sqr_latvecs, rec_latvecs, reci_latvecs, hxg_latvecs, + obl_latvecs, cub_latvecs, fcc_latvecs, bcc_latvecs, tet_latvecs, + bct_latvecs1, bct_latvecs2, orc_latvecs, orcf_latvecs1, orcf_latvecs2, + orcf_latvecs3, orci_latvecs, orcc_latvecs,hex_latvecs, rhl_latvecs1, + rhl_latvecs2, mcl_latvecs, mclc_latvecs1, mclc_latvecs2, mclc_latvecs3, + mclc_latvecs4, mclc_latvecs5, tri_latvecs1, tri_latvecs2, tri_latvecs3, + tri_latvecs4] + +pgsizes = [8, 4, 4, 12, 2, 48, 48, 48, 16, 16, 16, 8, 8, 8, 8, 8, 8, 24, 12, 12, + 4, 4, 4, 4, 4, 4, 2, 2, 2, 2] + +atom_types=[0] +coords="Cartesian" +convention="ordinary" +bzformat="convex hull" +ibzformat="convex hull" + +@testset "symmetry" begin + @testset "calc_pointgroup" begin + for (i,real_latvecs)=enumerate(listreal_latvecs) + pointgroup=calc_pointgroup(real_latvecs) + @test length(pointgroup) == pgsizes[i] + end + end + + @testset "calc_bz" begin + convention="ordinary" + bzformat="convex hull" + for real_latvecs=listreal_latvecs + bz=calc_bz(real_latvecs,convention,bzformat) + pointgroup=calc_pointgroup(real_latvecs) + # Rotate BZ vertices + bzverts=reduce(hcat,[op*(bz.points[i,:]) for op=pointgroup + for i=1:size(bz.points,1)]) + # BZ vertices should map to other BZ vertices under point group + # operations + @test size(IBZ.utilities.unique(bzverts),2) == size(bz.points,1) + end + end + + @testset "calc_ibz" begin + for real_latvecs=listreal_latvecs + if size(real_latvecs) == (2,2) + atom_pos=Array([0 0]') + else + atom_pos=Array([0 0 0]') + end + pointgroup=calc_pointgroup(real_latvecs) + bz=calc_bz(real_latvecs,convention,bzformat) + ibz=calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat, + convention); + + # Unfold IBZ + unfoldpts=reduce(hcat,[op*(ibz.points[i,:]) for op=pointgroup + for i=1:size(ibz.points,1)]) + unfold_chull = chull(Array(unfoldpts')) + unfoldpts=unfold_chull.points[unfold_chull.vertices,:] + @test size(unfoldpts,1) == size(bz.points[bz.vertices,:],1) + @test all([any([isapprox(i,j) for i=1:size(unfoldpts,1)]) for + j=1:size(bz.points,1)]) + end + end +end From 68c43bcb1399d7cf50a3fa7189ab0df5a6e1b0f2 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 13:16:31 -0600 Subject: [PATCH 05/21] Added .travis.yml --- .travis.yml | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 .travis.yml diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..7b238c2 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,42 @@ +language: julia +os: + - linux + - osx + - windows + +julia: + - 1.4 + +notifications: + email: false + +jobs: + include: + - stage: "Documentation" + julia: 1.4 + os: linux + script: + - julia --project=docs/ -e 'using Pkg; + Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - julia --project=docs/ docs/make.jl + after_success: skip + +script: + - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); + Pkg.test("IBZ"; coverage=true)' + +after_success: + # Install required packages for coverage and documentation + #- julia --project -e 'import Pkg; Pkg.add("Coverage"); Pkg.add("Documenter"); + # Pkg.add("DocumenterTools");' + + # Submit test coverage report + #- julia --project -e 'using Coverage; + # Coveralls.submit(Coveralls.process_folder())' + + # Build and deploy documentation + - julia --project ./docs/make.jl + + + +coveralls: true From 62bff25d5fe218dcaf9672e2c4a7196a348c6ee9 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 14:10:24 -0600 Subject: [PATCH 06/21] Fix unit test --- .travis.yml | 4 +--- Project.toml | 3 ++- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/.travis.yml b/.travis.yml index 7b238c2..c433b89 100644 --- a/.travis.yml +++ b/.travis.yml @@ -23,7 +23,7 @@ jobs: script: - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); - Pkg.test("IBZ"; coverage=true)' + Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' after_success: # Install required packages for coverage and documentation @@ -37,6 +37,4 @@ after_success: # Build and deploy documentation - julia --project ./docs/make.jl - - coveralls: true diff --git a/Project.toml b/Project.toml index 4540a83..455829e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "IBZ" uuid = "49a35663-c880-4242-bebb-1ec8c0fa8046" -authors = ["John Christensen "] +authors = ["Jeremy Jorgensen "] version = "0.1.0" [deps] @@ -17,3 +17,4 @@ PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" QHull = "a8468747-bd6f-53ef-9e5c-744dbc5c59e7" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 9b4ee3359c1477fdad2f21e321cad60c586f8ddf Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 15:44:48 -0600 Subject: [PATCH 07/21] Remove old modules --- src/lattices.jl | 849 ------------------------------------ src/plotting.jl | 406 ----------------- src/symmetry.jl | 590 ------------------------- test/make_bz.jl | 26 -- test/make_ibz.jl | 42 -- test/make_ibz_rand.jl | 45 -- test/symmetry.jl | 245 ----------- test/testHelperFunctions.jl | 61 --- 8 files changed, 2264 deletions(-) delete mode 100644 src/lattices.jl delete mode 100644 src/plotting.jl delete mode 100644 src/symmetry.jl delete mode 100644 test/make_bz.jl delete mode 100644 test/make_ibz.jl delete mode 100644 test/make_ibz_rand.jl delete mode 100644 test/symmetry.jl delete mode 100644 test/testHelperFunctions.jl diff --git a/src/lattices.jl b/src/lattices.jl deleted file mode 100644 index a89df67..0000000 --- a/src/lattices.jl +++ /dev/null @@ -1,849 +0,0 @@ -module lattices - -using LinearAlgebra, Combinatorics - -export get_recip_latvecs, minkowski_reduce, check_reduced -# 2D Bravais lattices -export genlat_SQR, genlat_HXG, genlat_REC, genlat_RECI, genlat_OBL -# 3D Bravais lattices -export genlat_CUB, genlat_FCC, genlat_BCC, genlat_TET, genlat_BCT, genlat_ORC, - genlat_ORCF, genlat_ORCI, genlat_ORCC, genlat_HEX, genlat_RHL, genlat_MCL, - genlat_MCLC, genlat_TRI - -@doc """ - get_recip_latvecs(real_latvecs, convention) - -Calculate the reciprocal lattice vectors. - -# Arguments -- `real_latvecs::AbstractArray{<:Real,2}`: the real-space lattice vectors or - primitive translation vectors as columns of a 2x2 or 3x3 array. -- `convention::String="ordinary"`: the convention used to go between real and - reciprocal space. The two conventions are ordinary (temporal) frequency and - angular frequency. The transformation from real to reciprocal space is - unitary if the convention is ordinary. - -# Returns -- `recip_latvecs::Array{<:Real,2}` the reciprocal lattice vectors (reciprocal - primitive translation vectors) as columns of a 2x2 or 3x3 array. - -# Examples -```jldoctest -using IBZ -real_latvecs=[1 0 0; 0 1 0; 0 0 1] -convention="angular" -get_recip_latvecs(real_latvecs,convention) -# output -3×3 Array{Float64,2}: - 6.28319 0.0 0.0 - 0.0 6.28319 0.0 - 0.0 0.0 6.28319 -``` -""" -function get_recip_latvecs(real_latvecs::AbstractArray{<:Real,2}, - convention::String="ordinary")::Array{Float64,2} - if convention == "ordinary" - recip_latvecs = Array(inv(real_latvecs)') - elseif convention == "angular" - recip_latvecs = Array(2π*inv(real_latvecs)') - else - throw(ArgumentError("The allowed conventions are \"ordinary\" and - \"angular\".")) - end - recip_latvecs -end - - -""" - reduce_basis!(basis,k) - -Reduces the `k`th lattice vector. This is accomplished by locating the -lattice point closest to the projection of the `k`th lattice vector onto -the line or plane given by the other lattice vectors, subtracting the -closest lattice point from the `k`th lattice vector, and reordering the -lattice vectors by increasing Euclidean norms. - -# Arguments --`basis::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. --`k::Int`: Keeps track of which lattice vector needs to be reduced. - -# Returns --`basis::AbstractArray{<:Real,2}`: the partially reduced lattice basis as - columns of an array. --`k::Int`: The index of the lattice vector that needs to be reduced next. - -# Examples -```jldoctest -using IBZ -basis = Array([1 2 0; 0 1 0; 3 2 1]') -k=2 -k=reduce_basis!(basis,k) -basis -# output -3×3 Array{Int64,2}: - 0 1 3 - 1 2 2 - 0 0 1 -""" -function reduce_basis!(basis::AbstractArray{<:Real,2},k::Int)::Int - if k == 2 - v1,v2=[basis[:,i] for i=1:k] - i = round(Int64,dot(v1,v2)/dot(v1,v1)) - vecs = [v2-j*v1 for j=i-1:i+1] - v2 = vecs[sortperm(norm.(vecs))[1]] - if norm(v1) <= norm(v2) - basis[:,k] = v2 - k=3 - else - k=2 - basis[:,2] = v1 - basis[:,1] = v2 - end - - elseif k==3 - v1,v2,v3=[basis[:,i] for i=1:k] - - i = round(Int64,dot(v3,v1)/dot(v1,v1)) - j = round(Int64,dot(v3,v2)/dot(v2,v2)) - - vecs = [[v3-m*v2-l*v1 for m=j-1:j+1,l=i-1:i+1]...] - v3 = vecs[sortperm(norm.(vecs))[1]] - if norm(v2) <= norm(v3) - basis[:,k] = v3 - k=4 - else - if norm(v1) <= norm(v3) - k=3 - basis[:,2] = v3 - basis[:,3] = v2 - else - k=2 - basis[:,1] = v3 - basis[:,2] = v1 - basis[:,3] = v2 - end - end - else - throw(ArgumentError("Argument `k` can be 2 or 3.")) - end - k -end - -@doc """ - minkowski_reduce(basis) - -Minkowski reduce a lattice basis. Follows the logic of Fig. 4 in -\"Low-Dimensional Lattice Basis Reduction Revisited\" by Nguyen, 2009. - -# Arguments -- `basis::AbstractArray{<:Real,2}`: the lattice basis given by the columns - of a 2x2 or 3x3 array. - -# Returns -- `rbasis:: the Minkowski reduced lattice basis as columns of an array. - -# Examples -```jldoctest -using IBZ -basis = [1 2 0; 0 1 0; 0 0 1] -minkowski_reduce(basis) -# output -3×3 Array{Int64,2}: - 0 1 0 - 0 0 1 - 1 0 0 -""" -function minkowski_reduce( - basis::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} - - # Sort the lattice vectors by increasing norm. - order = sortperm([basis[:,i] for i=1:size(basis,1)]) - rbasis = basis[:,order] - - k=2 - while k <= size(rbasis,1) - k=reduce_basis!(rbasis,k) - end - rbasis -end - -@doc """ - check_reduced(basis) - -Verify a lattice basis is Minkowski reduced - -# Arguments -- `basis::AbstractArray{<:Real,2}`: the lattice basis given by the columns - of a 2x2 or 3x3 matrix. - -# Returns -- `Bool`: a boolean that indicates if the lattice basis is reduced. - -# Examples -``` jldoctest -using IBZ -basis = [1 0; 0 1] -check_reduced(basis) -# output -true -``` -""" -function check_reduced(basis::AbstractArray{<:Real,2})::Bool - if size(basis,2) == 2 - (a,b) = [basis[:,i] for i=1:2] - all([norm(a) <= norm(b), - norm(b) <= norm(b+a), - norm(b) <= norm(b-a)]) - elseif size(basis,1) == 3 - (a,b,c) = [basis[:,i] for i=1:3] - all([norm(a) <= norm(b), - norm(b) <= norm(b+a), - norm(b) <= norm(b-a), - norm(b) <= norm(c), - norm(c) <= norm(c+a), - norm(c) <= norm(c-a), - norm(c) <= norm(c+b), - norm(c) <= norm(c-b), - norm(c) <= norm(c+b-a), - norm(c) <= norm(c-b-a), - norm(c) <= norm(c+b+a), - norm(c) <= norm(c-b+a)]) - else - throw(ArgumentError("Can only verify Minkowski reduction in 2D and 3D.")) - end -end - - -@doc """ - genlat_SQR(a) - -Generate a square lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -genlat_SQR(a) -# output -2×2 Array{Int64,2}: - 1 0 - 0 1 -``` -""" -function genlat_SQR(a::Real)::AbstractArray{<:Real,2} - [a 0; 0 a] -end - -@doc """ - genlat_HXG(a) - -Generate a 2D hexagonal lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -genlat_HXG(a) -# output -2×2 Array{Float64,2}: - 1.0 -0.5 - 0.0 0.866025 -``` -""" -function genlat_HXG(a::Real)::AbstractArray{<:Real,2} - Array([a 0; -a/2 a*√3/2]') -end - -@doc """ - genlat_REC(a,b) - -Generate a rectangular lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.2 -genlat_REC(a,b) -# output -2×2 Array{Float64,2}: - 1.0 0.0 - 0.0 1.2 -``` -""" -function genlat_REC(a::Real,b::Real)::AbstractArray{<:Real,2} - Array([a 0; 0 b]') -end - -@doc """ - genlat_RECI(a,b) - -Generate a body-centered rectangular lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.2 -genlat_RECI(a,b) -# output -2×2 Array{Float64,2}: - 0.5 0.5 - -0.6 0.6 -``` -""" -function genlat_RECI(a::Real,b::Real)::AbstractArray{<:Real,2} - if a ≈ b - throw(ArgumentError("The lattice constant must be different sizes for a - rectangular lattice.")) - end - Array([a/2 -b/2; a/2 b/2]') -end - -@doc """ - genlat_OBL(a,b) - -Generate a body-centered rectangular lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.2 -θ=π/3 -test=genlat_OBL(a,b,θ) -# output -2×2 Array{Float64,2}: - 1.0 0.6 - 0.0 1.03923 -``` -""" -function genlat_OBL(a::Real,b::Real,θ::Real)::AbstractArray{<:Real,2} - if a ≈ b - throw(ArgumentError("The lattice constant must be different sizes for an - oblique lattice.")) - end - if θ ≈ π/2 - throw(ArgumentError("The lattice angle must not equal π/2 for an oblique - lattice.")) - end - Array([a 0; b*cos(θ) b*sin(θ)]') -end - -@doc """ - genlat_CUB(a) - -Generate a simple cubic lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -genlat_CUB(a) -# output -3×3 Array{Int64,2}: - 1 0 0 - 0 1 0 - 0 0 1 -``` -""" -function genlat_CUB(a::Real)::AbstractArray{<:Real,2} - [a 0 0; 0 a 0; 0 0 a] -end - -@doc """ - genlat_FCC(a) - -Generate a face-centered cubic lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -genlat_CUB(a) -# output -3×3 Array{Float64,2}: - 0.0 0.5 0.5 - 0.5 0.0 0.5 - 0.5 0.5 0.0 -""" -function genlat_FCC(a::Real)::AbstractArray{<:Real,2} - Array([0 a/2 a/2; a/2 0 a/2; a/2 a/2 0]') -end - -@doc """ - genlat_BCC(a) - -Generate a body-centered cubic lattice. - -# Arguments -- `a::Real`: the lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -genlat_BCC(a) -# output -3×3 Array{Float64,2}: - 0.0 0.5 0.5 - 0.5 0.0 0.5 - 0.5 0.5 0.0 -""" -function genlat_BCC(a::Real)::AbstractArray{<:Real,2} - a/2*Array([-1 1 1; 1 -1 1; 1 1 -1]') -end - -@doc """ - genlat_TET(a) - -Generate a simple tetragonal lattice. - -# Arguments -- `a::Real`: a lattice constant -- `c::Real`: a lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -c=1.2; -genlat_TET(a,c) -# output -3×3 Array{Float64,2}: - 1.0 0.0 0.0 - 0.0 1.0 0.0 - 0.0 0.0 1.2 -""" -function genlat_TET(a::Real,c::Real)::AbstractArray{<:Real,2} - if a ≈ c - throw(ArgumentError("The lattices constants must be different for a - tetragonal lattice.")) - end - Array([a 0 0; 0 a 0; 0 0 c]') -end - -@doc """ - genlat_BCT(a) - -Generate a body-centered tetragonal lattice. - -# Arguments -- `a::Real`: a lattice constant -- `c::Real`: a lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -c=1.2; -genlat_BCT(a,c) -# output -3×3 Array{Float64,2}: - -0.5 0.5 0.5 - 0.5 -0.5 0.5 - 0.6 0.6 -0.6 -""" -function genlat_BCT(a::Real,c::Real)::AbstractArray{<:Real,2} - if a ≈ c - throw(ArgumentError("The lattices constants must be different for a - tetragonal lattice.")) - end - Array([-a/2 a/2 c/2; a/2 -a/2 c/2; a/2 a/2 -c/2]') -end - -@doc """ - genlat_ORC(a,b,c) - -Generate an orthorhombic lattice. - -# Arguments -- `a::Real`: a lattice constant -- `b::Real`: a lattice constant -- `c::Real`: a lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.4; -c=1.2; -genlat_ORC(a,b,c) -# output -3×3 Array{Float64,2}: - 1.0 0.0 0.0 - 0.0 1.2 0.0 - 0.0 0.0 1.4 -""" -function genlat_ORC(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} - if a ≈ b || a ≈ c || b ≈ c - throw(ArgumentError("The lattices constants must all be different for an - orthorhombic lattice.")) - end - Array([a 0 0; 0 b 0; 0 0 c]') -end - -@doc """ - genlat_ORCF(a,b,c) - -Generate a face-centered orthorhombic lattice. - -# Arguments -- `a::Real`: a lattice constant -- `b::Real`: a lattice constant -- `c::Real`: a lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.4; -c=1.2; -genlat_ORCF(a,b,c) -# output -3×3 Array{Float64,2}: - 0.0 0.5 0.5 - 0.6 0.0 0.6 - 0.7 0.7 0.0 -""" -function genlat_ORCF(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} - if a ≈ b || a ≈ c || b ≈ c - throw(ArgumentError("The lattices constants must all be different for an - orthorhombic lattice.")) - end - Array([0 b/2 c/2; a/2 0 c/2; a/2 b/2 0]') -end - -@doc """ - genlat_ORCI(a,b,c) - -Generate a body-centered orthorhombic lattice. - -# Arguments -- `a::Real`: a lattice constant -- `b::Real`: a lattice constant -- `c::Real`: a lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.4; -c=1.2; -genlat_ORCI(a,b,c) -# output -3×3 Array{Float64,2}: - -0.5 0.5 0.5 - 0.6 -0.6 0.6 - 0.7 0.7 -0.7 -""" -function genlat_ORCI(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} - if a ≈ b || a ≈ c || b ≈ c - throw(ArgumentError("The lattices constants must all be different for an - orthorhombic lattice.")) - end - Array([-a/2 b/2 c/2; a/2 -b/2 c/2; a/2 b/2 -c/2]') -end - -@doc """ - genlat_ORCC(a,b,c) - -Generate a base-centered orthorhombic lattice. - -# Arguments -- `a::Real`: a lattice constant -- `b::Real`: a lattice constant -- `c::Real`: a lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.2; -c=1.4; -genlat_ORCC(a,b,c) -# output -3×3 Array{Float64,2}: - 0.5 0.5 0.0 - -0.6 0.6 0.0 - 0.0 0.0 1.4 -""" -function genlat_ORCC(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} - if a ≈ b || a ≈ c || b ≈ c - throw(ArgumentError("The lattices constants must all be different for an - orthorhombic lattice.")) - end - Array([a/2 -b/2 0; a/2 b/2 0; 0 0 c]') -end - -@doc """ - genlat_HEX(a,c) - -Generate a hexagonal lattice. - -# Arguments -- `a::Real`: a lattice constant -- `c::Real`: a lattice constant - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -c=1.2; -genlat_HEX(a,c) -# output -3×3 Array{Float64,2}: - 0.5 0.5 0.0 - -0.866025 0.866025 0.0 - 0.0 0.0 1.2 -""" -function genlat_HEX(a::Real,c::Real)::AbstractArray{<:Real,2} - if a ≈ c - throw(ArgumentError("The lattices constants must be different for a - hexagonal lattice.")) - end - Array([a/2 -a*√3/2 0; a/2 a*√3/2 0; 0 0 c]') -end - -@doc """ - genlat_RHL(a,α) - -Generate a rhombohedral lattice. - -# Arguments -- `a::Real`: a lattice constant -- `α::Real`: a lattice angle in radians - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -α=π/6; -genlat_RHL(a,α) -# output -3×3 Array{Float64,2}: - 0.965926 0.965926 0.896575 - -0.258819 0.258819 0.0 - 0.0 0.0 0.442891 -""" -function genlat_RHL(a::Real,α::Real)::AbstractArray{<:Real,2} - if α ≈ π/2 - throw(ArgumentError("The lattice angle cannot be π/2 for a rhombohedral - lattice.")) - end - Array([a*cos(α/2) -a*sin(α/2) 0; a*cos(α/2) a*sin(α/2) 0; - a*cos(α)/cos(α/2) 0 a*√(1-cos(α)^2/cos(α/2)^2)]') -end - -@doc """ - genlat_MCL(a,b,c,α) - -Generate a monoclinic lattice. - -# Arguments -- `a::Real`: a lattice constant -- `b::Real`: a lattice constant -- `c::Real`: a lattice constant -- `α::Real`: a lattice angle in radians less than π/2 - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.2 -c=1.4 -α=π/6; -genlat_MCL(a,b,c,α) -# output -3×3 Array{Float64,2}: - 1.0 0.0 0.0 - 0.0 1.2 1.21244 - 0.0 0.0 0.7 -""" -function genlat_MCL(a::Real,b::Real,c::Real,α::Real)::AbstractArray{<:Real,2} - if a ≈ b || a ≈ c || b ≈ c - throw(ArgumentError("The lattices constants must all be different for a - monoclinic lattice.")) - end - if α ≈ π/2 || α > π/2 - throw(ArgumentError("The lattice angle must be less than π/2.")) - end - Array([a 0 0; 0 b 0; 0 c*cos(α) c*sin(α)]') -end - -@doc """ - genlat_MCLC(a,b,c,α) - -Generate a base-centered monoclinic lattice. - -# Arguments -- `a::Real`: a lattice constant -- `b::Real`: a lattice constant -- `c::Real`: a lattice constant -- `α::Real`: a lattice angle in radians less than π/2 - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.2 -c=1.4 -α=π/6; -genlat_MCLC(a,b,c,α) -# output -3×3 Array{Float64,2}: - 0.5 -0.5 0.0 - 0.6 0.6 1.21244 - 0.0 0.0 0.7 -""" -function genlat_MCLC(a::Real,b::Real,c::Real,α::Real)::AbstractArray{<:Real,2} - if a ≈ b || a ≈ c || b ≈ c - throw(ArgumentError("The lattices constants must all be different for a - monoclinic lattice.")) - end - if α ≈ π/2 || α > π/2 - throw(ArgumentError("The lattice angle must be less than π/2.")) - end - Array([a/2 b/2 0; -a/2 b/2 0; 0 c*cos(α) c*sin(α)]') -end - -@doc """ - genlat_TRI(a,b,c,α,β,γ) - -Generate a triclinic lattice. - -# Arguments -- `a::Real`: a lattice constant -- `b::Real`: a lattice constant -- `c::Real`: a lattice constant -- `α::Real`: a lattice angle in radians less than π/2 - -# Returns -- `AbstractArray{<:Real,2}`: the basis of the lattice as - columns of an array. - -# Examples -```jldoctest -using IBZ -a=1 -b=1.2 -c=1.4 -α=π/6; -β=π/3; -γ=π/4; -genlat_TRI(a,b,c,α,β,γ) -# output -3×3 Array{Float64,2}: - 1.0 0.848528 0.7 - 0.0 0.848528 1.01464 - 0.0 0.0 0.663702 -""" -function genlat_TRI(a::Real,b::Real,c::Real,α::Real,β::Real, - γ::Real)::AbstractArray{<:Real,2} - if a ≈ b || a ≈ c || b ≈ c - throw(ArgumentError("The lattices constants must all be different for a - triclinic lattice.")) - end - if α ≈ β || α ≈ γ || β ≈ γ || α ≈ π/2 || β ≈ π/2 || γ ≈ π/2 - throw(ArgumentError("No lattice angles can be the same or equal to π/2 - for a triclinic lattice.")) - end - Array([a 0 0; b*cos(γ) b*sin(γ) 0; - c*cos(β) c/sin(γ)*(cos(α)-cos(β)*cos(γ)) c/sin(γ)* - √(sin(γ)^2-cos(α)^2-cos(β)^2+2*cos(α)*cos(β)*cos(γ))]') -end - -end diff --git a/src/plotting.jl b/src/plotting.jl deleted file mode 100644 index 30f5874..0000000 --- a/src/plotting.jl +++ /dev/null @@ -1,406 +0,0 @@ -module plotting - -include("symmetry.jl") -include("lattices.jl") - -using .symmetry -using .lattices - -ENV["MPLBACKEND"]="qt5agg" -using PyCall, PyPlot, QHull, LinearAlgebra - -export plot_cells - -@doc """ - function sort_points_comparison(a,b,c) - -A "less than" function for sorting Cartesian points in 2D. - -# Arguments -- `a::AbstractArray{<:Real,1}`: a point in Cartesian coordinates. -- `b::AbstractArray{<:Real,1}`: a point in Cartesian coordinates. -- `c::AbstractArray{<:Real,1}`: the position of the origin in Cartesian - coordinates. - -# Returns -- `Bool`: a boolean that indicates if `a` is less than `b`. - - -# Examples -```jldoctest -using IBZ -a=[0,0] -b=[1,0] -c=[0.5,0.5] -sort_points_comparison(a,b,c) -# output -false -``` -""" -function sort_points_comparison(a::AbstractArray{<:Real,1}, - b::AbstractArray{<:Real,1},c::AbstractArray{<:Real,1})::Bool - (ax,ay)=a - (bx,by)=b - (cx,cy)=c - - if ((ay - cy) > 0) & ((by - cy) < 0) - return true - elseif ((ay - cy) < 0) & ((by - cy) > 0) - return false - elseif ((ay - cy) >= 0) & ((by - cy) >= 0) - return bx > ax - else #((ay - cy) < 0) & ((by - cy) < 0) - return ax > bx - end -end - -@doc """ - affine_trans(pts) - -Calculate the affine transformation that maps the points to the xy-plane. - -# Arguments -- `pts::AbstractArray{<:Real,2}`: an array of Cartesian points as the columns - of a 2D array. The points must all lie on a plane in 3D. - -# Returns -- `M::AbstractArray{<:Real,2}`: the affine transformation matrix that operates - on points in homogeneous coordinates from the left. - -# Examples -```jldoctest -using IBZ -pts = [0.5 0.5 0.5; 0.5 -0.5 0.5; -0.5 0.5 0.5; -0.5 -0.5 0.5]' -affine_trans(pts) -# output -4×4 Array{Float64,2}: - 0.0 -1.0 0.0 0.5 - -1.0 0.0 0.0 0.5 - 0.0 0.0 -1.0 0.5 - 0.0 0.0 0.0 1.0 -""" -function affine_trans(pts::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} - a,b,c = [pts[:,i] for i=1:3] - - # Create a coordinate system with two vectors lying on the plane the points - # lie on. - u = b-a - v = c-a - u = u/norm(u) - v = v - (u⋅v)*u/(u⋅u) - v = v/norm(v) - w = u×v - - # Augmented matrix of affine transform - inv(vcat(hcat([u v w],a),[0 0 0 1])) -end - -@doc """ - function mapto_xyplane(pts) - -Map Cartesian points embedded in 3D to the xy-plane embedded in 2D. - -# Arguments -- `pts::AbstractArray{<:Real,2}`: Cartesian points in 3D as columns of a 2D - array. - -# Returns -- `AbstractArray{<:Real,2}`: Cartesian points in2D as columns of a 2D array. - -# Examples -```jldoctest -using IBZ -pts = [0.5 -0.5 0.5; 0.5 -0.5 -0.5; 0.5 0.5 -0.5; 0.5 0.5 0.5]' -mapto_xyplane(pts) -# output -2×4 Array{Float64,2}: - 0.0 1.0 1.0 0.0 - 0.0 0.0 1.0 1.0 -``` -""" -function mapto_xyplane(pts::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} - - M = affine_trans(pts) - reduce(hcat,[(M*[pts[:,i]..., 1])[1:2] for i=1:size(pts,2)]) -end - -@doc """ - function sortslice_perm(xypts,sxypts) - -Return the permutation vector that maps Cartesian 2D points `xypts` to `sxypts`. - -# Arguments -- `xypts::AbstractArray{<:Real,2}`: 2D Cartesian points as columns of a 2D - array. -- `sxypts::AbstractArray{<:Real,2}`: sorted 2D Cartesian points as columns of a - 2D array. - -# Returns -- `perm::AbstractArray{<:Real,1}`: a permutation vector that sorts an array of - 2D Cartesian coordinates. - -# Examples -```jldoctest -usig IBZ -xypts = [0 0; 0 1; 1 0; 1 1]' -sxypts = [0 1; 1 1; 1 0; 0 0]' -perm=sortslice_perm(xypts,sxypts) -xypts[:,perm] -# output -2×4 Array{Int64,2}: - 0 1 1 0 - 1 1 0 0 -``` -""" -function sortslice_perm(xypts::AbstractArray{<:Real,2}, - sxypts::AbstractArray{<:Real,2})::AbstractArray{<:Real,1} - perm = zeros(Int64,size(xypts,2)) - for i=1:size(xypts,2) - for j=1:size(xypts,2) - if isapprox(xypts[:,i],sxypts[:,j]) - perm[j]=i - continue - end - end - end - perm -end - -@doc """ - function sortpts_perm(pts) - -Calculate the permutation vector that sorts Cartesian points embedded in 3D that - lie on a plane (counter)clockwise with respect to the average of all points. - -# Arguments -- `pts::AbstractArray{<:Real,2}`: Cartesian points embedded in 3D that all lie - on a plane. The points are columns of a 2D array. - -# Returns -- `perm::AbstractArray{<:Real,1}`: the permutation vector that orders the points - clockwise or counterclockwise. - -# Examples -```jldoctest -using IBZ -pts = [0.5 -0.5 0.5; 0.5 -0.5 -0.5; 0.5 0.5 -0.5; 0.5 0.5 0.5]' -perm=sortpts_perm(pts) -pts[:,perm] -# output -3×4 Array{Float64,2}: - 0.5 0.5 0.5 0.5 - 0.5 0.5 -0.5 -0.5 - 0.5 -0.5 -0.5 0.5 -``` -""" -function sortpts_perm(pts::AbstractArray{<:Real,2}) - xypts=mapto_xyplane(pts) - middlept=[sum(xypts[i,:])/size(pts,2) for i=1:2] - sxypts=sortslices(xypts, lt=(x,y)->sort_points_comparison(x,y,middlept), - dims=2) - perm=sortslice_perm(xypts,sxypts) - return perm -end - -@doc """ - get_uniquefacets(ch) - -Calculate the unique facets of a convex hull object. -""" -function get_uniquefacets(ch) - facets = ch.facets - unique_facets = [] - removed=zeros(Int64,size(facets,1)) - for i=1:size(facets,1) - if removed[i] == 0 - removed[i]=1 - face=ch.simplices[i] - for j=i+1:size(facets,1) - if isapprox(facets[i,:],facets[j,:]) - removed[j]=1 - append!(face,ch.simplices[j]) - end - end - face = unique(reduce(hcat,face)[:]) - # Order the corners of the face either clockwise or - # counterclockwise. - face = face[sortpts_perm(Array(ch.points[face,:]'))] - append!(unique_facets,[face]) - end - end - unique_facets -end - -@doc """ - plot_2Dcell(convexhull, axes, color) - -Plot a 2D convex hull - -# Arguments --`convexhull::Chull{<:Real}`: a convex hull object. --`axes::PyObject`: an axes object from matplotlib. --`color::String`: the face color of the convex hull. - -# Returns --`axes::PyObject`: updated `axes` that includes a plot of the convex hull. - -# Examples -```jldoctest -using IBZ -real_latvecs = [1 0; 0 1] -convention = "ordinary" -bzformat = "convex hull" -bz = calc_bz(real_latvecs,convention,bzformat) - -fig,axes=subplots(figsize=figaspect(1)*1.5) -axes = plot_2Dcell(bz,axes,"deepskyblue"); -color="deepskyblue" -plot_2Dcell(bz,axes,color) -# output -PyObject -``` -""" -function plot_2Dcell(convexhull::Chull{<:Real}, axes::PyObject, - color::String)::PyObject - - patch=pyimport("matplotlib.patches") - collections=pyimport("matplotlib.collections") - - bzpts=Array(convexhull.points') - middlept=[sum(bzpts[i,:])/size(bzpts,2) for i=1:2] - bzpts=sortslices(bzpts, lt=(x,y)->sort_points_comparison(x,y,middlept), - dims=2) - - (x,y)=[bzpts[i,:] for i=1:2] - axes.fill(x,y, facecolor=color,edgecolor="black",linewidth=3) - axes -end - -@doc """ - plot_3Dcell(convexhull, axes, color) - -Plot a 3D convex hull - -# Arguments --`convexhull::Chull{<:Real}`: a convex hull object. --`axes::PyObject`: an axes object from matplotlib. --`color::String`: the face color of the convex hull. - -# Returns --`axes::PyObject`: updated `axes` that includes a plot of the convex hull. - -# Examples -```jldoctest -using IBZ -real_latvecs = [1 0 0; 0 1 0; 0 0 1] -convention = "ordinary" -bzformat = "convex hull" -bz = calc_bz(real_latvecs,convention,bzformat) -fig = figure(figsize=figaspect(1)*1.5) -axes = fig.add_subplot(111, projection="3d") -axes = plot_3Dcell(bz,axes,"deepskyblue") -# output -PyObject -``` -""" -function plot_3Dcell(convexhull, axes, color, plotrange=false) - - art3d=pyimport("mpl_toolkits.mplot3d.art3d") - - if plotrange ==false - ϵ=0.1*convexhull.volume - plotrange=[[minimum(convexhull.points[:,i])-ϵ, - maximum(convexhull.points[:,i])+ϵ] for i=1:3] - end - - facesᵢ=get_uniquefacets(convexhull) - edgesᵢ=deepcopy(facesᵢ) - for i=1:length(edgesᵢ) - append!(edgesᵢ[i],edgesᵢ[i][1]) - end - - faces=[convexhull.points[i,:] for i in facesᵢ] - edges=[convexhull.points[i,:] for i in edgesᵢ] - - # Reshape faces and edges - faces=[[faces[j][i,:] for i=1:size(faces[j],1)] for j=1:length(faces)] - edges=[[edges[j][i,:] for i=1:size(edges[j],1)] for j=1:length(edges)] - - p=art3d.Poly3DCollection(faces, alpha=0.2, facecolors=color) - l=art3d.Line3DCollection(edges, colors="black", linewidths=1) - - axes.add_collection3d(p) - axes.add_collection3d(l) - - #ax.view_init(-45,-45) - axes.set_xlim3d(plotrange[1]...) - axes.set_ylim3d(plotrange[2]...) - axes.set_zlim3d(plotrange[3]...) - return axes -end - -@doc """ - plot_cells(real_latvecs,atom_types,atom_pos,coords,convention,rtol,atol) - -Plot the Brillouin and Irreducible Brillouin zone in 2D or 3D. - -# Arguments --`real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as - columns of an array. --`atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. --`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal - structure as columns of an array. --`coords::String`: indicates the positions of the atoms are in \"lattice\" or - \"Cartesian\" coordinates. --`convention::String="ordinary"`: the convention used to go between real and - reciprocal space. The two conventions are ordinary (temporal) frequency and - angular frequency. The transformation from real to reciprocal space is - unitary if the convention is ordinary. --`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for - floating point comparisons. --`atol::Real=0.0`: an absolute tolerance for floating point comparisons. - -# Returns --`(fig,axes)`:the figure and axes Python objects. - -# Examples -```jldoctest -real_latvecs = [1 0; .5 1] -atom_types=[0] -coords = "Cartesian" -atom_pos = Array([0 0]') -convention = "ordinary" -(fig,axes)=plot_cells(real_latvecs,atom_types,atom_pos,coords,convention) -``` -""" -function plot_cells(real_latvecs,atom_types,atom_pos,coords,convention, - rtol::Real=sqrt(eps(float(maximum(real_latvecs)))),atol::Real=0.0) - - bzformat = "convex hull" - bz = calc_bz(real_latvecs,convention,bzformat) - ibzformat = "convex hull" - ibz = calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention, - rtol,atol) - - dim = size(real_latvecs,1) - - if dim == 2 - fig,axes=subplots(figsize=figaspect(1)*1.5) - axes = plot_2Dcell(bz,axes,"deepskyblue") - axes = plot_2Dcell(ibz,axes,"coral") - elseif dim == 3 - fig = figure(figsize=figaspect(1)*1.5) - axes = fig.add_subplot(111, projection="3d") - ϵ=0.1*bz.volume - plotrange=[[minimum(bz.points[:,i])-ϵ, - maximum(bz.points[:,i])+ϵ] for i=1:3] - axes = plot_3Dcell(bz,axes,"blue",plotrange) - axes = plot_3Dcell(ibz,axes,"red",plotrange) - else - throw(ArgumentError("The lattice vectors must be in a 2x2 or 3x3 - array.")) - end - (fig,axes) -end - -end #module diff --git a/src/symmetry.jl b/src/symmetry.jl deleted file mode 100644 index ddc39ec..0000000 --- a/src/symmetry.jl +++ /dev/null @@ -1,590 +0,0 @@ -module symmetry - -using Distances, LinearAlgebra, CDDLib, QHull, Polyhedra, Combinatorics - -include("lattices.jl") -using .lattices - -export mapto_unitcell, calc_pointgroup, calc_spacegroup, sampleCircle, - sampleSphere, calc_bz, calc_ibz - -@doc """ - shoelace(vertices) - -Calculate the area of a polygon with the shoelace algorithm. - -# Arguments -- `vertices::AbstractArray{<:Real,2}`: the xy-coordinates of the vertices - of the polygon as the columns of a 2D array. - -# Returns -- `<:Real`: the area of the polygon. -""" -function shoelace(vertices) - xs = vertices[1,:] - ys = vertices[2,:] - abs(xs[end]*ys[1] - xs[1]*ys[end] + - sum([xs[i]*ys[i+1]-xs[i+1]*ys[i] for i=1:(size(vertices,2)-1)]))/2 -end - -@doc """ - edgelengths(basis,radius) - -Calculate the edge lengths of a parallelepiped circumscribed by a sphere. - -# Arguments -- `basis::Array{<:Real,2}`: a 2x2 or 3x3 matrix whose columns give the - parallelogram or parallelepiped directions, respectively. -- `radius::Real`: the radius of the sphere. -- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for - floating point comparisons. -- `atol::Real=0.0`: an absolute tolerance for floating point - comparisons. - -# Returns -- `[la,lb,lc]::Array{Float64,1}`: a list of parallelepiped lengths. - -# Examples -```jldoctest -using IBZ -basis=Array([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]) -radius=3.0 -edgelengths(basis,radius) -# output -3-element Array{Float64,1}: - 3.0 - 3.0 - 3.0 -""" -function edgelengths(basis::Array{<:Real,2}, radius::Real, - rtol::Real=sqrt(eps(float(radius))), atol::Real=0.0)::Array{Float64,1} - - if radius < 0 | isapprox(radius, 0, rtol=rtol, atol=atol) - error("The radius has to be a positive number.") - end - - if size(basis,1) == 2 - (a,b)=[basis[:,i] for i=1:2] - ax,ay=a - bx,by=b - la=2*abs(radius*sqrt(bx^2+by^2)/(ay*bx-ax*by)) - lb=2*abs(radius*sqrt(ax^2+ay^2)/(ay*bx-ax*by)) - return [la,lb] - - elseif size(basis,1) == 3 - (a,b,c)=[basis[:,i] for i=1:3] - ax,ay,az=a - bx,by,bz=b - cx,cy,cz=c - - la=abs(radius*sqrt((by*cx-bx*cy)^2+(bz*cx-bx*cz)^2+(bz*cy-by*cz)^2)/ - (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) - lb=abs((radius*sqrt((ay*cx-ax*cy)^2+(az*cx-ax*cz)^2+(az*cy-ay*cz)^2))/ - (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) - lc=abs(radius*sqrt((ay*bx-ax*by)^2+(az*bx-ax*bz)^2+(az*by-ay*bz)^2)/ - (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) - return [la,lb,lc] - else - throw(ArgumentError("Basis has to be a 2x2 or 3x3 matrix.")) - end -end - - -@doc """ - sampleCircle(basis,radius,offset,rtol,atol) - -Sample uniformly within a circle centered about a point. - -## Arguments -- `basis::Array{<:Real,2}`: a 2x2 matrix whose columns are the grid generating - vectors. -- `radius::Real`: the radius of the circle. -- `offset::Array{<:Real,1}=[0.,0.]`: the xy-coordinates of the center of the - circle. -- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for - floating point comparisons. -- `atol::Real=0.0`: an absolute tolerance for floating point - comparisons. - -## Returns -- `pts::Array{Float64,2}` a matrix whose columns are sample points in Cartesian - coordinates. - -## Examples -```jldoctest -using IBZ -basis=Array([1. 0.; 0. 1.]') -radius=1.0 -offset=[0.,0.] -sampleCircle(basis,radius,offset) -# output -2×5 Array{Float64,2}: - 0.0 -1.0 0.0 1.0 0.0 - -1.0 0.0 0.0 0.0 1.0 -``` -""" -function sampleCircle(basis::AbstractArray{<:Real,2}, radius::Real, - offset::AbstractArray{<:Real,1}=[0.,0.], - rtol::Real=sqrt(eps(float(radius))), atol::Real=0.0)::Array{Float64,2} - - # Put the offset in lattice coordinates and round. - (o1,o2)=round.(inv(basis)*offset) - lens=edgelengths(basis,radius) - n1,n2=round.(lens) .+ 1 - - l=0; - pt=Array{Float64,1}(undef,2) - pts=Array{Float64,2}(undef,2,Int((2*n1+1)*(2*n2+1))); - distances=Array{Float64,1}(undef,size(pts,2)) - for (i,j) in Iterators.product((-n1+o1):(n1+o1),(-n2+o2):(n2+o2)) - l+=1 - pt=BLAS.gemv('N',float(basis),[i,j]) - pts[:,l]=pt - distances[l]=euclidean(pt,offset) - end - - return pts[:,distances.<=radius] - -end - -@doc """ - sampleSphere(basis,radius,offset,rtol,atol) - -Sample uniformly within a circle centered about a point. - -# Arguments -- `basis::Array{<:Real,2}`: a 3x3 matrix whose columns are the grid generating - vectors. -- `radius::Real`: the radius of the sphere. -- `offset::Array{<:Real,1}=[0.,0.]`: the xy-coordinates of the center of the - circle. -- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for - floating point comparisons. -- `atol::Real=0.0`: an absolute tolerance for floating point - comparisons. - -# Returns -- `pts::Array{Float64,2}` a matrix whose columns are sample points in Cartesian - coordinates. - -# Examples -```jldoctest -using IBZ -basis=Array([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]) -radius=1.0 -offset=[0.,0.,0.] -sampleSphere(basis,radius,offset) -# output -3×7 Array{Float64,2}: - 0.0 0.0 -1.0 0.0 1.0 0.0 0.0 - 0.0 -1.0 0.0 0.0 0.0 1.0 0.0 - -1.0 0.0 0.0 0.0 0.0 0.0 1.0 -``` -""" -function sampleSphere(basis::AbstractArray{<:Real,2}, radius::Real, - offset::Array{<:Real,1}=[0.,0.,0.], rtol::Real=sqrt(eps(float(radius))), - atol::Real=0.0)::Array{Float64,2} - - # Put the offset in lattice coordinates and round. - (o1,o2,o3)=round.(inv(basis)*offset) - lens=edgelengths(basis,radius) - n1,n2,n3=round.(lens) .+ 1 - - l=0; - pt=Array{Float64,1}(undef,3) - pts=Array{Float64,2}(undef,3,Int((2*n1+1)*(2*n2+1)*(2*n3+1))); - distances=Array{Float64,1}(undef,size(pts,2)) - for (i,j,k) in Iterators.product((-n1+o1):(n1+o1),(-n2+o2):(n2+o2), - (-n3+o3):(n3+o3)) - l+=1 - pt=BLAS.gemv('N',float(basis),[i,j,k]) - pts[:,l]=pt - distances[l]=euclidean(pt,offset) - end - - return pts[:,distances.<=radius] -end - -@doc """ - calc_pointgroup(real_latvecs,rtol,atol) - -Calculate the point group of lattice in 2D or 3D. - -# Arguments --`real_latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns - of an array. --`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))`: a relative tolerance for - floating point comparisons. It is used to compare lengths of vectors and the - volumes of primitive cells. --`atol::Real=0.0`: an absolute tolerance for floating point comparisons. It is - used to compare lengths of vectors and the volumes of primitive cells. - -# Returns --`pointgroup::Array{Array{Float64,2},1}`: the point group of the lattice. The - operators operate on points in Cartesian coordinates. - -# Examples -```jldoctest -using IBZ -basis = [1 0; 0 1] -calc_pointgroup(basis) -# output -8-element Array{Array{Float64,2},1}: - [0.0 -1.0; -1.0 0.0] - [0.0 1.0; -1.0 0.0] - [-1.0 0.0; 0.0 -1.0] - [-1.0 0.0; 0.0 1.0] - [1.0 0.0; 0.0 -1.0] - [1.0 0.0; 0.0 1.0] - [0.0 -1.0; 1.0 0.0] - [0.0 1.0; 1.0 0.0] -``` -""" -function calc_pointgroup(real_latvecs::AbstractArray{<:Real,2}, - rtol::Real=sqrt(eps(float(maximum(real_latvecs)))), - atol::Real=0.0)::Array{Array{Float64,2},1} - - dim=size(real_latvecs,1) - real_latvecs = minkowski_reduce(real_latvecs) - radius = maximum([norm(real_latvecs[:,i]) for i=1:dim]) - if dim==2 - pts=sampleCircle(real_latvecs,radius,[0,0],rtol,atol) - elseif dim==3 - pts=sampleSphere(real_latvecs,radius,[0,0,0],rtol,atol) - else - throw(ArgumentError("The lattice basis must be a 2x2 or 3x3 array.")) - end - - normsᵢ=[norm(real_latvecs[:,i]) for i=1:dim] - sizeᵢ=abs(det(real_latvecs)) - pointgroup=Array{Array{Float64,2}}([]) - inv_latvecs=inv(real_latvecs) - for perm=permutations(1:size(pts,2),dim) - real_latvecsᵣ=pts[:,perm] - normsᵣ=[norm(real_latvecsᵣ[:,i]) for i=1:dim] - sizeᵣ=abs(det(real_latvecsᵣ)) - # Point operation preserves length of lattice vectors and - # volume on primitive cell. - if (isapprox(normsᵢ,normsᵣ,rtol=rtol,atol=atol) & - isapprox(sizeᵢ,sizeᵣ,rtol=rtol,atol=atol)) - op=inv_latvecs*real_latvecsᵣ - # Point operators are orthogonal. - if isapprox(op*op',I,rtol=rtol,atol=atol) - append!(pointgroup,[op]) - end - end - end - pointgroup -end - - -@doc """ - mapto_unitcell(pt,real_latvecs,inv_latvecs,coords,rtol) - -Map a point to the first unit (primitive) cell. - -# Arguments --`pt::AbstractArray{<:Real,1}`: a point in lattice or Cartesian coordinates. --`latvecs::AbstractArray{<:Real,2}`: the basis vectors of the lattice as columns - of an array. --`inv_latvecs::AbstractArray{<:Real,2}`: the inverse of the matrix of that - contains the lattice vectors. --`coords::String`: indicates whether `pt` is in \"Cartesian\" or \"lattice\" - coordinates. --`rtol::Real=sqrt(eps(float(maximum(inv_latvecs))))`: a relative tolerance for - floating point comparisons. Finite precision errors creep up when `pt` is - transformed to lattice coordinates because the transformation requires - calculating a matrix inverse. The components of the point in lattice - coordinates are checked to ensure that values close to 1 are equal to 1. - -# Returns --`AbstractArray{<:Real,1}`: a translationally equivalent point to `pt` in the - first unit cell in the same coordinates. - -#Examples -```jldoctest -real_latvecs = [0 1 2; 0 -1 1; 1 0 0] -inv_latvecs=inv(real_latvecs) -pt=[1,2,3.2] -coords = "Cartesian" -mapto_unitcell(pt,real_latvecs,inv_latvecs,coords) -# output -3-element Array{Float64,1}: - 0.0 - 0.0 - 0.20000000000000018 -``` -""" -function mapto_unitcell(pt::AbstractArray{<:Real,1}, - latvecs::AbstractArray{<:Real,2},inv_latvecs::AbstractArray{<:Real,2}, - coords::String,rtol::Real=sqrt(eps(float(maximum(inv_latvecs)))) - )::Array{Float64,1} - if coords == "Cartesian" - Array{Float64,1}(latvecs*[isapprox(mod(c,1),1,rtol=rtol) ? 0 : mod(c,1) - for c=inv_latvecs*pt]) - elseif coords == "lattice" - Array{Float64,1}([isapprox(mod(c,1),1,rtol=rtol) ? 0 : mod(c,1) for - c=pt]) - else - throw(ArgumentError("Allowed coordinates of the point are \"Cartesian\" - and \"lattice\".")) - end -end - -@doc """ - calc_spacegroup(real_latvecs,atom_types,atom_pos,coords,rtol,atol) - -Calculate the space group of a crystal structure. - -# Arguments --`real_latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns - of an array. --`atom_types::AbstractArray{<:Int,1}`: a list of atom types as integers. --`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal - structure as columns of an array. --`coords::String`: indicates the positions of the atoms are in \"lattice\" or - \"Cartesian\" coordinates. --`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for - floating point comparisons. --`atol::Real=0.0`: an absolute tolerance for floating point comparisons. - -# Returns --`spacegroup`: the space group of the crystal structure. The first element of - `spacegroup` is a list of fractional translations, and the second element is - a list of point operators. The translations are in Cartesian coordinates, - and the operators operate on points in Cartesian coordinates. - -# Examples -```jldoctest -real_latvecs = Array([1 0; 2 1]') -atom_types = [0, 1] -atom_pos = Array([0 0; 0.5 0.5]') -coords = "Cartesian" -spacegroup=calc_spacegroup(real_latvecs,atom_types,atom_pos,coords) -# output -(Any[[0.0, 0.0], [0.0, 0.0]], Any[[1.0 0.0; 0.0 1.0], [0.0 1.0; 1.0 0.0]]) -``` -""" -function calc_spacegroup(real_latvecs::AbstractArray{<:Real,2}, - atom_types::AbstractArray{<:Int,1},atom_pos::AbstractArray{<:Real,2}, - coords::String,rtol::Real=sqrt(eps(float(maximum(real_latvecs)))), - atol::Real=0.0) - - if length(atom_types) != size(atom_pos,2) - throw(ArgumentError("The number of atom types and positions must be the - same.")) - end - - real_latvecs = minkowski_reduce(real_latvecs) - inv_latvecs = inv(real_latvecs) - pointgroup = calc_pointgroup(real_latvecs,rtol,atol) - numatoms = length(atom_types) - - # Map points to the primitive cell. - atom_pos = reduce(hcat,[mapto_unitcell(atom_pos[:,i],real_latvecs, - inv_latvecs,coords,rtol) for i=1:numatoms]) - - # Place atom positions in Cartesian coordinates. - if coords == "lattice" - atom_pos = real_latvecs*atom_pos - end - - ops_spacegroup=[] - trans_spacegroup=[] - kindᵣ=atom_types[1] - - # Atoms of the same type as the first. - same_atoms=findall(x->x==kindᵣ,atom_types) - for op in pointgroup - # Use the first atom to calculate allowed fractional translations. - posᵣ=op*atom_pos[:,1] - - for atomᵢ=same_atoms - # Calculate a candidate fractional translation. - ftrans = mapto_unitcell(atom_pos[:,atomᵢ]-posᵣ,real_latvecs, - inv_latvecs,"Cartesian") - - # Check that all atoms land on the lattice when rotated and - # translated. - mapped = false - for atomⱼ = 1:numatoms - kindⱼ = atom_types[atomⱼ] - posⱼ = op*atom_pos[:,atomⱼ] + ftrans - mapped = any([kindⱼ == atom_types[i] && - isapprox(posⱼ,atom_pos[:,i],rtol=rtol,atol=atol) ? - true : false for i=1:numatoms]) - if !mapped continue end - end - if mapped - append!(ops_spacegroup,[op]) - append!(trans_spacegroup,[ftrans]) - end - end - end - (trans_spacegroup,ops_spacegroup) -end - -@doc """ - calc_bz(real_latvecs, convention, vertsOrHrep=true, eps=1e-9) - -Calculate the Brillouin zone for the given real-spcae lattice. - -# Arguments -- `real_latvecs::AbstractArray{<:Real,2}`: the real-space lattice vectors or - primitive translation vectors as columns of a 3x3 array. -- `convention::String="ordinary"`: the convention used to go between real and - reciprocal space. The two conventions are ordinary (temporal) frequency and - angular frequency. The transformation from real to reciprocal space is - unitary if the convention is ordinary. -- `bzformat::String`: the format of the Brillouin zone. Options include - \"convex-hull\" and \"half-space\". - -# Returns -- `bz`: the vertices or half-space representation of the Brillouin zone - depending on the value of `vertsOrHrep`. - -# Examples -```jldoctest -using IBZ -real_latvecs=[1 0 0; 0 1 0; 0 0 1] -convention="ordinary" -vertsOrHrep=true -make_bz(real_latvecs,convention,vertsOrHrep) -# output -8×3 Array{Float64,2}: - 0.5 -0.5 0.5 - 0.5 -0.5 -0.5 - 0.5 0.5 -0.5 - 0.5 0.5 0.5 - -0.5 0.5 -0.5 - -0.5 0.5 0.5 - -0.5 -0.5 -0.5 - -0.5 -0.5 0.5 -``` -""" -function calc_bz(real_latvecs::AbstractArray{<:Real,2},convention::String, - bzformat::String) - - recip_latvecs = get_recip_latvecs(real_latvecs,convention) - - if size(real_latvecs) == (2,2) - latpts = reduce(hcat,[recip_latvecs*[i,j] for - (i,j)=Iterators.product(-2:2,-2:2)]) - else - latpts = reduce(hcat,[recip_latvecs*[i,j,k] for - (i,j,k)=Iterators.product(-2:2,-2:2,-2:2)]) - end - - distances = [norm(latpts[:,i]) for i=1:size(latpts,2)] .^2 ./2 - bz = HalfSpace(latpts[:,2],distances[2]) - for i=3:size(latpts,2) - bz = bz ∩ HalfSpace(latpts[:,i],distances[i]) - end - verts = reduce(hcat,points(polyhedron(bz,CDDLib.Library()))) - bzvolume = chull(Array(verts')).volume - - if !(bzvolume ≈ abs(det(recip_latvecs))) - error("The area or volume of the Brillouin zone is incorrect.") - end - - if bzformat == "half-space" - bz - elseif bzformat == "convex hull" - chull(Array(verts')) - else - throw(ArgumentError("Formats for the BZ include \"half-space\" and - \"convex hull\".")) - end -end - -@doc """ - calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention,rtol, - atol) - -Calculate the irreducible Brillouin zone of a crystal structure in 2D or 3D. - -# Arguments --`real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as - columns of an array. --`atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. --`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal - structure as columns of an array. --`coords::String`: indicates the positions of the atoms are in \"lattice\" or - \"Cartesian\" coordinates. --`ibzformat::String`: the format of the irreducible Brillouin zone. Options - include \"convex-hull\" and \"half-space\". --`convention::String="ordinary"`: the convention used to go between real and - reciprocal space. The two conventions are ordinary (temporal) frequency and - angular frequency. The transformation from real to reciprocal space is - unitary if the convention is ordinary. --`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for - floating point comparisons. --`atol::Real=0.0`: an absolute tolerance for floating point comparisons. - -# Returns --`ibz`: the irreducible Brillouin zone as a convex hull or intersection of - half-spaces. - -# Examples -```jldoctest -real_latvecs = [1 0; 0 1] -convention="ordinary" -atom_types=[0] -atom_pos = Array([0 0]') -coords = "Cartesian" -ibzformat = "convex hull" -ibz = calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention) -# output -Convex Hull of 3 points in 2 dimensions -Hull segment vertex indices: -[1, 2, 3] -Points on convex hull in original order: - -[0.0 0.0; 0.5 0.0; 0.5 0.5] -``` -""" -function calc_ibz(real_latvecs::AbstractArray{<:Real,2}, - atom_types::AbstractArray{<:Int,1},atom_pos::AbstractArray{<:Real,2}, - coords::String,ibzformat::String,convention::String="ordinary", - rtol::Real=sqrt(eps(float(maximum(real_latvecs)))), - atol::Real=0.0)::Chull{<:Real} - - pointgroup = calc_spacegroup(real_latvecs,atom_types,atom_pos,coords)[2] - copy_pg = deepcopy(pointgroup) - - bzformat = "half-space" - bz = calc_bz(real_latvecs,convention,bzformat) - bz_vertices = collect(points(polyhedron(bz,CDDLib.Library()))) - ibz = bz - while length(copy_pg) > 0 - op = pop!(copy_pg) - for v=bz_vertices - vʳ=op*v - if vʳ != v - a = vʳ-v - ibz = ibz ∩ HalfSpace(a,0) - break - end - end - end - - ibz_vertices = reduce(hcat,points(polyhedron(ibz,CDDLib.Library()))) - bz_vertices = reduce(hcat,bz_vertices) - bzvolume = chull(Array(bz_vertices')).volume - ibzvolume = chull(Array(ibz_vertices')).volume - - if !(ibzvolume ≈ bzvolume/size(pointgroup,1)) - error("The area or volume of the irreducible Brillouin zone is - incorrect.") - end - if ibzformat == "half-space" - ibz - elseif ibzformat == "convex hull" - chull(Array(ibz_vertices')) - else - throw(ArgumentError("Formats for the IBZ include \"half-space\" and - \"convex hull\".")) - end -end - -end #module diff --git a/test/make_bz.jl b/test/make_bz.jl deleted file mode 100644 index e59b27b..0000000 --- a/test/make_bz.jl +++ /dev/null @@ -1,26 +0,0 @@ -using Test,IBZ,PyCall -include("testHelperFunctions.jl") -lats = lattices -labels =lattice_labels - -@testset "make_bz_test" begin - symmetry = pyimport("phenum.symmetry") - i=1 - at = [1] - atom_pos = [0 0 0] - convention = "ordinary" - @testset "$l" for l in labels - lat = lats[i] - rlat = make_recip_latvecs(lat,convention) - bz = make_bz(lat,convention) - ch = chull(bz) - @test ch.volume ≈ det(rlat) - - spaceGroup = symmetry.get_spaceGroup(lats[i]',at,atom_pos)[1] - ops = [spaceGroup[j,:,:] for j=1:size(spaceGroup)[1]] - unfoldedpts=unique([op*bz[j,:] for op in ops, j in 1:size(bz,1)]) - @test all([any([unfoldedpts[j] ≈ bz[k,:] for j in - 1:length(unfoldedpts)]) for k in 1:size(bz,1)]) - i+=1 - end -end diff --git a/test/make_ibz.jl b/test/make_ibz.jl deleted file mode 100644 index ffbe7bf..0000000 --- a/test/make_ibz.jl +++ /dev/null @@ -1,42 +0,0 @@ -using Test,IBZ,PyCall,QHull -include("testHelperFunctions.jl") -lats = lattices -labels =lattice_labels - -@testset "bz reduce fixed" begin - i = 0 - @testset "$l" for l in labels - symmetry = pyimport("phenum.symmetry") - at = [1] - atom_pos = [0 0 0] - convention = "ordinary" - i+=1 - println(labels[i]) - lat = lats[i] - println(lat) - println("Getting BZ...") - ch = make_bz(lat,convention) - #convert to format - ch = chull(ch) - # origionalPoints = ch.points - #need to transpose for phenum lib to give the correct space group - println("Getting symmetry...") - spaceGroupTranspose = symmetry.get_spaceGroup(transpose(lat),at,atom_pos)[1] - #spaceGroup = symmetry.get_spaceGroup(lat,at,atom_pos)[1] - - println("Getting IBZ...") - rch = reduce_bz(lat,at,atom_pos,convention) - #rch_old = reduce_bz_old(lat,at,atom_pos) - - @test expectedOrder[i] ≈ size(spaceGroupTranspose,1) - - @test ch.volume/rch.volume ≈ size(spaceGroupTranspose,1) - println("Testing unflold...") - @test comparePolygon(ch, unfold(rch,spaceGroupTranspose)) - println("Testing bz mapping...") - @test BZMaps(ch,spaceGroupTranspose) - - #@test rch_old.volume/ch.volume ≈ 1/size(spaceGroupTranspose)[1] - #@test comparePolygon(ch, unfold(rch_old,spaceGroupTranspose)) - end -end diff --git a/test/make_ibz_rand.jl b/test/make_ibz_rand.jl deleted file mode 100644 index 7672462..0000000 --- a/test/make_ibz_rand.jl +++ /dev/null @@ -1,45 +0,0 @@ -using Test,IBZ,PyCall -include("testHelperFunctions.jl") -lats = lattices -labels =lattice_labels - -@testset "bz reduce rand" begin - convention = "ordinary" - lats = Array[] - test_size =20 - for i in 1:test_size push!(lats,rand_sc()) end - for i in 1:test_size push!(lats,rand_fcc()) end - for i in 1:test_size push!(lats,rand_bcc()) end - for i in 1:test_size push!(lats,rand_hex()) end - for i in 1:test_size push!(lats,rand_rhom_a()) end - for i in 1:test_size push!(lats,rand_rhom_b()) end - for i in 1:test_size push!(lats,rand_st()) end - for i in 1:test_size push!(lats,rand_bct_a()) end - for i in 1:test_size push!(lats,rand_bct_b()) end - @testset "$i" for (i,lat) in enumerate(lats) - println(lat) - at = [1] - atom_pos = [0 0 0] - symmetry = pyimport("phenum.symmetry") - - ch = make_bz(lat,convention) - #convert to format - #ch = copy(hcat(ch...)') - ch = chull(ch) - - # origionalPoints = ch.points - #need to transpose for phenum lib to give the correct space group - spaceGroupTranspose = symmetry.get_spaceGroup(transpose(lat),at,atom_pos)[1] - #println(size(spaceGroupTranspose)[1]) - #println(size(spaceGroupTranspose)[1]) - - rch = reduce_bz(lat,at,atom_pos,convention) - #rch_old = reduce_bz_old(lat,at,atom_pos) - @test rch.volume/ch.volume ≈ 1/size(spaceGroupTranspose)[1] - @test comparePolygon(ch, unfold(rch,spaceGroupTranspose)) - @test BZMaps(ch,spaceGroupTranspose) - - #@test rch_old.volume/ch.volume ≈ 1/size(spaceGroupTranspose)[1] - #@test comparePolygon(ch, unfold(rch_old,spaceGroupTranspose)) - end -end diff --git a/test/symmetry.jl b/test/symmetry.jl deleted file mode 100644 index e2c6ea4..0000000 --- a/test/symmetry.jl +++ /dev/null @@ -1,245 +0,0 @@ -using Test,IBZ - -# Lattice vectors -# 2D -a=ℯ -b=π -θ=π/3.2 -sqr_latvecs=genlat_SQR(a) -rec_latvecs=genlat_REC(a,b) -reci_latvecs=genlat_RECI(a,b) -hxg_latvecs=genlat_HXG(a) -obl_latvecs=genlat_OBL(a,b,θ) - -# 3D -a=π -cub_latvecs=genlat_CUB(a) -fcc_latvecs=genlat_FCC(a) -bcc_latvecs=genlat_BCC(a) - -c=ℯ -#1 c < a -#2 c > a -tet_latvecs=genlat_TET(a,c) -bct_latvecs1=genlat_BCT(a,c) -bct_latvecs2=genlat_BCT(c,a) - -b=5/7+(π+ℯ)/3 -orc_latvecs=genlat_ORC(a,b,c) - -# Three face-centered orthorhombic cases -# 1/a² <,=,> 1/b²+1/c² -a=π -b=5/7+(π+ℯ)/3 -c=ℯ -orcf_latvecs1=genlat_ORCF(a,b,c) -a=π/2 -b=5/7+(π+ℯ)/3 -c=ℯ -orcf_latvecs2=genlat_ORCF(a,b,c) -b=5/7 -c=ℯ -a=b*c/√(b^2 + c^2) -orcf_latvecs3=genlat_ORCF(a,b,c) - -# One body-centered orthorhombic case -a=ℯ -b=π -c=4.0+1/3 -orci_latvecs=genlat_ORCI(a,b,c) -# One base-centered orthorhombic case -orcc_latvecs=genlat_ORCC(a,b,c) -# One hexagonal case -hex_latvecs=genlat_HEX(a,c) - -# Two rhombohedral case α < π/2 and α > π/2 -α=π/7+π/2 -rhl_latvecs1=genlat_RHL(a,α) -α=π/7 -rhl_latvecs2=genlat_RHL(a,α) - -# One monoclinic case -mcl_latvecs=genlat_MCL(a,b,c,α) - -# Five base-centered monoclinic cases -# kᵧ > π/2 -a=1+1/3 -b=a+1/3 -c=b+1/3 -α=7*π/20 -mclc_latvecs1=genlat_MCLC(a,b,c,α) -# rlatvecs=get_recip_latvecs(mclc_latvecs1) -# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show γ>π/2 - -# kᵧ == π/2 -a=ℯ -b=π -c=4+1/3 -α=1.0456607472366295 -mclc_latvecs2=genlat_MCLC(a,b,c,α) -# rlatvecs=get_recip_latvecs(mclc_latvecs2) -# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show γ==π/2 - -# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 <1 -a=ℯ -b=π -c=4+1/3 -α=.3 -mclc_latvecs3=genlat_MCLC(a,b,c,α) -#@show b*cos(α)/c + b^2*sin(α)^2/a^2 <1 -# rlatvecs=get_recip_latvecs(mclc_latvecs3) -# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show γ<π/2 - -# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 -a=1.2 -b=1.3 -c=1.4 -α=0.6 -mclc_latvecs4=genlat_MCLC(a,b,c,α) -# @show b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 -# rlatvecs=get_recip_latvecs(mclc_latvecs4) -# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show γ<π/2 - -# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 -b=1.1 -c=1.2 -α=1.1 -a=(b*√(c)*sin(α))/√(c - b*cos(α)); -mclc_latvecs5=genlat_MCLC(a,b,c,α) -# @show b*cos(α)/c + b^2*sin(α)^2/a^2 == 1 -# rlatvecs=get_recip_latvecs(mclc_latvecs5) -# ((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show γ==π/2 - -# Four triclinic cases -# kₐ>π/2,kᵦ>π/2,kᵧ>π/2 and kₐ = min(kₐ,kᵦ,kᵧ) -a=1.1 -b=1.2 -c=1.3 -α=3*π/20 -β=3.5*π/20 -γ=4*π/20 -tri_latvecs1=genlat_TRI(a,b,c,α,β,γ) -rlatvecs=get_recip_latvecs(tri_latvecs1) -((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show (α,β,γ) -# @show α>π/2 -# @show β>π/2 -# @show γ>π/2 -# @show γ==minimum([α,β,γ]) - -# kₐ>π/2,kᵦ>π/2,kᵧ>π/2 and kₐ = max(kₐ,kᵦ,kᵧ) -a=1.1 -b=1.2 -c=1.3 -α=3.5*π/20 + π/3 -β=5*π/20 + π/3 -γ=3.3*π/20 + π/3 -tri_latvecs2=genlat_TRI(a,b,c,α,β,γ) -rlatvecs=get_recip_latvecs(tri_latvecs2) -((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show (α,β,γ) -# @show α<π/2 -# @show β<π/2 -# @show γ<π/2 -# @show γ==maximum([α,β,γ]) - -# kₐ>π/2,kᵦ>π/2,kᵧ==π/2 -a=1.1 -b=1.2 -c=1.3 -α=3.5*π/20 -β=5*π/20 -γ=6*π/20-0.0188221060748125 -tri_latvecs3=genlat_TRI(a,b,c,α,β,γ) -rlatvecs=get_recip_latvecs(tri_latvecs3) -((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show γ-π/2 -# @show (α,β,γ) -# @show α>π/2 -# @show β>π/ -# @show γ≈π/2 - -# kₐ<π/2,kᵦ<π/2,kᵧ==π/2 -a=1.1 -b=1.2 -c=1.3 -α=2.5*π/20 -β=2*π/20 -γ=3.1*π/20 + 0.01079768761229742 -tri_latvecs4=genlat_TRI(a,b,c,α,β,γ) -rlatvecs=get_recip_latvecs(tri_latvecs4) -((a,b,c),(α,β,γ))=IBZ.lattices.get_latparams(rlatvecs) -# @show γ-π/2 -# @show (α,β,γ) -# @show α>π/2 -# @show β>π/2 -# @show γ≈π/2 - -listreal_latvecs = [sqr_latvecs, rec_latvecs, reci_latvecs, hxg_latvecs, - obl_latvecs, cub_latvecs, fcc_latvecs, bcc_latvecs, tet_latvecs, - bct_latvecs1, bct_latvecs2, orc_latvecs, orcf_latvecs1, orcf_latvecs2, - orcf_latvecs3, orci_latvecs, orcc_latvecs,hex_latvecs, rhl_latvecs1, - rhl_latvecs2, mcl_latvecs, mclc_latvecs1, mclc_latvecs2, mclc_latvecs3, - mclc_latvecs4, mclc_latvecs5, tri_latvecs1, tri_latvecs2, tri_latvecs3, - tri_latvecs4] - -pgsizes = [8, 4, 4, 12, 2, 48, 48, 48, 16, 16, 16, 8, 8, 8, 8, 8, 8, 24, 12, 12, - 4, 4, 4, 4, 4, 4, 2, 2, 2, 2] - -atom_types=[0] -coords="Cartesian" -convention="ordinary" -bzformat="convex hull" -ibzformat="convex hull" - -@testset "symmetry" begin - @testset "calc_pointgroup" begin - for (i,real_latvecs)=enumerate(listreal_latvecs) - pointgroup=calc_pointgroup(real_latvecs) - @test length(pointgroup) == pgsizes[i] - end - end - - @testset "calc_bz" begin - convention="ordinary" - bzformat="convex hull" - for real_latvecs=listreal_latvecs - bz=calc_bz(real_latvecs,convention,bzformat) - pointgroup=calc_pointgroup(real_latvecs) - # Rotate BZ vertices - bzverts=reduce(hcat,[op*(bz.points[i,:]) for op=pointgroup - for i=1:size(bz.points,1)]) - # BZ vertices should map to other BZ vertices under point group - # operations - @test size(IBZ.utilities.unique(bzverts),2) == size(bz.points,1) - end - end - - @testset "calc_ibz" begin - for real_latvecs=listreal_latvecs - if size(real_latvecs) == (2,2) - atom_pos=Array([0 0]') - else - atom_pos=Array([0 0 0]') - end - pointgroup=calc_pointgroup(real_latvecs) - bz=calc_bz(real_latvecs,convention,bzformat) - ibz=calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat, - convention); - - # Unfold IBZ - unfoldpts=reduce(hcat,[op*(ibz.points[i,:]) for op=pointgroup - for i=1:size(ibz.points,1)]) - unfold_chull = chull(Array(unfoldpts')) - unfoldpts=unfold_chull.points[unfold_chull.vertices,:] - @test size(unfoldpts,1) == size(bz.points[bz.vertices,:],1) - @test all([any([isapprox(i,j) for i=1:size(unfoldpts,1)]) for - j=1:size(bz.points,1)]) - end - end -end diff --git a/test/testHelperFunctions.jl b/test/testHelperFunctions.jl deleted file mode 100644 index b35bf45..0000000 --- a/test/testHelperFunctions.jl +++ /dev/null @@ -1,61 +0,0 @@ -export unfold, BZMaps,comparePolygon - -function unfold(IBZ,spaceGroup) - ops = [spaceGroup[i,:,:] for i=1:size(spaceGroup)[1]] - points = IBZ.points - for m in ops - #println(m) - newPoints = Array{Float64,2}(undef,0,3) - - for i in 1:size(points)[1] - #println(points[i,:]) - newVert = transpose(points[i,:]) * m - #println(newVert) - newPoints = vcat(newPoints,newVert) - - - end - points = vcat(points, newPoints) - points = reducePoints(points) - - end - return chull(points) -end -function BZMaps(BZ,spaceGroup) - ops = [spaceGroup[i,:,:] for i=1:size(spaceGroup)[1]] - for m in ops - for i in 1:size(BZ.points)[1] - #newVert = transpose(transpose(BZ.points[i,:]) * m) - newVert = m*BZ.points[i,:] - #check if the newVert maps to another point - point_maps = false - for j in 1:size(BZ.points)[1] - if newVert ≈ BZ.points[j,:] - point_maps = true - end - end - if !point_maps - return false - end - end - end - return true -end -function comparePolygon(bz1,bz2) - num_points_in_2 = 0 - for i in bz1.vertices - point = bz1.points[i,:] - point_in_2 = false - for j in bz2.vertices - compare_point = bz2.points[j,:] - if compare_point ≈ point - point_in_2 = true - num_points_in_2 +=1 - end - end - if !point_in_2 - return false - end - end - return num_points_in_2 == size(bz2.vertices)[1] -end From 40ada1f078e82601f9382ded6c94186a4826fa16 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 15:49:42 -0600 Subject: [PATCH 08/21] Add unit tests --- src/Lattices.jl | 903 +++++++++++++++++++++++++++++++++++++++++++++++ src/Plotting.jl | 402 +++++++++++++++++++++ src/Symmetry.jl | 597 +++++++++++++++++++++++++++++++ test/symmetry.jl | 245 +++++++++++++ 4 files changed, 2147 insertions(+) create mode 100644 src/Lattices.jl create mode 100644 src/Plotting.jl create mode 100644 src/Symmetry.jl create mode 100644 test/symmetry.jl diff --git a/src/Lattices.jl b/src/Lattices.jl new file mode 100644 index 0000000..6c600b0 --- /dev/null +++ b/src/Lattices.jl @@ -0,0 +1,903 @@ +module Lattices + +using LinearAlgebra, Combinatorics + +export get_recip_latvecs, minkowski_reduce, check_reduced, get_latparams +# 2D Bravais lattices +export genlat_SQR, genlat_HXG, genlat_REC, genlat_RECI, genlat_OBL +# 3D Bravais lattices +export genlat_CUB, genlat_FCC, genlat_BCC, genlat_TET, genlat_BCT, genlat_ORC, + genlat_ORCF, genlat_ORCI, genlat_ORCC, genlat_HEX, genlat_RHL, genlat_MCL, + genlat_MCLC, genlat_TRI + +@doc """ + get_recip_latvecs(real_latvecs, convention) + +Calculate the reciprocal lattice vectors. + +# Arguments +- `real_latvecs::AbstractArray{<:Real,2}`: the real-space lattice vectors or + primitive translation vectors as columns of a 2x2 or 3x3 array. +- `convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. + +# Returns +- `recip_latvecs::Array{<:Real,2}` the reciprocal lattice vectors (reciprocal + primitive translation vectors) as columns of a 2x2 or 3x3 array. + +# Examples +```jldoctest +using IBZ +real_latvecs=[1 0 0; 0 1 0; 0 0 1] +convention="angular" +IBZ.Lattices.get_recip_latvecs(real_latvecs,convention) +# output +3×3 Array{Float64,2}: + 6.28319 0.0 0.0 + 0.0 6.28319 0.0 + 0.0 0.0 6.28319 +``` +""" +function get_recip_latvecs(real_latvecs::AbstractArray{<:Real,2}, + convention::String="ordinary")::Array{Float64,2} + if convention == "ordinary" + recip_latvecs = Array(inv(real_latvecs)') + elseif convention == "angular" + recip_latvecs = Array(2π*inv(real_latvecs)') + else + throw(ArgumentError("The allowed conventions are \"ordinary\" and + \"angular\".")) + end + recip_latvecs +end + +@doc """ + get_latparams(latvecs) + +Calculate the lattice constants and angles of a lattice basis. + +# Arguments +-`latvecs::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. + +# Returns +- A list where the first element is a list lattice constants `(a,b,c)` and second + lattice angles in radians `(α,β,γ)`. + +# Examples +```jldoctest +using IBZ +latvecs = [1 0; 0 1] +IBZ.Lattices.get_latparams(latvecs) +# output +2-element Array{Array{Float64,1},1}: + [1.0, 1.0] + [1.5707963267948966, 1.5707963267948966] +``` +""" +function get_latparams(latvecs) + + if size(latvecs) == (2,2) + (a,b)=[latvecs[:,i] for i=1:2] + θ=acos(dot(a,b)/(norm(a)*norm(b))) + [[norm(a),norm(b)],[θ,θ]] + elseif size(latvecs) == (3,3) + (a,b,c)=[latvecs[:,i] for i=1:3] + α=acos(dot(b,c)/(norm(b)*norm(c))) + β=acos(dot(a,c)/(norm(a)*norm(c))) + γ=acos(dot(a,b)/(norm(a)*norm(b))) + [[norm(a),norm(b),norm(c)],[α,β,γ]] + else + throw(ArgumentError("The lattice vectors must be a 2x2 or 3x3 array.")) + end +end + +""" + reduce_basis!(basis,k) + +Reduces the `k`th lattice vector. This is accomplished by locating the +lattice point closest to the projection of the `k`th lattice vector onto +the line or plane given by the other lattice vector(s), subtracting the +closest lattice point from the `k`th lattice vector, and reordering the +lattice vectors by increasing Euclidean norms. + +# Arguments +-`basis::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. +-`k::Int`: Keeps track of which lattice vector needs to be reduced. + +# Returns +-`basis::AbstractArray{<:Real,2}`: the partially reduced lattice basis as + columns of an array. +-`k::Int`: The index of the lattice vector that needs to be reduced next. + +# Examples +```jldoctest +using IBZ +basis = Array([1 2 0; 0 1 0; 3 2 1]') +k=2 +IBZ.Lattices.reduce_basis!(basis,k) +basis +# output +3×3 Array{Int64,2}: + 0 1 3 + 1 2 2 + 0 0 1 +``` +""" +function reduce_basis!(basis::AbstractArray{<:Real,2},k::Int)::Int + if k == 2 + v1,v2=[basis[:,i] for i=1:k] + i = round(Int64,dot(v1,v2)/dot(v1,v1)) + vecs = [v2-j*v1 for j=i-1:i+1] + v2 = vecs[sortperm(norm.(vecs))[1]] + if norm(v1) <= norm(v2) + basis[:,k] = v2 + k=3 + else + k=2 + basis[:,2] = v1 + basis[:,1] = v2 + end + + elseif k==3 + v1,v2,v3=[basis[:,i] for i=1:k] + + i = round(Int64,dot(v3,v1)/dot(v1,v1)) + j = round(Int64,dot(v3,v2)/dot(v2,v2)) + + vecs = [[v3-m*v2-l*v1 for m=j-1:j+1,l=i-1:i+1]...] + v3 = vecs[sortperm(norm.(vecs))[1]] + if norm(v2) <= norm(v3) + basis[:,k] = v3 + k=4 + else + if norm(v1) <= norm(v3) + k=3 + basis[:,2] = v3 + basis[:,3] = v2 + else + k=2 + basis[:,1] = v3 + basis[:,2] = v1 + basis[:,3] = v2 + end + end + else + throw(ArgumentError("Argument `k` can be 2 or 3.")) + end + k +end + +@doc """ + minkowski_reduce(basis) + +Minkowski reduce a lattice basis. Follows the logic of Fig. 4 in +\"Low-Dimensional Lattice Basis Reduction Revisited\" by Nguyen, 2009. + +# Arguments +- `basis::AbstractArray{<:Real,2}`: the lattice basis given by the columns + of a 2x2 or 3x3 array. + +# Returns +- `rbasis:: the Minkowski reduced lattice basis as columns of an array. + +# Examples +```jldoctest +using IBZ +basis = [1 2 0; 0 1 0; 0 0 1] +IBZ.Lattices.minkowski_reduce(basis) +# output +3×3 Array{Int64,2}: + 0 1 0 + 0 0 1 + 1 0 0 +``` +""" +function minkowski_reduce( + basis::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} + + # Sort the lattice vectors by increasing norm. + order = sortperm([basis[:,i] for i=1:size(basis,1)]) + rbasis = basis[:,order] + + k=2 + while k <= size(rbasis,1) + k=reduce_basis!(rbasis,k) + end + rbasis +end + +@doc """ + check_reduced(basis) + +Verify a lattice basis is Minkowski reduced + +# Arguments +- `basis::AbstractArray{<:Real,2}`: the lattice basis given by the columns + of a 2x2 or 3x3 matrix. + +# Returns +- `Bool`: a boolean that indicates if the lattice basis is reduced. + +# Examples +``` jldoctest +using IBZ +basis = [1 0; 0 1] +IBZ.Lattices.check_reduced(basis) +# output +true +``` +""" +function check_reduced(basis::AbstractArray{<:Real,2})::Bool + if size(basis,2) == 2 + (a,b) = [basis[:,i] for i=1:2] + all([norm(a) <= norm(b), + norm(b) <= norm(b+a), + norm(b) <= norm(b-a)]) + elseif size(basis,1) == 3 + (a,b,c) = [basis[:,i] for i=1:3] + all([norm(a) <= norm(b), + norm(b) <= norm(b+a), + norm(b) <= norm(b-a), + norm(b) <= norm(c), + norm(c) <= norm(c+a), + norm(c) <= norm(c-a), + norm(c) <= norm(c+b), + norm(c) <= norm(c-b), + norm(c) <= norm(c+b-a), + norm(c) <= norm(c-b-a), + norm(c) <= norm(c+b+a), + norm(c) <= norm(c-b+a)]) + else + throw(ArgumentError("Can only verify Minkowski reduction in 2D and 3D.")) + end +end + + +@doc """ + genlat_SQR(a) + +Generate a square lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +IBZ.Lattices.genlat_SQR(a) +# output +2×2 Array{Int64,2}: + 1 0 + 0 1 +``` +""" +function genlat_SQR(a::Real)::AbstractArray{<:Real,2} + [a 0; 0 a] +end + +@doc """ + genlat_HXG(a) + +Generate a 2D hexagonal lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +IBZ.Lattices.genlat_HXG(a) +# output +2×2 Array{Float64,2}: + 1.0 -0.5 + 0.0 0.866025 +``` +""" +function genlat_HXG(a::Real)::AbstractArray{<:Real,2} + Array([a 0; -a/2 a*√3/2]') +end + +@doc """ + genlat_REC(a,b) + +Generate a rectangular lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +IBZ.Lattices.genlat_REC(a,b) +# output +2×2 Array{Float64,2}: + 1.0 0.0 + 0.0 1.2 +``` +""" +function genlat_REC(a::Real,b::Real)::AbstractArray{<:Real,2} + Array([a 0; 0 b]') +end + +@doc """ + genlat_RECI(a,b) + +Generate a body-centered rectangular lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +IBZ.Lattices.genlat_RECI(a,b) +# output +2×2 Array{Float64,2}: + 0.5 0.5 + -0.6 0.6 +``` +""" +function genlat_RECI(a::Real,b::Real)::AbstractArray{<:Real,2} + if a ≈ b + throw(ArgumentError("The lattice constant must be different sizes for a + rectangular lattice.")) + end + Array([a/2 -b/2; a/2 b/2]') +end + +@doc """ + genlat_OBL(a,b) + +Generate a body-centered rectangular lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +θ=π/3 +IBZ.Lattices.genlat_OBL(a,b,θ) +# output +2×2 Array{Float64,2}: + 1.0 0.6 + 0.0 1.03923 +``` +""" +function genlat_OBL(a::Real,b::Real,θ::Real)::AbstractArray{<:Real,2} + if a ≈ b + throw(ArgumentError("The lattice constant must be different sizes for an + oblique lattice.")) + end + if θ ≈ π/2 + throw(ArgumentError("The lattice angle must not equal π/2 for an oblique + lattice.")) + end + Array([a 0; b*cos(θ) b*sin(θ)]') +end + +@doc """ + genlat_CUB(a) + +Generate a simple cubic lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +IBZ.Lattices.genlat_CUB(a) +# output +3×3 Array{Int64,2}: + 1 0 0 + 0 1 0 + 0 0 1 +``` +""" +function genlat_CUB(a::Real)::AbstractArray{<:Real,2} + [a 0 0; 0 a 0; 0 0 a] +end + +@doc """ + genlat_FCC(a) + +Generate a face-centered cubic lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +IBZ.Lattices.genlat_FCC(a) +# output +3×3 Array{Float64,2}: + 0.0 0.5 0.5 + 0.5 0.0 0.5 + 0.5 0.5 0.0 +``` +""" +function genlat_FCC(a::Real)::AbstractArray{<:Real,2} + Array([0 a/2 a/2; a/2 0 a/2; a/2 a/2 0]') +end + +@doc """ + genlat_BCC(a) + +Generate a body-centered cubic lattice. + +# Arguments +- `a::Real`: the lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +IBZ.Lattices.genlat_BCC(a) +# output +3×3 Array{Float64,2}: + -0.5 0.5 0.5 + 0.5 -0.5 0.5 + 0.5 0.5 -0.5 +``` +""" +function genlat_BCC(a::Real)::AbstractArray{<:Real,2} + a/2*Array([-1 1 1; 1 -1 1; 1 1 -1]') +end + +@doc """ + genlat_TET(a) + +Generate a simple tetragonal lattice. + +# Arguments +- `a::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +c=1.2; +IBZ.Lattices.genlat_TET(a,c) +# output +3×3 Array{Float64,2}: + 1.0 0.0 0.0 + 0.0 1.0 0.0 + 0.0 0.0 1.2 +``` +""" +function genlat_TET(a::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ c + throw(ArgumentError("The lattices constants must be different for a + tetragonal lattice.")) + end + Array([a 0 0; 0 a 0; 0 0 c]') +end + +@doc """ + genlat_BCT(a) + +Generate a body-centered tetragonal lattice. + +# Arguments +- `a::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +c=1.2; +IBZ.Lattices.genlat_BCT(a,c) +# output +3×3 Array{Float64,2}: + -0.5 0.5 0.5 + 0.5 -0.5 0.5 + 0.6 0.6 -0.6 +``` +""" +function genlat_BCT(a::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ c + throw(ArgumentError("The lattices constants must be different for a + tetragonal lattice.")) + end + Array([-a/2 a/2 c/2; a/2 -a/2 c/2; a/2 a/2 -c/2]') +end + +@doc """ + genlat_ORC(a,b,c) + +Generate an orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.4; +c=1.2; +IBZ.Lattices.genlat_ORC(a,b,c) +# output +3×3 Array{Float64,2}: + 1.0 0.0 0.0 + 0.0 1.4 0.0 + 0.0 0.0 1.2 +``` +""" +function genlat_ORC(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([a 0 0; 0 b 0; 0 0 c]') +end + +@doc """ + genlat_ORCF(a,b,c) + +Generate a face-centered orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.4; +c=1.2; +IBZ.Lattices.genlat_ORCF(a,b,c) +# output +3×3 Array{Float64,2}: + 0.0 0.5 0.5 + 0.7 0.0 0.7 + 0.6 0.6 0.0 +``` +""" +function genlat_ORCF(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([0 b/2 c/2; a/2 0 c/2; a/2 b/2 0]') +end + +@doc """ + genlat_ORCI(a,b,c) + +Generate a body-centered orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.4; +c=1.2; +IBZ.Lattices.genlat_ORCI(a,b,c) +# output +3×3 Array{Float64,2}: + -0.5 0.5 0.5 + 0.7 -0.7 0.7 + 0.6 0.6 -0.6 +``` +""" +function genlat_ORCI(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([-a/2 b/2 c/2; a/2 -b/2 c/2; a/2 b/2 -c/2]') +end + +@doc """ + genlat_ORCC(a,b,c) + +Generate a base-centered orthorhombic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2; +c=1.4; +IBZ.Lattices.genlat_ORCC(a,b,c) +# output +3×3 Array{Float64,2}: + 0.5 0.5 0.0 + -0.6 0.6 0.0 + 0.0 0.0 1.4 +``` +""" +function genlat_ORCC(a::Real,b::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for an + orthorhombic lattice.")) + end + Array([a/2 -b/2 0; a/2 b/2 0; 0 0 c]') +end + +@doc """ + genlat_HEX(a,c) + +Generate a hexagonal lattice. + +# Arguments +- `a::Real`: a lattice constant +- `c::Real`: a lattice constant + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +c=1.2; +IBZ.Lattices.genlat_HEX(a,c) +# output +3×3 Array{Float64,2}: + 0.5 0.5 0.0 + -0.866025 0.866025 0.0 + 0.0 0.0 1.2 +``` +""" +function genlat_HEX(a::Real,c::Real)::AbstractArray{<:Real,2} + if a ≈ c + throw(ArgumentError("The lattices constants must be different for a + hexagonal lattice.")) + end + Array([a/2 -a*√3/2 0; a/2 a*√3/2 0; 0 0 c]') +end + +@doc """ + genlat_RHL(a,α) + +Generate a rhombohedral lattice. + +# Arguments +- `a::Real`: a lattice constant +- `α::Real`: a lattice angle in radians + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +α=π/6; +IBZ.Lattices.genlat_RHL(a,α) +# output +3×3 Array{Float64,2}: + 0.965926 0.965926 0.896575 + -0.258819 0.258819 0.0 + 0.0 0.0 0.442891 +``` +""" +function genlat_RHL(a::Real,α::Real)::AbstractArray{<:Real,2} + if α ≈ π/2 + throw(ArgumentError("The lattice angle cannot be π/2 for a rhombohedral + lattice.")) + end + Array([a*cos(α/2) -a*sin(α/2) 0; a*cos(α/2) a*sin(α/2) 0; + a*cos(α)/cos(α/2) 0 a*√(1-(cos(α)/cos(α/2))^2)]') +end + +@doc """ + genlat_MCL(a,b,c,α) + +Generate a monoclinic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant +- `α::Real`: a lattice angle in radians less than π/2 + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +c=1.4 +α=π/6; +IBZ.Lattices.genlat_MCL(a,b,c,α) +# output +3×3 Array{Float64,2}: + 1.0 0.0 0.0 + 0.0 1.2 1.21244 + 0.0 0.0 0.7 +``` +""" +function genlat_MCL(a::Real,b::Real,c::Real,α::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for a + monoclinic lattice.")) + end + if α ≈ π/2 || α > π/2 + throw(ArgumentError("The lattice angle must be less than π/2.")) + end + Array([a 0 0; 0 b 0; 0 c*cos(α) c*sin(α)]') +end + +@doc """ + genlat_MCLC(a,b,c,α) + +Generate a base-centered monoclinic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant +- `α::Real`: a lattice angle in radians less than π/2 + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +c=1.4 +α=π/6; +IBZ.Lattices.genlat_MCLC(a,b,c,α) +# output +3×3 Array{Float64,2}: + 0.5 -0.5 0.0 + 0.6 0.6 1.21244 + 0.0 0.0 0.7 +``` +""" +function genlat_MCLC(a::Real,b::Real,c::Real,α::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for a + monoclinic lattice.")) + end + if α ≈ π/2 || α > π/2 + throw(ArgumentError("The lattice angle must be less than π/2.")) + end + Array([a/2 b/2 0; -a/2 b/2 0; 0 c*cos(α) c*sin(α)]') +end + +@doc """ + genlat_TRI(a,b,c,α,β,γ) + +Generate a triclinic lattice. + +# Arguments +- `a::Real`: a lattice constant +- `b::Real`: a lattice constant +- `c::Real`: a lattice constant +- `α::Real`: a lattice angle in radians less than π/2 + +# Returns +- `AbstractArray{<:Real,2}`: the basis of the lattice as + columns of an array. + +# Examples +```jldoctest +using IBZ +a=1 +b=1.2 +c=1.4 +α=π/6; +β=π/3; +γ=π/4; +IBZ.Lattices.genlat_TRI(a,b,c,α,β,γ) +# output +3×3 Array{Float64,2}: + 1.0 0.848528 0.7 + 0.0 0.848528 1.01464 + 0.0 0.0 0.663702 +``` +""" +function genlat_TRI(a::Real,b::Real,c::Real,α::Real,β::Real, + γ::Real)::AbstractArray{<:Real,2} + if a ≈ b || a ≈ c || b ≈ c + throw(ArgumentError("The lattices constants must all be different for a + triclinic lattice.")) + end + if α ≈ β || α ≈ γ || β ≈ γ + throw(ArgumentError("No lattice angles can be the same for a triclinic + lattice.")) + end + Array([a 0 0; b*cos(γ) b*sin(γ) 0; + c*cos(β) c/sin(γ)*(cos(α)-cos(β)*cos(γ)) c/sin(γ)* + √(sin(γ)^2-cos(α)^2-cos(β)^2+2*cos(α)*cos(β)*cos(γ))]') +end + +end diff --git a/src/Plotting.jl b/src/Plotting.jl new file mode 100644 index 0000000..5027f9d --- /dev/null +++ b/src/Plotting.jl @@ -0,0 +1,402 @@ +module Plotting + +include("Symmetry.jl") +include("Lattices.jl") + +using .Symmetry +using .Lattices + +ENV["MPLBACKEND"]="qt5agg" +using PyCall, PyPlot, QHull, LinearAlgebra + +export plot_cells + +@doc """ + function sort_points_comparison(a,b,c) + +A "less than" function for sorting Cartesian points in 2D. + +# Arguments +- `a::AbstractArray{<:Real,1}`: a point in Cartesian coordinates. +- `b::AbstractArray{<:Real,1}`: a point in Cartesian coordinates. +- `c::AbstractArray{<:Real,1}`: the position of the origin in Cartesian + coordinates. + +# Returns +- `Bool`: a boolean that indicates if `a` is less than `b`. + + +# Examples +```jldoctest +using IBZ +a=[0,0] +b=[1,0] +c=[0.5,0.5] +IBZ.Plotting.sort_points_comparison(a,b,c) +# output +false +``` +""" +function sort_points_comparison(a::AbstractArray{<:Real,1}, + b::AbstractArray{<:Real,1},c::AbstractArray{<:Real,1})::Bool + (ax,ay)=a + (bx,by)=b + (cx,cy)=c + + if ((ay - cy) > 0) & ((by - cy) < 0) + return true + elseif ((ay - cy) < 0) & ((by - cy) > 0) + return false + elseif ((ay - cy) >= 0) & ((by - cy) >= 0) + return bx > ax + else #((ay - cy) < 0) & ((by - cy) < 0) + return ax > bx + end +end + +@doc """ + affine_trans(pts) + +Calculate the affine transformation that maps the points to the xy-plane. + +# Arguments +- `pts::AbstractArray{<:Real,2}`: an array of Cartesian points as the columns + of a 2D array. The points must all lie on a plane in 3D. + +# Returns +- `M::AbstractArray{<:Real,2}`: the affine transformation matrix that operates + on points in homogeneous coordinates from the left. + +# Examples +```jldoctest +using IBZ +pts = [0.5 0.5 0.5; 0.5 -0.5 0.5; -0.5 0.5 0.5; -0.5 -0.5 0.5]' +IBZ.Plotting.affine_trans(pts) +# output +4×4 Array{Float64,2}: + 0.0 -1.0 0.0 0.5 + -1.0 0.0 0.0 0.5 + 0.0 0.0 -1.0 0.5 + 0.0 0.0 0.0 1.0 +``` +""" +function affine_trans(pts::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} + a,b,c = [pts[:,i] for i=1:3] + + # Create a coordinate system with two vectors lying on the plane the points + # lie on. + u = b-a + v = c-a + u = u/norm(u) + v = v - (u⋅v)*u/(u⋅u) + v = v/norm(v) + w = u×v + + # Augmented matrix of affine transform + inv(vcat(hcat([u v w],a),[0 0 0 1])) +end + +@doc """ + function mapto_xyplane(pts) + +Map Cartesian points embedded in 3D to the xy-plane embedded in 2D. + +# Arguments +- `pts::AbstractArray{<:Real,2}`: Cartesian points in 3D as columns of a 2D + array. + +# Returns +- `AbstractArray{<:Real,2}`: Cartesian points in2D as columns of a 2D array. + +# Examples +```jldoctest +using IBZ +pts = [0.5 -0.5 0.5; 0.5 -0.5 -0.5; 0.5 0.5 -0.5; 0.5 0.5 0.5]' +IBZ.Plotting.mapto_xyplane(pts) +# output +2×4 Array{Float64,2}: + 0.0 1.0 1.0 0.0 + 0.0 0.0 1.0 1.0 +``` +""" +function mapto_xyplane(pts::AbstractArray{<:Real,2})::AbstractArray{<:Real,2} + + M = affine_trans(pts) + reduce(hcat,[(M*[pts[:,i]..., 1])[1:2] for i=1:size(pts,2)]) +end + +@doc """ + function sortslice_perm(xypts,sxypts) + +Return the permutation vector that maps Cartesian 2D points `xypts` to `sxypts`. + +# Arguments +- `xypts::AbstractArray{<:Real,2}`: 2D Cartesian points as columns of a 2D + array. +- `sxypts::AbstractArray{<:Real,2}`: sorted 2D Cartesian points as columns of a + 2D array. + +# Returns +- `perm::AbstractArray{<:Real,1}`: a permutation vector that sorts an array of + 2D Cartesian coordinates. + +# Examples +```jldoctest +using IBZ +xypts = [0 0; 0 1; 1 0; 1 1]' +sxypts = [0 1; 1 1; 1 0; 0 0]' +perm=IBZ.Plotting.sortslice_perm(xypts,sxypts) +xypts[:,perm] +# output +2×4 Array{Int64,2}: + 0 1 1 0 + 1 1 0 0 +``` +""" +function sortslice_perm(xypts::AbstractArray{<:Real,2}, + sxypts::AbstractArray{<:Real,2})::AbstractArray{<:Real,1} + perm = zeros(Int64,size(xypts,2)) + for i=1:size(xypts,2) + for j=1:size(xypts,2) + if isapprox(xypts[:,i],sxypts[:,j]) + perm[j]=i + continue + end + end + end + perm +end + +@doc """ + function sortpts_perm(pts) + +Calculate the permutation vector that sorts Cartesian points embedded in 3D that + lie on a plane (counter)clockwise with respect to the average of all points. + +# Arguments +- `pts::AbstractArray{<:Real,2}`: Cartesian points embedded in 3D that all lie + on a plane. The points are columns of a 2D array. + +# Returns +- `perm::AbstractArray{<:Real,1}`: the permutation vector that orders the points + clockwise or counterclockwise. + +# Examples +```jldoctest +using IBZ +pts = [0.5 -0.5 0.5; 0.5 -0.5 -0.5; 0.5 0.5 -0.5; 0.5 0.5 0.5]' +perm=IBZ.Plotting.sortpts_perm(pts) +pts[:,perm] +# output +3×4 Array{Float64,2}: + 0.5 0.5 0.5 0.5 + 0.5 0.5 -0.5 -0.5 + 0.5 -0.5 -0.5 0.5 +``` +""" +function sortpts_perm(pts::AbstractArray{<:Real,2}) + xypts=mapto_xyplane(pts) + middlept=[sum(xypts[i,:])/size(pts,2) for i=1:2] + sxypts=sortslices(xypts, lt=(x,y)->sort_points_comparison(x,y,middlept), + dims=2) + perm=sortslice_perm(xypts,sxypts) + return perm +end + +@doc """ + get_uniquefacets(ch) + +Calculate the unique facets of a convex hull object. +""" +function get_uniquefacets(ch) + facets = ch.facets + unique_facets = [] + removed=zeros(Int64,size(facets,1)) + for i=1:size(facets,1) + if removed[i] == 0 + removed[i]=1 + face=ch.simplices[i] + for j=i+1:size(facets,1) + if isapprox(facets[i,:],facets[j,:],rtol=1e-6) + removed[j]=1 + append!(face,ch.simplices[j]) + end + end + face = unique(reduce(hcat,face)[:]) + # Order the corners of the face either clockwise or + # counterclockwise. + face = face[sortpts_perm(Array(ch.points[face,:]'))] + append!(unique_facets,[face]) + end + end + unique_facets +end + +@doc """ + plot_2Dcell(convexhull, ax, color) + +Plot a 2D convex hull + +# Arguments +- `convexhull::Chull{<:Real}`: a convex hull object. +- `ax::PyObject`: an axes object from matplotlib. +- `color::String`: the face color of the convex hull. + +# Returns +- `ax::PyObject`: updated `ax` that includes a plot of the convex hull. + +# Examples +``` +using IBZ +real_latvecs = [1 0; 0 1] +convention = "ordinary" +bzformat = "convex hull" +bz = calc_bz(real_latvecs,convention,bzformat) + +fig,ax=subplots(figsize=figaspect(1)*1.5) +ax = IBZ.Plotting.plot_2Dcell(bz,ax,"deepskyblue"); +color="deepskyblue" +IBZ.Plotting.plot_2Dcell(bz,ax,color) +``` +""" +function plot_2Dcell(convexhull::Chull{<:Real}, ax::PyObject, + color::String)::PyObject + + patch=pyimport("matplotlib.patches") + collections=pyimport("matplotlib.collections") + + bzpts=Array(convexhull.points') + middlept=[sum(bzpts[i,:])/size(bzpts,2) for i=1:2] + bzpts=sortslices(bzpts, lt=(x,y)->sort_points_comparison(x,y,middlept), + dims=2) + + (x,y)=[bzpts[i,:] for i=1:2] + ax.fill(x,y, facecolor=color,edgecolor="black",linewidth=3) + ax +end + +@doc """ + plot_3Dcell(convexhull, ax, color) + +Plot a 3D convex hull + +# Arguments +- `convexhull::Chull{<:Real}`: a convex hull object. +- `ax::PyObject`: an axes object from matplotlib. +- `color::String`: the face color of the convex hull. + +# Returns +- `ax::PyObject`: updated `ax` that includes a plot of the convex hull. + +# Examples +``` +using IBZ +real_latvecs = [1 0 0; 0 1 0; 0 0 1] +convention = "ordinary" +bzformat = "convex hull" +bz = calc_bz(real_latvecs,convention,bzformat) +fig = figure(figsize=figaspect(1)*1.5) +ax = fig.add_subplot(111, projection="3d") +ax = IBZ.Plotting.plot_3Dcell(bz,ax,"deepskyblue") +``` +""" +function plot_3Dcell(convexhull, ax, color, plotrange=false) + + art3d=pyimport("mpl_toolkits.mplot3d.art3d") + + if plotrange ==false + ϵ=0.1*convexhull.volume + plotrange=[[minimum(convexhull.points[:,i])-ϵ, + maximum(convexhull.points[:,i])+ϵ] for i=1:3] + end + + facesᵢ=get_uniquefacets(convexhull) + edgesᵢ=deepcopy(facesᵢ) + for i=1:length(edgesᵢ) + append!(edgesᵢ[i],edgesᵢ[i][1]) + end + + faces=[convexhull.points[i,:] for i in facesᵢ] + edges=[convexhull.points[i,:] for i in edgesᵢ] + + # Reshape faces and edges + faces=[[faces[j][i,:] for i=1:size(faces[j],1)] for j=1:length(faces)] + edges=[[edges[j][i,:] for i=1:size(edges[j],1)] for j=1:length(edges)] + + p=art3d.Poly3DCollection(faces, alpha=0.2, facecolors=color) + l=art3d.Line3DCollection(edges, colors="black", linewidths=1) + + ax.add_collection3d(p) + ax.add_collection3d(l) + + #ax.view_init(-45,-45) + ax.set_xlim3d(plotrange[1]...) + ax.set_ylim3d(plotrange[2]...) + ax.set_zlim3d(plotrange[3]...) + return ax +end + +@doc """ + plot_cells(real_latvecs,atom_types,atom_pos,coords,convention,rtol,atol) + +Plot the Brillouin and Irreducible Brillouin zone in 2D or 3D. + +# Arguments +- `real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as + columns of an array. +- `atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. +- `atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal + structure as columns of an array. +- `coords::String`: indicates the positions of the atoms are in \"lattice\" or + \"Cartesian\" coordinates. +- `convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. +- `rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for + floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point comparisons. + +# Returns +- `(fig,ax)`: the figure and axes Python objects. + +# Examples +``` +using IBZ +real_latvecs = [1 0; .5 1] +atom_types=[0] +coords = "Cartesian" +atom_pos = Array([0 0]') +convention = "ordinary" +(fig,ax)=plot_cells(real_latvecs,atom_types,atom_pos,coords,convention) +``` +""" +function plot_cells(real_latvecs,atom_types,atom_pos,coords,convention, + rtol::Real=sqrt(eps(float(maximum(real_latvecs)))),atol::Real=0.0) + bzformat = "convex hull" + bz = calc_bz(real_latvecs,convention,bzformat) + ibzformat = "convex hull" + ibz = calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention, + rtol,atol) + + dim = size(real_latvecs,1) + if dim == 2 + fig,ax=subplots(figsize=figaspect(1)*1.5) + ax = plot_2Dcell(bz,ax,"deepskyblue") + ax = plot_2Dcell(ibz,ax,"coral") + elseif dim == 3 + fig = figure(figsize=figaspect(1)*1.5) + ax = fig.add_subplot(111, projection="3d") + ϵ=0.1*bz.volume + plotrange=[[minimum(bz.points[:,i])-ϵ, + maximum(bz.points[:,i])+ϵ] for i=1:3] + ax = plot_3Dcell(bz,ax,"blue",plotrange) + ax = plot_3Dcell(ibz,ax,"red",plotrange) + else + throw(ArgumentError("The lattice vectors must be in a 2x2 or 3x3 + array.")) + end + (fig,ax) +end + +end #module diff --git a/src/Symmetry.jl b/src/Symmetry.jl new file mode 100644 index 0000000..6d3958b --- /dev/null +++ b/src/Symmetry.jl @@ -0,0 +1,597 @@ +module Symmetry + +using Distances, LinearAlgebra, CDDLib, QHull, Polyhedra, Combinatorics + +include("Lattices.jl") +using .Lattices + +export mapto_unitcell, calc_pointgroup, calc_spacegroup, sampleCircle, + sampleSphere, calc_bz, calc_ibz + +@doc """ + shoelace(vertices) + +Calculate the area of a polygon with the shoelace algorithm. + +# Arguments +- `vertices::AbstractArray{<:Real,2}`: the xy-coordinates of the vertices + of the polygon as the columns of a 2D array. + +# Returns +- `<:Real`: the area of the polygon. +""" +function shoelace(vertices) + xs = vertices[1,:] + ys = vertices[2,:] + abs(xs[end]*ys[1] - xs[1]*ys[end] + + sum([xs[i]*ys[i+1]-xs[i+1]*ys[i] for i=1:(size(vertices,2)-1)]))/2 +end + +@doc """ + edgelengths(basis,radius) + +Calculate the edge lengths of a parallelepiped circumscribed by a sphere. + +# Arguments +- `basis::Array{<:Real,2}`: a 2x2 or 3x3 matrix whose columns give the + parallelogram or parallelepiped directions, respectively. +- `radius::Real`: the radius of the sphere. +- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for + floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point + comparisons. + +# Returns +- `[la,lb,lc]::Array{Float64,1}`: a list of parallelepiped lengths. + +# Examples +```jldoctest +using IBZ +basis=Array([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]) +radius=3.0 +IBZ.Symmetry.edgelengths(basis,radius) +# output +3-element Array{Float64,1}: + 3.0 + 3.0 + 3.0 +``` +""" +function edgelengths(basis::Array{<:Real,2}, radius::Real, + rtol::Real=sqrt(eps(float(radius))), atol::Real=0.0)::Array{Float64,1} + + if radius < 0 | isapprox(radius, 0, rtol=rtol, atol=atol) + error("The radius has to be a positive number.") + end + + if size(basis,1) == 2 + (a,b)=[basis[:,i] for i=1:2] + ax,ay=a + bx,by=b + la=2*abs(radius*sqrt(bx^2+by^2)/(ay*bx-ax*by)) + lb=2*abs(radius*sqrt(ax^2+ay^2)/(ay*bx-ax*by)) + return [la,lb] + + elseif size(basis,1) == 3 + (a,b,c)=[basis[:,i] for i=1:3] + ax,ay,az=a + bx,by,bz=b + cx,cy,cz=c + + la=abs(radius*sqrt((by*cx-bx*cy)^2+(bz*cx-bx*cz)^2+(bz*cy-by*cz)^2)/ + (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) + lb=abs((radius*sqrt((ay*cx-ax*cy)^2+(az*cx-ax*cz)^2+(az*cy-ay*cz)^2))/ + (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) + lc=abs(radius*sqrt((ay*bx-ax*by)^2+(az*bx-ax*bz)^2+(az*by-ay*bz)^2)/ + (az*by*cx-ay*bz*cx-az*bx*cy+ax*bz*cy+ay*bx*cz-ax*by*cz)) + return [la,lb,lc] + else + throw(ArgumentError("Basis has to be a 2x2 or 3x3 matrix.")) + end +end + + +@doc """ + sampleCircle(basis,radius,offset,rtol,atol) + +Sample uniformly within a circle centered about a point. + +## Arguments +- `basis::Array{<:Real,2}`: a 2x2 matrix whose columns are the grid generating + vectors. +- `radius::Real`: the radius of the circle. +- `offset::Array{<:Real,1}=[0.,0.]`: the xy-coordinates of the center of the + circle. +- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for + floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point + comparisons. + +## Returns +- `pts::Array{Float64,2}` a matrix whose columns are sample points in Cartesian + coordinates. + +## Examples +```jldoctest +using IBZ +basis=Array([1. 0.; 0. 1.]') +radius=1.0 +offset=[0.,0.] +IBZ.Symmetry.sampleCircle(basis,radius,offset) +# output +2×5 Array{Float64,2}: + 0.0 -1.0 0.0 1.0 0.0 + -1.0 0.0 0.0 0.0 1.0 +``` +""" +function sampleCircle(basis::AbstractArray{<:Real,2}, radius::Real, + offset::AbstractArray{<:Real,1}=[0.,0.], + rtol::Real=sqrt(eps(float(radius))), atol::Real=0.0)::Array{Float64,2} + + # Put the offset in lattice coordinates and round. + (o1,o2)=round.(inv(basis)*offset) + lens=edgelengths(basis,radius) + n1,n2=round.(lens) .+ 1 + + l=0; + pt=Array{Float64,1}(undef,2) + pts=Array{Float64,2}(undef,2,Int((2*n1+1)*(2*n2+1))); + distances=Array{Float64,1}(undef,size(pts,2)) + for (i,j) in Iterators.product((-n1+o1):(n1+o1),(-n2+o2):(n2+o2)) + l+=1 + pt=BLAS.gemv('N',float(basis),[i,j]) + pts[:,l]=pt + distances[l]=euclidean(pt,offset) + end + + return pts[:,distances.<=radius] + +end + +@doc """ + sampleSphere(basis,radius,offset,rtol,atol) + +Sample uniformly within a circle centered about a point. + +# Arguments +- `basis::Array{<:Real,2}`: a 3x3 matrix whose columns are the grid generating + vectors. +- `radius::Real`: the radius of the sphere. +- `offset::Array{<:Real,1}=[0.,0.]`: the xy-coordinates of the center of the + circle. +- `rtol::Real=sqrt(eps(float(radius)))`: a relative tolerace for + floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point + comparisons. + +# Returns +- `pts::Array{Float64,2}` a matrix whose columns are sample points in Cartesian + coordinates. + +# Examples +```jldoctest +using IBZ +basis=Array([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]) +radius=1.0 +offset=[0.,0.,0.] +IBZ.Symmetry.sampleSphere(basis,radius,offset) +# output +3×7 Array{Float64,2}: + 0.0 0.0 -1.0 0.0 1.0 0.0 0.0 + 0.0 -1.0 0.0 0.0 0.0 1.0 0.0 + -1.0 0.0 0.0 0.0 0.0 0.0 1.0 +``` +""" +function sampleSphere(basis::AbstractArray{<:Real,2}, radius::Real, + offset::Array{<:Real,1}=[0.,0.,0.], rtol::Real=sqrt(eps(float(radius))), + atol::Real=0.0)::Array{Float64,2} + + # Put the offset in lattice coordinates and round. + (o1,o2,o3)=round.(inv(basis)*offset) + lens=edgelengths(basis,radius) + n1,n2,n3=round.(lens) .+ 1 + + l=0; + pt=Array{Float64,1}(undef,3) + pts=Array{Float64,2}(undef,3,Int((2*n1+1)*(2*n2+1)*(2*n3+1))); + distances=Array{Float64,1}(undef,size(pts,2)) + for (i,j,k) in Iterators.product((-n1+o1):(n1+o1),(-n2+o2):(n2+o2), + (-n3+o3):(n3+o3)) + l+=1 + pt=BLAS.gemv('N',float(basis),[i,j,k]) + pts[:,l]=pt + distances[l]=euclidean(pt,offset) + end + + pts[:,findall(x->(xx==kindᵣ,atom_types) + for op in pointgroup + # Use the first atom to calculate allowed fractional translations. + posᵣ=op*atom_pos[:,1] + + for atomᵢ=same_atoms + # Calculate a candidate fractional translation. + ftrans = mapto_unitcell(atom_pos[:,atomᵢ]-posᵣ,real_latvecs, + inv_latvecs,"Cartesian") + + # Check that all atoms land on the lattice when rotated and + # translated. + mapped = false + for atomⱼ = 1:numatoms + kindⱼ = atom_types[atomⱼ] + posⱼ = op*atom_pos[:,atomⱼ] + ftrans + mapped = any([kindⱼ == atom_types[i] && + isapprox(posⱼ,atom_pos[:,i],rtol=rtol,atol=atol) ? + true : false for i=1:numatoms]) + if !mapped continue end + end + if mapped + append!(ops_spacegroup,[op]) + append!(trans_spacegroup,[ftrans]) + end + end + end + (trans_spacegroup,ops_spacegroup) +end + +@doc """ + calc_bz(real_latvecs, convention, vertsOrHrep=true, eps=1e-9) + +Calculate the Brillouin zone for the given real-spcae lattice. + +# Arguments +- `real_latvecs::AbstractArray{<:Real,2}`: the real-space lattice vectors or + primitive translation vectors as columns of a 3x3 array. +- `convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. +- `bzformat::String`: the format of the Brillouin zone. Options include + \"convex-hull\" and \"half-space\". + +# Returns +- `bz`: the vertices or half-space representation of the Brillouin zone + depending on the value of `vertsOrHrep`. + +# Examples +```jldoctest +using IBZ +real_latvecs=[1 0; 0 1] +convention="ordinary" +bzformat = "convex hull" +IBZ.Symmetry.calc_bz(real_latvecs,convention,bzformat) +# output +Convex Hull of 4 points in 2 dimensions +Hull segment vertex indices: +[3, 2, 1, 4] +Points on convex hull in original order: + +[0.5 0.5; 0.5 -0.5; -0.5 -0.5; -0.5 0.5] +``` +""" +function calc_bz(real_latvecs::AbstractArray{<:Real,2},convention::String, + bzformat::String) + + recip_latvecs = minkowski_reduce(get_recip_latvecs(real_latvecs,convention)) + + if size(real_latvecs) == (2,2) + latpts = reduce(hcat,[recip_latvecs*[i,j] for + (i,j)=Iterators.product(-2:2,-2:2)]) + else + latpts = reduce(hcat,[recip_latvecs*[i,j,k] for + (i,j,k)=Iterators.product(-2:2,-2:2,-2:2)]) + end + + distances = [norm(latpts[:,i]) for i=1:size(latpts,2)] .^2 ./2 + bz = HalfSpace(latpts[:,2],distances[2]) + for i=3:size(latpts,2) + bz = bz ∩ HalfSpace(latpts[:,i],distances[i]) + end + verts = reduce(hcat,points(polyhedron(bz,CDDLib.Library()))) + bzvolume = chull(Array(verts')).volume + + if !(bzvolume ≈ abs(det(recip_latvecs))) + error("The area or volume of the Brillouin zone is incorrect.") + end + + if bzformat == "half-space" + bz + elseif bzformat == "convex hull" + chull(Array(verts')) + else + throw(ArgumentError("Formats for the BZ include \"half-space\" and + \"convex hull\".")) + end +end + +@doc """ + calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat,convention,rtol, + atol) + +Calculate the irreducible Brillouin zone of a crystal structure in 2D or 3D. + +# Arguments +-`real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as + columns of an array. +-`atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. +-`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal + structure as columns of an array. +-`coords::String`: indicates the positions of the atoms are in \"lattice\" or + \"Cartesian\" coordinates. +-`ibzformat::String`: the format of the irreducible Brillouin zone. Options + include \"convex-hull\" and \"half-space\". +-`convention::String="ordinary"`: the convention used to go between real and + reciprocal space. The two conventions are ordinary (temporal) frequency and + angular frequency. The transformation from real to reciprocal space is + unitary if the convention is ordinary. +-`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for + floating point comparisons. +-`atol::Real=0.0`: an absolute tolerance for floating point comparisons. + +# Returns +-`ibz`: the irreducible Brillouin zone as a convex hull or intersection of + half-spaces. + +# Examples +```jldoctest +using IBZ +real_latvecs = [1 0; 0 1] +convention="ordinary" +atom_types=[0] +atom_pos = Array([0 0]') +coords = "Cartesian" +ibzformat = "convex hull" +IBZ.Symmetry.calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat, + convention) +# output +Convex Hull of 3 points in 2 dimensions +Hull segment vertex indices: +[1, 2, 3] +Points on convex hull in original order: + +[0.0 0.0; 0.5 0.0; 0.5 0.5] +``` +""" +function calc_ibz(real_latvecs::AbstractArray{<:Real,2}, + atom_types::AbstractArray{<:Int,1},atom_pos::AbstractArray{<:Real,2}, + coords::String,ibzformat::String,convention::String="ordinary", + rtol::Real=sqrt(eps(float(maximum(real_latvecs)))), + atol::Real=0.0)::Chull{<:Real} + pointgroup = calc_spacegroup(real_latvecs,atom_types,atom_pos,coords)[2] + copy_pg = deepcopy(pointgroup) + bzformat = "half-space" + bz = calc_bz(real_latvecs,convention,bzformat) + bz_vertices = collect(points(polyhedron(bz,CDDLib.Library()))) + ibz = bz + while length(copy_pg) > 0 + op = pop!(copy_pg) + for v=bz_vertices + vʳ=op*v + if !isapprox(vʳ,v,rtol=rtol,atol=atol) + a = vʳ-v + ibz = ibz ∩ HalfSpace(a,0) + break + end + end + end + ibz_vertices = reduce(hcat,points(polyhedron(ibz,CDDLib.Library()))) + bz_vertices = reduce(hcat,bz_vertices) + bzvolume = chull(Array(bz_vertices')).volume + ibzvolume = chull(Array(ibz_vertices')).volume + reduction = bzvolume/ibzvolume + if !(ibzvolume ≈ bzvolume/size(pointgroup,1)) + error("The area or volume of the irreducible Brillouin zone is + incorrect.") + end + if ibzformat == "half-space" + ibz + elseif ibzformat == "convex hull" + chull(Array(ibz_vertices')) + else + throw(ArgumentError("Formats for the IBZ include \"half-space\" and + \"convex hull\".")) + end +end + +num_ops = [8, 4, 12, 2, 2, 4, 8, 12, 16, 24, 48] +lat_types = ["square", "rectangular", "triangular", "oblique", "triclinic", + "monoclinic", "orthorhombic", "rhombohedral", "tetragonal", "hexagonal", + "cubic"] + +"""Give the size of the point group of a Bravais lattice.""" +pointgroup_size = Dict{String,Integer}(i=>j for (i,j)=zip(lat_types, num_ops)) + +end #module diff --git a/test/symmetry.jl b/test/symmetry.jl new file mode 100644 index 0000000..a0f4f96 --- /dev/null +++ b/test/symmetry.jl @@ -0,0 +1,245 @@ +using Test,IBZ + +# Lattice vectors +# 2D +a=ℯ +b=π +θ=π/3.2 +sqr_latvecs=genlat_SQR(a) +rec_latvecs=genlat_REC(a,b) +reci_latvecs=genlat_RECI(a,b) +hxg_latvecs=genlat_HXG(a) +obl_latvecs=genlat_OBL(a,b,θ) + +# 3D +a=π +cub_latvecs=genlat_CUB(a) +fcc_latvecs=genlat_FCC(a) +bcc_latvecs=genlat_BCC(a) + +c=ℯ +#1 c < a +#2 c > a +tet_latvecs=genlat_TET(a,c) +bct_latvecs1=genlat_BCT(a,c) +bct_latvecs2=genlat_BCT(c,a) + +b=5/7+(π+ℯ)/3 +orc_latvecs=genlat_ORC(a,b,c) + +# Three face-centered orthorhombic cases +# 1/a² <,=,> 1/b²+1/c² +a=π +b=5/7+(π+ℯ)/3 +c=ℯ +orcf_latvecs1=genlat_ORCF(a,b,c) +a=π/2 +b=5/7+(π+ℯ)/3 +c=ℯ +orcf_latvecs2=genlat_ORCF(a,b,c) +b=5/7 +c=ℯ +a=b*c/√(b^2 + c^2) +orcf_latvecs3=genlat_ORCF(a,b,c) + +# One body-centered orthorhombic case +a=ℯ +b=π +c=4.0+1/3 +orci_latvecs=genlat_ORCI(a,b,c) +# One base-centered orthorhombic case +orcc_latvecs=genlat_ORCC(a,b,c) +# One hexagonal case +hex_latvecs=genlat_HEX(a,c) + +# Two rhombohedral case α < π/2 and α > π/2 +α=π/7+π/2 +rhl_latvecs1=genlat_RHL(a,α) +α=π/7 +rhl_latvecs2=genlat_RHL(a,α) + +# One monoclinic case +mcl_latvecs=genlat_MCL(a,b,c,α) + +# Five base-centered monoclinic cases +# kᵧ > π/2 +a=1+1/3 +b=a+1/3 +c=b+1/3 +α=7*π/20 +mclc_latvecs1=genlat_MCLC(a,b,c,α) +# rlatvecs=get_recip_latvecs(mclc_latvecs1) +# ((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show γ>π/2 + +# kᵧ == π/2 +a=ℯ +b=π +c=4+1/3 +α=1.0456607472366295 +mclc_latvecs2=genlat_MCLC(a,b,c,α) +# rlatvecs=get_recip_latvecs(mclc_latvecs2) +# ((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show γ==π/2 + +# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 <1 +a=ℯ +b=π +c=4+1/3 +α=.3 +mclc_latvecs3=genlat_MCLC(a,b,c,α) +#@show b*cos(α)/c + b^2*sin(α)^2/a^2 <1 +# rlatvecs=get_recip_latvecs(mclc_latvecs3) +# ((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show γ<π/2 + +# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 +a=1.2 +b=1.3 +c=1.4 +α=0.6 +mclc_latvecs4=genlat_MCLC(a,b,c,α) +# @show b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 +# rlatvecs=get_recip_latvecs(mclc_latvecs4) +# ((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show γ<π/2 + +# kᵧ < π/2, b*cos(α)/c + b^2*sin(α)^2/a^2 > 1 +b=1.1 +c=1.2 +α=1.1 +a=(b*√(c)*sin(α))/√(c - b*cos(α)); +mclc_latvecs5=genlat_MCLC(a,b,c,α) +# @show b*cos(α)/c + b^2*sin(α)^2/a^2 == 1 +# rlatvecs=get_recip_latvecs(mclc_latvecs5) +# ((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show γ==π/2 + +# Four triclinic cases +# kₐ>π/2,kᵦ>π/2,kᵧ>π/2 and kₐ = min(kₐ,kᵦ,kᵧ) +a=1.1 +b=1.2 +c=1.3 +α=3*π/20 +β=3.5*π/20 +γ=4*π/20 +tri_latvecs1=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs1) +((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show (α,β,γ) +# @show α>π/2 +# @show β>π/2 +# @show γ>π/2 +# @show γ==minimum([α,β,γ]) + +# kₐ>π/2,kᵦ>π/2,kᵧ>π/2 and kₐ = max(kₐ,kᵦ,kᵧ) +a=1.1 +b=1.2 +c=1.3 +α=3.5*π/20 + π/3 +β=5*π/20 + π/3 +γ=3.3*π/20 + π/3 +tri_latvecs2=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs2) +((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show (α,β,γ) +# @show α<π/2 +# @show β<π/2 +# @show γ<π/2 +# @show γ==maximum([α,β,γ]) + +# kₐ>π/2,kᵦ>π/2,kᵧ==π/2 +a=1.1 +b=1.2 +c=1.3 +α=3.5*π/20 +β=5*π/20 +γ=6*π/20-0.0188221060748125 +tri_latvecs3=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs3) +((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show γ-π/2 +# @show (α,β,γ) +# @show α>π/2 +# @show β>π/ +# @show γ≈π/2 + +# kₐ<π/2,kᵦ<π/2,kᵧ==π/2 +a=1.1 +b=1.2 +c=1.3 +α=2.5*π/20 +β=2*π/20 +γ=3.1*π/20 + 0.01079768761229742 +tri_latvecs4=genlat_TRI(a,b,c,α,β,γ) +rlatvecs=get_recip_latvecs(tri_latvecs4) +((a,b,c),(α,β,γ))=IBZ.Lattices.get_latparams(rlatvecs) +# @show γ-π/2 +# @show (α,β,γ) +# @show α>π/2 +# @show β>π/2 +# @show γ≈π/2 + +listreal_latvecs = [sqr_latvecs, rec_latvecs, reci_latvecs, hxg_latvecs, + obl_latvecs, cub_latvecs, fcc_latvecs, bcc_latvecs, tet_latvecs, + bct_latvecs1, bct_latvecs2, orc_latvecs, orcf_latvecs1, orcf_latvecs2, + orcf_latvecs3, orci_latvecs, orcc_latvecs,hex_latvecs, rhl_latvecs1, + rhl_latvecs2, mcl_latvecs, mclc_latvecs1, mclc_latvecs2, mclc_latvecs3, + mclc_latvecs4, mclc_latvecs5, tri_latvecs1, tri_latvecs2, tri_latvecs3, + tri_latvecs4] + +pgsizes = [8, 4, 4, 12, 2, 48, 48, 48, 16, 16, 16, 8, 8, 8, 8, 8, 8, 24, 12, 12, + 4, 4, 4, 4, 4, 4, 2, 2, 2, 2] + +atom_types=[0] +coords="Cartesian" +convention="ordinary" +bzformat="convex hull" +ibzformat="convex hull" + +@testset "symmetry" begin + @testset "calc_pointgroup" begin + for (i,real_latvecs)=enumerate(listreal_latvecs) + pointgroup=calc_pointgroup(real_latvecs) + @test length(pointgroup) == pgsizes[i] + end + end + + @testset "calc_bz" begin + convention="ordinary" + bzformat="convex hull" + for real_latvecs=listreal_latvecs + bz=calc_bz(real_latvecs,convention,bzformat) + pointgroup=calc_pointgroup(real_latvecs) + # Rotate BZ vertices + bzverts=reduce(hcat,[op*(bz.points[i,:]) for op=pointgroup + for i=1:size(bz.points,1)]) + # BZ vertices should map to other BZ vertices under point group + # operations + @test size(IBZ.Utilities.unique(bzverts),2) == size(bz.points,1) + end + end + + @testset "calc_ibz" begin + for real_latvecs=listreal_latvecs + if size(real_latvecs) == (2,2) + atom_pos=Array([0 0]') + else + atom_pos=Array([0 0 0]') + end + pointgroup=calc_pointgroup(real_latvecs) + bz=calc_bz(real_latvecs,convention,bzformat) + ibz=calc_ibz(real_latvecs,atom_types,atom_pos,coords,ibzformat, + convention); + + # Unfold IBZ + unfoldpts=reduce(hcat,[op*(ibz.points[i,:]) for op=pointgroup + for i=1:size(ibz.points,1)]) + unfold_chull = chull(Array(unfoldpts')) + unfoldpts=unfold_chull.points[unfold_chull.vertices,:] + @test size(unfoldpts,1) == size(bz.points[bz.vertices,:],1) + @test all([any([isapprox(unfoldpts[i,:],bz.points[j,:]) + for i=1:size(unfoldpts,1)]) for j=1:size(bz.points,1)]) + end + end +end From 034f5af464c94916a4fb804d5d118b4e3d826ae1 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 16:16:12 -0600 Subject: [PATCH 09/21] update .travis.yml --- .travis.yml | 49 +++++++++++++++++++++++-------------------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/.travis.yml b/.travis.yml index c433b89..f96bc52 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,14 +1,4 @@ language: julia -os: - - linux - - osx - - windows - -julia: - - 1.4 - -notifications: - email: false jobs: include: @@ -16,25 +6,32 @@ jobs: julia: 1.4 os: linux script: + # Install required packages for coverage and documentation + - julia --project -e 'import Pkg; Pkg.add("Coverage"); + Pkg.add("Documenter"); Pkg.add("DocumenterTools");' - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - julia --project=docs/ docs/make.jl - after_success: skip - -script: - - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); - Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' -after_success: - # Install required packages for coverage and documentation - #- julia --project -e 'import Pkg; Pkg.add("Coverage"); Pkg.add("Documenter"); - # Pkg.add("DocumenterTools");' + # Submit test coverage report + - julia --project -e 'using Coverage; + Coveralls.submit(Coveralls.process_folder())' - # Submit test coverage report - #- julia --project -e 'using Coverage; - # Coveralls.submit(Coveralls.process_folder())' + after_success: + - julia --project ./docs/make.jl + + coveralls: true - # Build and deploy documentation - - julia --project ./docs/make.jl - -coveralls: true + - stage: "Tests" + install: + - pip3 install scipy + julia: 1.4 + os: + - linux + - osx + - windows + notifications: + email: false + script: + - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); + Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' From 2e734b6858734e3db1ab996bca4ab17675c58ec8 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 19:01:08 -0600 Subject: [PATCH 10/21] Add documentation --- .travis.yml | 32 +++++++++++------------ docs/make.jl | 12 +++++++++ docs/src/Documentation.md | 53 +++++++++++++++++++++++++++++++++++++++ docs/src/IBZ.md | 0 docs/src/Usage.md | 38 ++++++++++++++++++++++++++++ docs/src/index.md | 9 +++++++ src/Lattices.jl | 10 ++++---- src/Symmetry.jl | 52 +++++++++++++++++++------------------- 8 files changed, 159 insertions(+), 47 deletions(-) create mode 100644 docs/make.jl create mode 100644 docs/src/Documentation.md create mode 100644 docs/src/IBZ.md create mode 100644 docs/src/Usage.md create mode 100644 docs/src/index.md diff --git a/.travis.yml b/.travis.yml index f96bc52..590e119 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,10 +1,24 @@ language: julia +jobs: + include: + - stage: "Tests" + install: + - pip3 install scipy + julia: 1.4 + os: + - osx + notifications: + email: false + script: + - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); + Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' + jobs: include: - stage: "Documentation" julia: 1.4 - os: linux + os: osx script: # Install required packages for coverage and documentation - julia --project -e 'import Pkg; Pkg.add("Coverage"); @@ -19,19 +33,5 @@ jobs: after_success: - julia --project ./docs/make.jl - - coveralls: true - - stage: "Tests" - install: - - pip3 install scipy - julia: 1.4 - os: - - linux - - osx - - windows - notifications: - email: false - script: - - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); - Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' + coveralls: true diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..cd8b780 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,12 @@ +push!(LOAD_PATH,"../src/") +using Documenter, IBZ + +makedocs(sitename="IBZ Documentation", + format = Documenter.HTML(prettyurls = false), + modules = [IBZ], + authors = "Jeremy Jorgensen", + pages=["index.md", "Documentation.md", "Usage.md"]) + +# deploydocs( +# repo = "github.com/jerjorg/IBZ.jl.git", +# ) diff --git a/docs/src/Documentation.md b/docs/src/Documentation.md new file mode 100644 index 0000000..3c78b14 --- /dev/null +++ b/docs/src/Documentation.md @@ -0,0 +1,53 @@ +# Documentation + +## Index + +### Lattices +```@index +Modules = [IBZ.Lattices] +Order = [:function, :type] +``` + +### Plotting +```@index +Modules = [IBZ.Plotting] +Order = [:function, :type] +``` + +### Symmetry +```@index +Modules = [IBZ.Symmetry] +Order = [:function, :type] +``` + +### Utilites +```@index +Modules = [IBZ.Utilities] +Order = [:function, :type] +``` + +## Functions + +### Lattices +```@autodocs +Modules = [IBZ.Lattices] +Order = [:function, :type] +``` + +### Plotting +```@autodocs +Modules = [IBZ.Plotting] +Order = [:function, :type] +``` + +### Symmetry +```@autodocs +Modules = [IBZ.Symmetry] +Order = [:function, :type] +``` + +### Utilities +```@autodocs +Modules = [IBZ.Utilities] +Order = [:function, :type] +``` diff --git a/docs/src/IBZ.md b/docs/src/IBZ.md new file mode 100644 index 0000000..e69de29 diff --git a/docs/src/Usage.md b/docs/src/Usage.md new file mode 100644 index 0000000..550ab72 --- /dev/null +++ b/docs/src/Usage.md @@ -0,0 +1,38 @@ +# User guide + +`IBZ.jl` calculates the irreducible Brillouin zone (IBZ) of a crystal structure +in 2D and 3D. It also contains functions related to the symmetry of lattices and +lattice reduction. Later on, it will be able to calculate high symmetry points +and paths within the IBZ. + +To calculate the IBZ, simply provide the lattice and atomic basis to `calc_ibz`. +The IBZ will be returned as either a convex hull or intersection of half spaces. +```@example +using IBZ # hide +a = 1.0; +real_latvecs = genlat_CUB(a) +atom_types = [0,0] +atom_pos = Array([0 0; 0.5 0.5]') +ibzformat = "convex hull" +coordinates = "Cartesian" +convention = "ordinary" +ibz = calc_ibz(real_latvecs,atom_types,atom_pos,coordinates,ibzformat, + convention) +``` +The columns of `real_latvecs` are the lattice generating vectors, the columns +of `atom_pos` are the positions of the atoms in Cartesian coordinates (in this +case), `coordinates` are the coordinates of the atom positions, and `convention` +gives the convention for going from real to reciprocal space (whether or not to +multiply by 2π). There is a simple function for visualizing the IBZ along with +the Brillouin zone (BZ). +```@example +using ibz # hide +a = 1.0 # hide +real_latvecs = genlat_CUB(a) # hide +atom_types = [0,0] # hide +atom_pos = Array([0 0; 0.5 0.5]') # hide +ibzformat = "convex hull" # hide +coordinates = "Cartesian" # hide +convention = "ordinary" # hide +plot_cells(real_latvecs,atom_types,atom_pos,coords,convention) +``` diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..cce2ac7 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,9 @@ +# Home + +`IBZ.jl` is a julia package for calculating the irreducible Brillouin zone. It +also contains methods related to the symmetry of lattices and lattice reduction. + +```@contents +Pages = ["Usage.md", "Documentation.md"] +depth = 5 +``` diff --git a/src/Lattices.jl b/src/Lattices.jl index 6c600b0..f2e0e84 100644 --- a/src/Lattices.jl +++ b/src/Lattices.jl @@ -59,7 +59,7 @@ end Calculate the lattice constants and angles of a lattice basis. # Arguments --`latvecs::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. +- `latvecs::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. # Returns - A list where the first element is a list lattice constants `(a,b,c)` and second @@ -103,13 +103,13 @@ closest lattice point from the `k`th lattice vector, and reordering the lattice vectors by increasing Euclidean norms. # Arguments --`basis::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. --`k::Int`: Keeps track of which lattice vector needs to be reduced. +- `basis::AbstractArray{<:Real,2}`: the lattice basis as columns of an array. +- `k::Int`: Keeps track of which lattice vector needs to be reduced. # Returns --`basis::AbstractArray{<:Real,2}`: the partially reduced lattice basis as +- `basis::AbstractArray{<:Real,2}`: the partially reduced lattice basis as columns of an array. --`k::Int`: The index of the lattice vector that needs to be reduced next. +- `k::Int`: The index of the lattice vector that needs to be reduced next. # Examples ```jldoctest diff --git a/src/Symmetry.jl b/src/Symmetry.jl index 6d3958b..f5fc843 100644 --- a/src/Symmetry.jl +++ b/src/Symmetry.jl @@ -212,16 +212,16 @@ end Calculate the point group of lattice in 2D or 3D. # Arguments --`latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns of an +- `latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns of an array. --`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))`: a relative tolerance for +- `rtol::Real=sqrt(eps(float(maximum(real_latvecs))))`: a relative tolerance for floating point comparisons. It is used to compare lengths of vectors and the volumes of primitive cells. --`atol::Real=0.0`: an absolute tolerance for floating point comparisons. It is +- `atol::Real=0.0`: an absolute tolerance for floating point comparisons. It is used to compare lengths of vectors and the volumes of primitive cells. # Returns --`pointgroup::Array{Array{Float64,2},1}`: the point group of the lattice. The +- `pointgroup::Array{Array{Float64,2},1}`: the point group of the lattice. The operators operate on points in Cartesian coordinates. # Examples @@ -285,21 +285,21 @@ end Map a point to the first unit (primitive) cell. # Arguments --`pt::AbstractArray{<:Real,1}`: a point in lattice or Cartesian coordinates. --`latvecs::AbstractArray{<:Real,2}`: the basis vectors of the lattice as columns +- `pt::AbstractArray{<:Real,1}`: a point in lattice or Cartesian coordinates. +- `latvecs::AbstractArray{<:Real,2}`: the basis vectors of the lattice as columns of an array. --`inv_latvecs::AbstractArray{<:Real,2}`: the inverse of the matrix of that +- `inv_latvecs::AbstractArray{<:Real,2}`: the inverse of the matrix of that contains the lattice vectors. --`coords::String`: indicates whether `pt` is in \"Cartesian\" or \"lattice\" +- `coords::String`: indicates whether `pt` is in \"Cartesian\" or \"lattice\" coordinates. --`rtol::Real=sqrt(eps(float(maximum(inv_latvecs))))`: a relative tolerance for +- `rtol::Real=sqrt(eps(float(maximum(inv_latvecs))))`: a relative tolerance for floating point comparisons. Finite precision errors creep up when `pt` is transformed to lattice coordinates because the transformation requires calculating a matrix inverse. The components of the point in lattice coordinates are checked to ensure that values close to 1 are equal to 1. # Returns --`AbstractArray{<:Real,1}`: a translationally equivalent point to `pt` in the +- `AbstractArray{<:Real,1}`: a translationally equivalent point to `pt` in the first unit cell in the same coordinates. #Examples @@ -339,19 +339,19 @@ end Calculate the space group of a crystal structure. # Arguments --`real_latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns +- `real_latvecs::AbstractArray{<:Real,2}`: the basis of the lattice as columns of an array. --`atom_types::AbstractArray{<:Int,1}`: a list of atom types as integers. --`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal +- `atom_types::AbstractArray{<:Int,1}`: a list of atom types as integers. +- `atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal structure as columns of an array. --`coords::String`: indicates the positions of the atoms are in \"lattice\" or +- `coords::String`: indicates the positions of the atoms are in \"lattice\" or \"Cartesian\" coordinates. --`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for +- `rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for floating point comparisons. --`atol::Real=0.0`: an absolute tolerance for floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point comparisons. # Returns --`spacegroup`: the space group of the crystal structure. The first element of +- `spacegroup`: the space group of the crystal structure. The first element of `spacegroup` is a list of fractional translations, and the second element is a list of point operators. The translations are in Cartesian coordinates, and the operators operate on points in Cartesian coordinates. @@ -504,25 +504,25 @@ end Calculate the irreducible Brillouin zone of a crystal structure in 2D or 3D. # Arguments --`real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as +- `real_latvecs::AbstractArray{<:Real,2}`: the basis of a real-space lattice as columns of an array. --`atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. --`atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal +- `atom_types:AbstractArray{<:Int,1}`: a list of atom types as integers. +- `atom_pos::AbstractArray{<:Real,2}`: the positions of atoms in the crystal structure as columns of an array. --`coords::String`: indicates the positions of the atoms are in \"lattice\" or +- `coords::String`: indicates the positions of the atoms are in \"lattice\" or \"Cartesian\" coordinates. --`ibzformat::String`: the format of the irreducible Brillouin zone. Options +- `ibzformat::String`: the format of the irreducible Brillouin zone. Options include \"convex-hull\" and \"half-space\". --`convention::String="ordinary"`: the convention used to go between real and +- `convention::String="ordinary"`: the convention used to go between real and reciprocal space. The two conventions are ordinary (temporal) frequency and angular frequency. The transformation from real to reciprocal space is unitary if the convention is ordinary. --`rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for +- `rtol::Real=sqrt(eps(float(maximum(real_latvecs))))` a relative tolerance for floating point comparisons. --`atol::Real=0.0`: an absolute tolerance for floating point comparisons. +- `atol::Real=0.0`: an absolute tolerance for floating point comparisons. # Returns --`ibz`: the irreducible Brillouin zone as a convex hull or intersection of +- `ibz`: the irreducible Brillouin zone as a convex hull or intersection of half-spaces. # Examples From cee206c2009aa4b0150677ccba830bd0db2b5752 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 20:26:47 -0600 Subject: [PATCH 11/21] Fixed documentation --- .gitignore | 2 ++ .travis.yml | 16 +++------------- docs/make.jl | 2 +- docs/src/Usage.md | 9 ++++----- src/Plotting.jl | 1 + 5 files changed, 11 insertions(+), 19 deletions(-) diff --git a/.gitignore b/.gitignore index 29126e4..0bc7f24 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,5 @@ docs/site/ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml + +notebooks diff --git a/.travis.yml b/.travis.yml index 590e119..237be32 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,26 +11,16 @@ jobs: notifications: email: false script: + name: "Unit tests" - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' -jobs: - include: - - stage: "Documentation" - julia: 1.4 - os: osx - script: # Install required packages for coverage and documentation - - julia --project -e 'import Pkg; Pkg.add("Coverage"); - Pkg.add("Documenter"); Pkg.add("DocumenterTools");' - - julia --project=docs/ -e 'using Pkg; + - julia --project=docs/ -e 'import Pkg; Pkg.add("Coverage"); + Pkg.add("Documenter"); Pkg.add("DocumenterTools"); Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - julia --project=docs/ docs/make.jl - # Submit test coverage report - - julia --project -e 'using Coverage; - Coveralls.submit(Coveralls.process_folder())' - after_success: - julia --project ./docs/make.jl diff --git a/docs/make.jl b/docs/make.jl index cd8b780..bc3ce09 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,7 @@ push!(LOAD_PATH,"../src/") using Documenter, IBZ -makedocs(sitename="IBZ Documentation", +makedocs(sitename="IBZ", format = Documenter.HTML(prettyurls = false), modules = [IBZ], authors = "Jeremy Jorgensen", diff --git a/docs/src/Usage.md b/docs/src/Usage.md index 550ab72..91ecca9 100644 --- a/docs/src/Usage.md +++ b/docs/src/Usage.md @@ -12,7 +12,7 @@ using IBZ # hide a = 1.0; real_latvecs = genlat_CUB(a) atom_types = [0,0] -atom_pos = Array([0 0; 0.5 0.5]') +atom_pos = Array([0 0 0; 0.5 0.5 0.5]') ibzformat = "convex hull" coordinates = "Cartesian" convention = "ordinary" @@ -26,13 +26,12 @@ gives the convention for going from real to reciprocal space (whether or not to multiply by 2π). There is a simple function for visualizing the IBZ along with the Brillouin zone (BZ). ```@example -using ibz # hide +using IBZ # hide a = 1.0 # hide real_latvecs = genlat_CUB(a) # hide atom_types = [0,0] # hide -atom_pos = Array([0 0; 0.5 0.5]') # hide -ibzformat = "convex hull" # hide +atom_pos = Array([0 0 0; 0.5 0.5 0.5]') # hide coordinates = "Cartesian" # hide convention = "ordinary" # hide -plot_cells(real_latvecs,atom_types,atom_pos,coords,convention) +plot_cells(real_latvecs,atom_types,atom_pos,coordinates,convention) ``` diff --git a/src/Plotting.jl b/src/Plotting.jl index 5027f9d..1eac4d8 100644 --- a/src/Plotting.jl +++ b/src/Plotting.jl @@ -373,6 +373,7 @@ convention = "ordinary" """ function plot_cells(real_latvecs,atom_types,atom_pos,coords,convention, rtol::Real=sqrt(eps(float(maximum(real_latvecs)))),atol::Real=0.0) + mplot3d=pyimport("mpl_toolkits.mplot3d") bzformat = "convex hull" bz = calc_bz(real_latvecs,convention,bzformat) ibzformat = "convex hull" From ccd81ee7d7382d36cc354ef8773b7598f71421cd Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 20:30:38 -0600 Subject: [PATCH 12/21] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 455829e..18c5b48 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0" [deps] CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" From 8265ddd77abc206812e464ff0b6e9092277523b5 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 20:40:43 -0600 Subject: [PATCH 13/21] Update .travis.yml --- .travis.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 237be32..1c7715d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,7 +19,8 @@ jobs: - julia --project=docs/ -e 'import Pkg; Pkg.add("Coverage"); Pkg.add("Documenter"); Pkg.add("DocumenterTools"); Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - - julia --project=docs/ docs/make.jl + - julia --project=docs/ -e 'import Pkg; Pkg.add("Documenter"); + docs/make.jl' after_success: - julia --project ./docs/make.jl From 1cc1bc5a6bc36ee73ceda915b528549659d020b6 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 21:03:43 -0600 Subject: [PATCH 14/21] Update .travis.yml --- .travis.yml | 9 --------- 1 file changed, 9 deletions(-) diff --git a/.travis.yml b/.travis.yml index 1c7715d..b81492a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -14,15 +14,6 @@ jobs: name: "Unit tests" - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' - - # Install required packages for coverage and documentation - - julia --project=docs/ -e 'import Pkg; Pkg.add("Coverage"); - Pkg.add("Documenter"); Pkg.add("DocumenterTools"); - Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - - julia --project=docs/ -e 'import Pkg; Pkg.add("Documenter"); - docs/make.jl' - after_success: - julia --project ./docs/make.jl - coveralls: true From 14ad4b43b43774f466c941de772a8fb9b2ed493c Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 21:15:16 -0600 Subject: [PATCH 15/21] Update .travis.yml --- .travis.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index b81492a..23372ad 100644 --- a/.travis.yml +++ b/.travis.yml @@ -13,7 +13,6 @@ jobs: script: name: "Unit tests" - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); - Pkg.activate("."); Pkg.test("IBZ"; coverage=true)' + Pkg.activate("."); Pkg.test("IBZ")' after_success: - julia --project ./docs/make.jl - coveralls: true From 78fe851689099e7a57e502247d35aa8b74cf24a7 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 21:42:39 -0600 Subject: [PATCH 16/21] Update .travis.yml --- .travis.yml | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/.travis.yml b/.travis.yml index 23372ad..751b637 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,18 +1,20 @@ language: julia -jobs: - include: - - stage: "Tests" - install: - - pip3 install scipy - julia: 1.4 - os: - - osx - notifications: - email: false - script: - name: "Unit tests" - - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); - Pkg.activate("."); Pkg.test("IBZ")' - after_success: - - julia --project ./docs/make.jl +install: + - pip3 install scipy + +os: + - osx + +julia: + - 1.4 + +notifications: + email: false + +script: + - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); + Pkg.activate("."); Pkg.test("IBZ")' + +after_success: + - julia --project ./docs/make.jl From 58bed9c6335d077177123d05e0b2903bb9841408 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 22:01:40 -0600 Subject: [PATCH 17/21] Updated README.md --- README.md | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index a1f7a7e..5663f9d 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,5 @@ -# IBZ -A julia package for calculating the Irreducible Brillouin zone (IBZ) for any lattice. - -Utilizes [phenumlib]( https://github.com/wsmorgan/phonon-enumeration) through [PyCall](https://github.com/JuliaPy/PyCall.jl) to get the point group. +[![Build Status](https://travis-ci.org/jerjorg/IBZ.svg?branch=master)](https://travis-ci.org/jerjorg/IBZ) -Phenumlib and PyCall need to be installed. +# IBZ +A julia package for calculating the Irreducible Brillouin zone (IBZ) of a 2D or +3D crystal structure. From 52eeefc1035d7f1ea275f4f75d1ce8a6109ced8e Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 22:05:45 -0600 Subject: [PATCH 18/21] Updated .travis.yml --- .travis.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 751b637..6f82464 100644 --- a/.travis.yml +++ b/.travis.yml @@ -14,7 +14,9 @@ notifications: script: - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); - Pkg.activate("."); Pkg.test("IBZ")' + Pkg.activate("."); Pkg.test("IBZ", coverage=true)' after_success: - julia --project ./docs/make.jl + +coveralls: true From c3991c53fe3a81cda0b49c26bc45a11412b3f656 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 22:29:30 -0600 Subject: [PATCH 19/21] Update README.md --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 5663f9d..5668780 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ [![Build Status](https://travis-ci.org/jerjorg/IBZ.svg?branch=master)](https://travis-ci.org/jerjorg/IBZ) +[![Coverage Status](https://coveralls.io/repos/github/jerjorg/IBZ/badge.svg?branch=master)](https://coveralls.io/github/jerjorg/IBZ?branch=master) + # IBZ A julia package for calculating the Irreducible Brillouin zone (IBZ) of a 2D or 3D crystal structure. From dc7e96593a903a7db0048c5d9f28917f4b5600d6 Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 22:31:53 -0600 Subject: [PATCH 20/21] Update .travis.yml --- .travis.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.travis.yml b/.travis.yml index 6f82464..e8fb10f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -18,5 +18,7 @@ script: after_success: - julia --project ./docs/make.jl + - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; + Coveralls.submit(process_folder())' coveralls: true From 5f27cc2c748a1c1128ce5dd5300d87aca222c42b Mon Sep 17 00:00:00 2001 From: Jeremy Jorgensen Date: Wed, 22 Jul 2020 22:47:44 -0600 Subject: [PATCH 21/21] update README.md --- README.md | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/README.md b/README.md index 5668780..0449935 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,4 @@ -[![Build Status](https://travis-ci.org/jerjorg/IBZ.svg?branch=master)](https://travis-ci.org/jerjorg/IBZ) - -[![Coverage Status](https://coveralls.io/repos/github/jerjorg/IBZ/badge.svg?branch=master)](https://coveralls.io/github/jerjorg/IBZ?branch=master) +[![Build Status](https://travis-ci.org/jerjorg/IBZ.svg?branch=master)](https://travis-ci.org/jerjorg/IBZ) [![Coverage Status](https://coveralls.io/repos/github/jerjorg/IBZ/badge.svg?branch=master)](https://coveralls.io/github/jerjorg/IBZ?branch=master) # IBZ A julia package for calculating the Irreducible Brillouin zone (IBZ) of a 2D or