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 new file mode 100644 index 0000000..e8fb10f --- /dev/null +++ b/.travis.yml @@ -0,0 +1,24 @@ +language: julia + +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", coverage=true)' + +after_success: + - julia --project ./docs/make.jl + - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; + Coveralls.submit(process_folder())' + +coveralls: true diff --git a/Project.toml b/Project.toml index f1a975e..18c5b48 100644 --- a/Project.toml +++ b/Project.toml @@ -1,14 +1,21 @@ name = "IBZ" uuid = "49a35663-c880-4242-bebb-1ec8c0fa8046" -authors = ["John Christensen "] +authors = ["Jeremy Jorgensen "] 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" +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" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index a1f7a7e..0449935 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) [![Coverage Status](https://coveralls.io/repos/github/jerjorg/IBZ/badge.svg?branch=master)](https://coveralls.io/github/jerjorg/IBZ?branch=master) -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. diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..bc3ce09 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,12 @@ +push!(LOAD_PATH,"../src/") +using Documenter, IBZ + +makedocs(sitename="IBZ", + 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..91ecca9 --- /dev/null +++ b/docs/src/Usage.md @@ -0,0 +1,37 @@ +# 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; 0.5 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; 0.5 0.5 0.5]') # hide +coordinates = "Cartesian" # hide +convention = "ordinary" # hide +plot_cells(real_latvecs,atom_types,atom_pos,coordinates,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/IBZ.jl b/src/IBZ.jl index e2b76ad..eacf5dc 100644 --- a/src/IBZ.jl +++ b/src/IBZ.jl @@ -1,692 +1,26 @@ 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 - -@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) - -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 - -""" -function reduce_bz(lat, at, atom_pos) - - bz = make_bz(lat,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) - - bz = make_bz_2d(lat,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,vertsOrHrep = true) - # Minkowski reduce the lattice. - vector_utils = pyimport("phenum.vector_utils") - eps=10^-9; - lat = vector_utils._minkowski_reduce_basis(lat',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] - #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,vertsOrHrep = true) - #@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] - #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 refrence - - - -!!! 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) - - bz = make_bz(lat) - 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) - @show polygons - 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 - #sometimes total is 0? TEMP FIX - if total < min_distance #&& total !=0 - @show min_distance - @show total - min_distance = total - min_order = order_num - if total == 0 - @show points - @show order - @show possible_orderings - @show polygon - 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 = [] - 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 - if size(polygon)[1] != 0 - @show polygon - push!(polygons,polygon) - end - 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(.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(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 - - -#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) -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 +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 +# 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..f2e0e84 --- /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..1eac4d8 --- /dev/null +++ b/src/Plotting.jl @@ -0,0 +1,403 @@ +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) + mplot3d=pyimport("mpl_toolkits.mplot3d") + 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..f5fc843 --- /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/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/make_bz.jl b/test/make_bz.jl deleted file mode 100644 index d516faa..0000000 --- a/test/make_bz.jl +++ /dev/null @@ -1,22 +0,0 @@ -using Test,IBZ,PyCall -include("testHelperFunctions.jl") -lats = lattices -labels =lattice_labels - -@testset "make_bz_test" begin - i=1 - @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) - i+=1 - end -end - diff --git a/test/make_ibz.jl b/test/make_ibz.jl deleted file mode 100644 index 47fe7f6..0000000 --- a/test/make_ibz.jl +++ /dev/null @@ -1,41 +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 - 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) - #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) - #rch_old = reduce_bz_old(lat,at,atom_pos) - - @test expectedOrder[i] ≈ size(spaceGroupTranspose)[1] - - @test rch.volume/ch.volume ≈ 1/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 7bd457f..0000000 --- a/test/make_ibz_rand.jl +++ /dev/null @@ -1,43 +0,0 @@ -using Test,IBZ,PyCall -include("testHelperFunctions.jl") -lats = lattices -labels =lattice_labels - -@testset "bz reduce rand" begin - 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) - #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) - #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/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..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 diff --git a/test/testHelperFunctions.jl b/test/testHelperFunctions.jl deleted file mode 100644 index 5b422ab..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