蜂の巣格子をJulia+PyPlotで描く

#序
初投稿です。
ほぼメモです。

蜂の巣格子は物性物理学において大流行しているグラフェンの結晶構造である。
ヘリカルエッジ状態などトポロジカルな物性が予言されたり，2層ずらして重ねることで超伝導が発現したりして，大変楽しい物質である。

そんな蜂の巣格子の画像を座標の情報から作成し，ベクター形式で保存したいと思った。
(本当はカーボンナノチューブを書きたい。本記事はその準備である。)

#基本並進ベクトルと格子ベクトル
選び方に任意性がある。
ここでは

\vec{a}_1 = \left(\frac{\sqrt{3}}{2}a,\frac{1}{2}a\right) \\
\vec{a}_2 = \left(0,a\right) 

としておく。$a$は格子定数で，原子間距離$d$と$a=\sqrt{3}d$の関係にある。
基本並進ベクトルによって格子ベクトルは

\vec{R} = m\vec{a}_{1} + n\vec{a}_{2} \quad (m,n \in \mathbb{Z})

と表すことができる。
ところで蜂の巣格子は，隣り合うサイトに格子ベクトルによって移動することができない。
したがって格子ベクトルの上には二つのサイトを乗せる。

#環境
Windows 10 Home
julia 1.4.0
vscode
shift + enterで実行可能ということを最近知りました。

#実装
2点間の距離が原子間距離に一致するときに線を引く。
なお，格子を書くだけならInkscape等を使う方がずっと簡単に書くことができる。

HoneycombLattice.jl

using PyPlot
pygui(true)#画像をポップアップで表示

d = 1.4 #隣接原子間の距離
a = sqrt(3) * d #格子定数

#基本並進ベクトル
a1 = [sqrt(3) * a / 2,a / 2] 
a2 = [0.0,a]

#格子ベクトル
function graphane_lattice(m, n)
    return m * a1 + n * a2
end

#適当に座標を作成
N = 12
M = 12
Avec = [graphane_lattice(m, n) for n = -N:N,m = 0:M ] #Aサイトの座標
Bvec = [([d,0] + graphane_lattice(m, n)) for n = -N:N,m = 0:M ] #Bサイトの座標
posvec = vcat(Avec, Bvec) #全サイトの座標
num = length(posvec) #サイトの総数

#格子をプロットする
function LatticePlot()
    fig = figure(figsize=(2, 2))
    ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
    ax.set_xlim(0, 25)
    ax.set_ylim(0, 25)
    ax.set_xticks([])
    ax.set_yticks([])
    
    for i = 1:num
        for j = 1:num
            if i < j #重複なく数える
                r = posvec[i] - posvec[j] #2点間のベクトル
                x = [posvec[i][1],posvec[j][1]]
                y = [posvec[i][2],posvec[j][2]]
                if sqrt(r[1]^2 + r[2]^2) < d + 1e-2 #2点間の距離が隣接原子間距離(+少し)に一致するときにプロット
                    ax.plot(x, y, color="black", lw=1)
                end
            end
        end
    end
    savefig("HoneycombLattice.svg") #保存,pdfやpngでも保存できる
end

LatticePlot() #実行

#参考文献
VSCode で Julia-1.4 を動かすまで

#追記
@tenfu2tea　様からご指摘いただきまして，コードの改善を図りました。

HoneycombLatticev2.jl

using PyPlot
pygui(true)#画像をポップアップで表示

d = 1.4 #隣接原子間の距離
a = sqrt(3) * d #格子定数


a1 = [sqrt(3) * a / 2,a / 2] 
a2 = [0.0,a]

#格子ベクトル
function graphane_lattice(m, n)
    return m * a1 + n * a2
end

#格子をプロットする
function LatticePlot()
    N = 12
    M = 12
    #座標の生成方法を改善
    Avec = hcat((graphane_lattice(m, n) for n = -N:N,m = 0:M )...) #Aサイトの座標
    Bvec = Avec .+ [d,0] #Bサイトの座標
    posvec = hcat(Avec, Bvec) #全サイトの座標
    num = size(posvec, 2) #サイトの総数
    

    fig = figure(figsize=(3, 3))
    ax = fig.add_axes((0.1, 0.1, 0.8, 0.8))
    ax.set_aspect("equal")
    ax.set_xlim(0, 20)
    ax.set_ylim(0, 20)
    ax.set_xticks([])
    ax.set_yticks([])
    
    for i = 1:num
        for j =  i + 1:num
            r = posvec[:,i] - posvec[:,j] #2点間のベクトル
            x = [posvec[1,i],posvec[1,j]]
            y = [posvec[2,i],posvec[2,j]]
            if sqrt(r'r) < d + 1e-2      #2点間の距離が隣接原子間距離(+少し)に一致するときにプロット
                ax.plot(x, y, color="black", lw=1)
            end
            
        end
    end
    savefig("HoneycombLattice.png", dpi=720)
end

LatticePlot() #実行

基本並進ベクトルの選び方であったり，AサイトからBサイトへの移動などまだ改善点はありそうです。