Tips

Lapack

線形計算ライブラリです。詳しいことはintelのマニュアルを参考にしてください。

行列の対角化

zheev関数

エルミート行列の対角化を行う関数です。ここでは単位行列+σ_yを対角化させています。C言語についてこのURLが参考になる

プログラム置き場

六角格子系の状態密度

状態密度をデルタ関数で定義した。六角格子系の状態密度をdensity.txtに出力するプログラム

出力されるdensity.txtをgnuplotでグラフにすると次の図を得る。横軸はエネルギー、縦軸が六角格子系の状態密度である。

Hofstadter's butterfly

六方格子系に垂直に磁場を与えることでHofstadter's butterflyが現れることを確認できます。

出力されるzheev.txtをgnuplotでグラフにすると次の図を得る。横軸が磁場,縦軸が許されるエネルギー準位である。周期境界条件により許される磁場は離散的である。

Energy disperstion of zigzag graphen nanoribbons

グラフェンのジグザグナノリボンのバンド図を求めます。Nはユニットセル内の原子の数です。最近接のとびうつりのみを考えたタイトバインディングで計算しています。計算結果は/dataにdatファイルとして出力します。

N=70の計算結果をgnuplotっでプロットさせたものです。横軸は波数、縦軸はエネルギー(eV)を表しています

メモ

ほぼ備忘録です。いちいち調べるのが面倒なのでメモしているだけです。

Gnuplot

◆軸のメモリの文字の大きさを変える
set tics font "Helvetica,20"
◆y軸のメモリの間隔をかえる
set ytics 0.5
◆凡例を消す.
unset key
◆data plotの色と太さを変える。lwはlinewidth,lcはline color
plot "~.dat" w l lw 5 lc "red"
◆data plotを点線にする。ltはline type
plot "~.dat" w l lt 0
◆余白を左にいれる。右ならrmargin。同様に上下はtとb
set lmargin 8

C言語

◆long double型の精度の桁数を知るためにはfloat.hを読み込んで次のようにする。float,double型の精度を知りたいときはそれぞれFIT_DIG,DBL_DIGにかえる。
printf("%d",LDBL_DIG)
◆コマンドラインからの入力。実行するときには ./a.out 4 12 … などとする。argcには入力した変数の個数-1が代入される。argv[1],argv[2]…には4,12…と入力した値が順番にchar型で代入される。atof()でchar型を実数に直している。 ◆ファイル名に変数を含むファイルを出力したい場合は、ファイル名をchar型の変数で予め定義しておけばよい。

並列計算(Open MP)

◆C言語のプログラム中で,AとBについての演算の結果を最後に足しつつfor文を並列計算させたいときは次のようにする。dynamicは各コアが計算が終わり次第、次の命令を受け取るという意味。ネストしたfor文では最も内側のfor文を並列処理させるとトラブルがすくない。
#pragma omp parallel for reduction(+:A,B) schedule(dynamic)
for(;;){}
test.cをコンパイルするとき
icc -openmp test.c
./a.outをコアを10個つかって実行するとき
env OMP_NUM_THREADS=10 ./a.out

Mathematica

◆G[x]だけ半透明にする方法
Plot[{F[x],G[x]},{x,-1,1},PlotStyle->{Opacity[1],Opacity[0.7]}]
◆行列の固有値・固有ベクトルの取り出し方
行列Hの固有値(q=1)、固有ベクトル(q=2)を取り出す。n番目の固有ベクトルの第m成分を取り出す
EigenSystem[H][[q]]
EigenSystem[H][[2]][[n]][[m]]
◆Eigenvector[A]は規格化されていない固有ベクトルを得る。規格化された固有ベクトルを得るにはNormalize /@が便利。/@はMapの省略記号で、リストの要素それぞれに関数を適用する。
Normalize /@ Eigenvectors[A]
◆正の実数xとyに関する不定方程式x^2+y^2=4の解を1000個求めて、解をx.datとy.datに出力する。
Export["x.dat", List[x/.FindInstance[ x^2 + y^2 == 4 && x>0 && y>0 ,{x, y},Reals,1000]//N],"Data"]
Export["y.dat", List[y/.FindInstance[ x^2 + y^2 == 4 && x>0 && y>0 ,{x, y},Reals,1000]//N],"Data"]
◆3Dplotで枠なし、軸なし、メッシュなし、縦横比1を実現するオプション
Plot3D[{Sqrt[x^2 + y^2], -Sqrt[x^2 + y^2]}, {x, -1, 1}, {y, -1, 1}, Axes -> None, Mesh -> None, AspectRatio -> 1, Boxed -> False]
◆3Dplotでmeshを指定した高度に赤太字でいれる。下の例ではz=6.5,-6.5に入れる。
Plot3D[{Sqrt[2 + x^2 + y^2], -Sqrt[2 + x^2 + y^2]}, {x, -10, 10}, {y, -10, 10}, Mesh -> {{6.5, -6.5}}, MeshStyle -> {Red, Thick}, MeshFunctions -> {#3 &}]
◆指定したcolor dataで等高線を引く。下の例ではColorFunction -> ColorData["LightTemperatureMap"]を変えることでテーマを変えられる。テーマはここで確認できる。
ContourPlot[x^2 + y^2, {x, 0, 3}, {y, 0, 3},ColorFunction -> ColorData["LightTemperatureMap"]]
◆実数x,yに関する多項式の実部だけを取り出す。下の例ではy*Exp[I*x + 2]の実部だけを取り出す。
Assuming[{x \[Element] Reals, y \[Element] Reals}, FullSimplify[Re[y*Exp[I*x + 2]]]]

emacs

◆すべて選択
C-x h
◆選択範囲を削除
C-w
◆選択範囲をコピー
Esc-w
◆コピーしたものを貼り付け
C-y