Juliaで任意の結晶構造をプロット

使用するライブラリ

• LinerAlgebra : 線形代数用
• GLMakie : Plot用

Library.jl

using LinearAlgebra, GLMakie,GeometryBasics

与えるデータ ;
unit cellに含まれる格子点の座標
基本並進ベクトルのデータ (3d)

data.jl

#unit cell内格子点のデータ
o = [0.0,0.0,0.0] #原点は必須
t1 = [0.5,0.0,0.5] #c軸に並行な点は上手くいかない
t2 = [0.2,0.1,0.0]
t3 = [0.3,0.3,0.0]


#基本並進ベクトルのデータ
a1 = [1/2,-√3/2,0.0]
b1 = [1/2,√3/2,0.0]
c1 = [0.0,0.0,1.0]

step1 : 格子点のデータを生成
与えるデータ；
unit cellの数 (xy平面内 $\geq$ 1)
Layer数 (c軸方向に並べる数 $\geq$ 1 )

Lattice.jl

using LinearAlgebra,GLMakie,GeometryBasics



#make data of lattice point
#make_lattice([lattice in unit cell], number of unit cell, [primitive vectors])
function make_lattice(cell,num_cell::Int64,Layer::Int64,primitive)
    vects = []
    dist = []
    

    #格子点プロットの大きさ
    rad = minimum(norm.(primitive))/20
    
    #Lattce data
    Lattice = []

    #線を引かない点のデータ
    rem_num = []
    for i in 0:num_cell
        for j in 0:num_cell
            for k in 0:Layer-1
                vec = i*primitive[1] + j*primitive[2] + k*primitive[3]
                if i == num_cell || j == num_cell || k == Layer-1 
                    push!(Lattice,cell[1] + vec)
                    push!(rem_num,length(Lattice))
                else
                    for n in 1:length(cell)
                        push!(Lattice,cell[n] + vec)
                        if n == 1
                            push!(rem_num,length(Lattice))
                        end
                    end
                end
            end
        end
    end

    #site number
    site = size(Lattice)[1]

    #primitive vectorと等しい距離、かつ、unit cell境界の点を記録
    edge_mat = zeros(Int64,site,site)
    for i in 1:site
        for j in 1:site
            d = Lattice[i] - Lattice[j]
            θ = zeros(length(primitive))
            for N in 1:length(primitive)
                prod = dot(d,primitive[N])/(norm(d)*norm(primitive[N]))
                if abs(prod) <= 1.0
                    θ[N] = acos(prod)
                else
                    θ[N] = 0.0
                end
            end
            
            edge = norm(d)

            cond = minimum((edge .- norm.(primitive)))<= 0.01 && abs(minimum(θ)) <= 0.01
            
            if in(i,rem_num) && in(j,rem_num) && cond
                edge_mat[i,j] = 1
            end
        end
    end

    println(length(findall(edge_mat.!= 0 )))
    println(length(rem_num))

    return Lattice,edge_mat,rad
end

#格子点間に引く線のデータ
function get_Graph_Edges3D(Lattice,edge_mat)
    #格子点間に引く線のデータ
    xyzos = []


    for i in 1:length(Lattice)
        for j in 1:length(Lattice)
            if edge_mat[i, j] != 0
                push!(xyzos, Lattice[i])
                push!(xyzos, Lattice[j])
            end
        end
    end
    return (Point3f.(xyzos))
end


#plot用の関数
function plotGraph3D(cell,num_cell::Int64,Layer::Int64,primitive)
    Lattice,edge_mat,rad = make_lattice(cell,num_cell,Layer,primitive)
    
    segm = get_Graph_Edges3D(Lattice,edge_mat)

    #格子点のデータを meshscatter用に変換
    position_t = []
    xmin = 0.0
    ymin = 0.0
    zmin = 0.0
    xmax = 0.0
    ymax = 0.0
    zmax = 0.0
    for i in 1:length(Lattice)
        push!(position_t,tuple(Lattice[i]...))
        if Lattice[i][1] > xmax
            xmax = Lattice[i][1]
        elseif Lattice[i][1] < xmin
            xmin = Lattice[i][1]
        end
        if Lattice[i][2] > ymax
            ymax = Lattice[i][2]
        elseif Lattice[i][2] < ymin
            ymin = Lattice[i][2]
        end
        if Lattice[i][3] > zmax
            zmax = Lattice[i][3]
        elseif Lattice[i][3] < zmin
            zmin = Lattice[i][3]
        end 
    end
    println("x: ",xmin,"~",xmax)
    println("y: ",ymin,"~",ymax)
    println("z: ",zmin,"~",zmax)
    xmin = xmin - 5rad 
    ymin = ymin - 5rad
    zmin = zmin - 5rad
    xmax = xmax + 5rad
    ymax = ymax + 5rad
    zmax = zmax + 5rad
    position_t=convert(Vector{Tuple{Float64,Float64,Float64}},position_t)

    #plot
    #格子点間の線をプロット
    fig = Figure(resolution=(800, 800), dpi=700)
    ax = Axis3(fig[1,1], aspect=:data, 
        perspectiveness=1.0,
        azimuth = 0.2 * pi,elevation = 0.2 * pi,
        limits = (xmin,xmax,ymin,ymax,zmin,zmax)
        )
    #azimuth, elevation : set view point
    
    
    linesegments!(ax,segm)
    
    
    #格子点上に半径radの球をplot
    meshscatter!(ax,position_t,
            markersize = rad
            )   
    hidedecorations!(ax)
    hidespines!(ax)

    
    
    fig
    
    #pngに出力
    #save("Lattice.png",fig)
end

実行例

data.jl

#unit cell 内の格子点
o = [0.0,0.0,0.0]
t1 = [0.5,0.0,0.75]

#基本並進ベクトルのデータ
a1 = [1/2,-√3/2,0.0]
b1 = [1/2,√3/2,0.0]
c1 = [0.0,0.0,1.5]

#unit cell = 3
#layer = 2 
plotGraph3D([o,t1],3,2,[a1,b1,c1])

結果