ある三次元の配列があるとします。つまり、x,y,zを指定すると一つの値f(x,y,z)があるような状況です。この時、fがある値となるような点の集合を描きたいとします。
これは固体物理学では「フェルミ面を描く」ことに対応しています。Juliaの色々なライブラリを使うことでできることがわかりましたので、記録として残しておきます。
バージョン
Julia 1.8.1
インストール
必要なパッケージは以下の通りです。addでインストールしてください。
- Interpolations
- Meshing
- FileIO
- MeshIO
- GeometryBasics
サンプルデータの生成
サンプルデータとして、
E(k_x,k_y,k_z) = -(\cos(k_x) + \cos(k_y) + \cos(k_z))
とします。この関数がゼロとなるような点の集合を集めると、面になります。固体物理学ではこれはフェルミ面と呼ばれるものです。
この関数は解析的な形で求まっていますが、一般には数値的な値しかないかもしれません。そこで、この関数から適当に生成したサンプルデータを作成します。
using Interpolations
using Meshing
using FileIO #MeshIOもインストールしないといけない
using GeometryBasics
using LinearAlgebra
function make_data(nkx,nky,nkz)
kxs = range(-pi,pi,length=nkx)
kys = range(-pi,pi,length=nky)
kzs = range(-pi,pi,length=nkz)
energy= zeros(Float64,nkx,nky,nkz)
for iz=1:nkz
kz = kzs[iz]
for iy=1:nky
ky = kys[iy]
for ix=1:nkx
kx = kxs[ix]
energy[ix,iy,iz] = -(cos(kx) + cos(ky) + cos(kz))
end
end
end
return energy
end
データの補間
次に、メッシュで切ったデータ点以外の値を計算できるようにするために、スプライン補間のパッケージInterpolationsを用います。
function energy_dispersion(energymesh)
nka,nkb,nkc = size(energymesh)
kas = range(0,2pi,length=nka)
kbs = range(0,2pi,length=nkb)
kcs = range(0,2pi,length=nkc)
interp_cubic = cubic_spline_interpolation((kas, kbs,kcs), energymesh)
function evaluate(ks)
#println(typeof(ks))
vecka = [ks[1],ks[2],ks[3]]
for i=1:3
while vecka[i] < 0
vecka[i] += 2pi
end
while vecka[i] > 2pi
vecka[i] -= 2pi
end
end
#println(vecka)
return interp_cubic(vecka...)
end
return ks -> evaluate(ks),interp_cubic
end
この関数は、データのメッシュをもらうと、関数evaluate(ks)と補間用データを返してくれます。
ですので、
nkx = 21
nky = 21
nkz = 21
energy_mesh = make_data(nkx,nky,nkz)
evaluate,interp_cubic = energy_dispersion(energy_mesh)
のようにすると、データのメッシュを21x21x21で作成した後に、補間用の関数evaluateを作成したことになります。
ポリゴンデータの作成
最後に、ポリゴンを作成します。これは、
function main()
nkx = 21
nky = 21
nkz = 21
energy_mesh = make_data(nkx,nky,nkz)
evaluate,interp_cubic = energy_dispersion(energy_mesh)
μ = 0.0
println(evaluate([0.2,0.1,0.3]))
km = pi
cosmesh = Mesh(evaluate, Rect(Vec(-km,-km,-km),Vec(2km,2km,2km)),
MarchingCubes(iso=-μ), samples=(50,50,50))
save("cos.ply", cosmesh)
a = zeros(Float64,3)
b = zero(a)
c = zero(a)
g = zero(a)
fp = open("fermi.txt","w")
for tri in cosmesh
kcount = 0
for r in tri
kcount += 1
if kcount == 1
a[:] = r[:]
elseif kcount == 2
b[:] = r[:]
else
c[:] = r[:]
end
end
g[:] = (a+b+c)/3 #三角形の中心座標
a[:] = a-c
b[:] = b-c
if norm(a-b) < 1e-10
continue
end
s = sqrt(norm(a)^2*norm(b)^2 - dot(a,b)^2)/2 #三角形の面積
println(fp,g[1],"\t",g[2],"\t",g[3],"\t",s)
end
close(fp)
end
main()
で可能です。
ここで、
km = pi
cosmesh = Mesh(evaluate, Rect(Vec(-km,-km,-km),Vec(2km,2km,2km)),
MarchingCubes(iso=-μ), samples=(50,50,50))
がポリゴンデータを作成しているところです。Rect(Vec(-km,-km,-km),Vec(2km,2km,2km))は、始点が(-km,-km,-km)で、そこからベクトル(2km,2km,2km)を足した点を対角線の点とした時の直方体領域を指定しています。
このデータは
save("cos.ply", cosmesh)
でply形式で保存ができます。ply形式は、MeshLabなどで読み込むことができます。
例えば、
こんな感じです。
また、Makie.jlも対応していますので、綺麗なプロットを描くこともできます。
サンプルコードではデータとして書き出していますから、gnuplotなどでもプロット可能で、この場合には、
みたいにプロットできます。
コードではポリゴンを構成する三角形の面積を求めて、それもファイルに書き出しています。これは、表面積分を行うときの面積素として使うことができますから、表面積の計算や表面での何らかの量の計算が可能となります。